Pore pressure prediction using seismic velocity modeling: case study, Sefid-Zakhor gas field in Southern Iran

Pore pressure estimation is important for both exploration and drilling projects. During the exploration phase, a prediction of pore pressure can be used to evaluate exploration risk factors including the migration of formation fluids and seal integrity. To optimize drilling decisions and well planning in abnormal pressured areas, it is essential to carry out pore pressure predictions before drilling. Mud weight and fracture gradient are essential parameters to have wellbore stability, prevent blowout, lost circulation, kick, sand production and reservoir damages. Predrill pore pressure accurate prediction allows the appropriate mud weight to be selected and allows the casing program to be optimized, thus enabling safety by preventing wellbore collapse and economic drilling by reducing the cost. The goal of this study is to estimate pore pressure relation with subsurface velocity in the Sefid-Zakhor gas field. Manufactured sonic logs are modified using the check shot interval velocity of Sefid-Zakhor well No. 1. The final acoustic impedance model is converted to the velocity model by removing density. Finally, the velocity model is converted to pore pressure using Bowers (in: IADC/SPE drilling conference proceedings, 1995) relation. The results of the pore pressure model are validated by pore pressure data obtained by the MDT well test tool. Generally, the results show the normal trend for pore pressure in the area, except in the left side of the anticline in the 2D seismic section, because of tectonic uplifting.


Introduction
The methods used for obtaining the interval velocity model for an area depend on the complexity of subsurface structure in the area (Mavko et al. 2009). Sometimes maybe a simple interval velocity model has good results and matches with well data. Pore pressure estimation is important for both exploration and drilling projects (Jamali and Sokooti 2008). During the exploration phase, a prediction of pore pressure can be used to evaluate exploration risk factors including the migration of formation fluids and seal integrity (Hearn and Meulenbroek 2011;Swarbrick et al. 2005). To have optimizing decisions for drilling decisions and well planning in abnormal pressured areas, it is essential to carry out pore pressure predictions before drilling (Hosseini et al. 2019). Mud weight and fracture gradient are essential parameters to have wellbore stability, prevent blowout, lost circulation, kick, sand production and reservoir damages (Baofu et al. 2008;Taillandier et al. 2008). Predrill pore pressure accurate prediction allows the appropriate mud weight to be selected and allows the casing program to be optimized, thus enabling safe by preventing wellbore collapse and economic drilling by reducing the cost (Cibin et al. 2008;Coevering et al. 2007). However, in this research modified the seismic interval velocity model is used for pore pressure prediction, but the correct velocity model is the most important seismic parameter (Zhu et al. 2008;Sheehan et al. 2005). Interval velocity is an important parameter in determining lithology, depth migration, reservoir characterization and even variation of fluid type in the oil reservoir. Interpretations that need correct depth such as fault are more accurate and easier where correct distributions of velocity in the area are available. Buitrago et al. (2010) employed REVEL™ residual velocity analysis to refine the input gathers and velocity field for pressure prediction in Cretaceous carbonates and further processed to produce an inverted velocity cube. Comparison between Eaton's method, originally developed and applicable for the Gulf of Mexico, shows that this method cannot be applied in the present carbonate-dominated environment, while Bowers method seems more suitable to predict pore pressure. Yu (2010) introduced two main approaches for estimating pore pressure: geological using basing modeling and geophysical using seismic velocity. SNR and resolution can significantly impact velocity analysis and the quality of pore pressure prediction. In his article, the impact of the input data quality in terms of SNR and frequency bandwidth on the accuracy of seismic velocity analysis, and ultimately, the reliability of pore pressure estimation are reviewed. Etminan estimates the pore pressure of the Sefid-Zakhor gas field by converting velocity to pore pressure using Bowers equation. He obtained the velocity model using model-based inversion and well No. 1 of the Sefid-Zakhor field. Bowers constant obtained in his research is different from our study. Babu and Sircar (2011) predicted pore pressure from seismic velocities at the synclinal and flank part of the Atharamura anticline, and overpressured zones are identified. Pore pressures were predicted using both modified equivalent depth and modified Eaton's methods to select the best prediction method. When the predicted pore pressures were compared with offset well pore pressures, an excellent match is observed with the offset wells measured pore pressures. Sompotan et al. (2011) estimated pore pressure cube for a southwestern oil field of Iran. Firstly he creates a velocity cube of the field. This velocity cube is created by using the combination of interval velocity cub, obtained by inversion, for reservoir depth and interval velocity obtained by converting stacking velocity to interval velocity-by Dix equation for depth above the reservoir region to surface. He used the geostatistical method for the spatial distribution of interval velocity in depth from the surface to the reservoir. Finally, he obtained a pore pressure cube using Bowers pressure-velocity relation. Li et al. (2012) presented a paper about pore pressure and wellbore stability. In this paper, commonly used methods, for analyses of pore pressure, in situ stress and borehole shear failure, are evaluated for their strengths and weaknesses. Examples are provided that demonstrate that integration of predrill pore pressure and geomechanics analyses with real-time monitoring consistently provides an effective way to mitigate predrill uncertainties and improve well construction efficiency. Etminan (2010) used a conventional method (stacking velocity) and acoustic impedance inversion method in the reservoir region to obtain the velocity model for the area. But in this study velocity model is obtained based on velocity updating using refraction tomography in near surface. Also, some manufactured sonic logs using refraction results and converted stacking velocity to interval velocity are created. Using these manufactured sonic logs and acoustic impedance inversion, a high-resolution velocity model is obtained for the area. In summary, the difference between this study with the previous study is the procedures of generating the velocity model and the resolution of the final velocity model of the area. In this study, the methods of prediction and estimation of pore pressure are initially investigated. Then, pore pressure is calculated for part of the Sefid-Zakhor reservoir using seismic data. Seismic data are converted to the impedance model using Hampson-Russell software (2006), and the velocity model is obtained by removing the density effect. By finding the parameters of the Bowers method for the region, the velocity model is converted to the pore pressure model. The Sefid-Zakhor structure is an anticline with 40 km length and 8 km width at the surface, where Bakhtiary, Aghajari, Mishan, Gachsaran and Asmari formations have outcrop. It is surrounded by "Halegan" and "Lar" anticline which are elongated NW-SE. Sefid-Zakhor is an asymmetry structure with two crests (Fig. 1).

Pore pressure estimation methods
Generally, there are two methods for pore pressure estimation; direct and indirect methods. Direct methods are such as well test which measures pressure directly. Indirect methods are those which measure pore pressure from physical-geophysical properties deflection of formation with respect to their normal state. This deflection would be calibrated with pressure variation, and then pressure is calculated (Chopra and Huffman 2006).

Overpressure detection by well logs
Pore pressure variation changes the physical properties of rock such as compaction, porosity and electrical properties. These changes observed as anomalies in the overpressure zone rather than their normal state. Almost overpressure detection by well logs is a qualitative method.

Pressure effect on sonic log
The sonic log is sensitive to pressure changes since there is a relation between porosity, compaction, acoustic parameters and pore pressure (Chopra and Huffman 2006;Li et al. 2000). An increase in burial depth leads to an increase in compaction and a decrease in porosity (Zhu et al. 1992;Keary et al. 2002). Hence, shear and bulk modulus increases with depth leading to a decrease in rock compressibility and an increase in velocity. In overpressure zone, the compaction rate is decreased, and sometimes it even ceases. The acoustic velocity decreases in the overpressure zone because of the decrease in compaction. Hence, the sonic log is a powerful tool to detect overpressure zones. The sonic log follows the same physics' principles as a seismic survey, but the difference is a range of frequency. To predict pore pressure using the sonic log, a number of relations are developed. The most basic and reliance relations are Eaton (1969Eaton ( , 1972 and Bowers (1995). These are used for both sonic and seismic data to convert velocity to pore pressure, especially for sandstone formations. Both relations have several constants that are determined by well test.

Pressure effects on density log
An increase in compaction leads to an increase in bulk density since pore volume decreases with respect to bulk volume (Dutta and Khazanedari 2006;Kelly et al. 2005). The increase in bulk density with depth is steady. Several theories are represented for increasing density by depth. In the overpressure zone by considering naught decreasing porosity and naught increasing rock compaction, density increasing trend with depth would be slow down or stop. But in overpressure zone rock density is smaller normal state due to a decrease in compaction. The density log by companionship with seismic is used to predict pore pressure quantitatively.
The density log is used to calculate the overburden pressure (Eq. 1).
where g m s 2 . is acceleration gravity, and b kg m 3 . is bulk density calculated by (Dutta 2002a, b); . is porosity, f . is fluid density, and g is grain density. In the absence of a density log, the sediment density can be estimated from the depth below the seafloor using Traugott relations such as (Eq. 3): But having the density data within the study area increases the accuracy of calculating the reservoir overburden pressure, and the results of calculating the conversion of velocity to pore pressure will be reliable. In unconsolidated and overpressure zone the slope of both density and the sonic log would be negative, and the trend of the curve would be slow down.

Pressure effect on resistivity log
Increasing pore pressure causes decreasing grain contacts and not decreasing in pore space, so electrical conductivity in rock is depended on fluid. Resistivity changes are extremely dependent on fluid type and lithology, so resistivity log should be used cautiously as a detecting index for the (Abtahi 2008) overpressure zone. Eaton relation (1975) is used to convert resistivity to pore pressure (Eq. 4): where P pore . is estimated pore pressure, P ovb . is overburden pressure isl pore pressure, R a . is observed resistivity, and R n is the normal resistivity.

Pressure effects on porosity logs
Effective pressure on rock is one of the most important controlling factors of rock porosity.
An increase in effective pressure with depth decreases porosity. This trend is reduced or ceased in the overpressure zone. Basic methods to estimate pressure from porosity are called porosity-based method (Swarbrick 2002). In these methods, it is assumed that there is a direct relation between measured porosity and pressure. Terzaghi's relationship (Chopra and Huffman 2006) represents the effect of effective pressure on porosity by depth (Eq. 5).
The complexity of petrology and lithology introduces inaccuracies in basic porosity methods since the pressure is not the only controlling factor of porosity. Physical and chemical factors also affect porosity (Swarbrick 2002).

Advantage and disadvantage of well logs in pressure estimation
Logs are usually used for the qualitative detection of overpressure zone since their quantitative pressure predictions are generally inaccurate. Also, the developed relationships that are used for pressure estimations are empirical. This estimation is only used in well location, not between the wells; however, they have appropriate accuracy. Presently, methods such as neural networks are used to estimate pressure from well log quantitatively. These methods are used after the well has been drilled or logged. The high accuracy is only at the well location, and it is possible to have uncertainty between the wells because of the lack of information removal between the wells. Pressure estimation from well logs has proper accuracy, since high precision and log accurate variation are based on wellbore properties. In addition to all drilled well logging is done. The disadvantage is that pressure is estimated after drilling the well, and the scope is limited to the area around the well. Also, pressure interpolation between the well errors is possible, since drilling mud penetrates into the formation making an invasion area; the wells logs are strongly influenced by this invasion area.

VSP and pressure estimation
In recent years seismic methods within the well have various applications, such as surface seismic calibration and detecting overpressure zone below the bit. Nowadays new well equipment can obtain three-component compression and shear waves, so complete information of the formation is obtained. The transient time for the seismic wave is lower than the surface seismic survey. The signal-to-noise ratio is high since in well seismic method receivers are more near to reflector tops. Vertical seismic profile (VSP) resolution is higher than surface seismic data. VSP is also used for deeper seismic data where quality is poor. But zero-offset VSP data have small spread around the well, and usually, its data are below the bit (Blackburn et al. 2007;Badri et al. 2005). Seismic inversion, seismic wave velocity calculation and velocity anomaly detection help to detect the overpressure zone. Also, VSP is almost used as a qualitative tool to determine pressure below the bit and rarely used to estimate pressure quantitatively. VSP operations need drilling to be stopped completely, but the frequency bandwidth of VSP gives a high resolution of the overpressure zone.

Well tests
Pressure well test is the most accurate method to measure the pressure within and near the well. Like previous method, pressure well test is done after drilling the well and has little applicability to prevent drilling risk. DST and RFT are two common tests of this method. Both RFT and DST results are used to estimate constant values in related equations and estimated pressure validation.

Gravity data in pore pressure estimation
Among the non-seismic methods, gravity methods have particular importance for use in estimating the pore pressure, because shale has lower bulk density. In high-pressure shale, the density varies about 2100 to 2300 kg/m 3 , so the pressure anomalies could have prospected in the large mass of shale using this approach. However, gravity method has a good potential in prospecting abnormal pressure, but the method is inherently low resolution. In most marine areas, the high-pressure transition region that can be thousands of feet, gravity cannot be predicted. Musgrave and Hicks in 1966 to separate under pressure shale and salt used reflection operation and gravity method.

Pore pressure estimation using seismic survey
Pore pressure could be estimated by seismic data. The basis of this method is based on the premise that any physical change in formation effects on seismic wave properties. Seismic survey is usually done before the main drilling phase of the field, so this method is the only way whose results can be used in drilling phase (Badri et al. 2000). Porosity and rock compaction are two parameters that control the subsurface seismic responses. These parameters are influenced by the effective stress. Increasing the effective stress decreases the porosity. Decreasing the effective stress increases the pore pressure and porosity, so grain contact decreases and the result is decreasing the velocity wave passing through the rock. The differences between seismic and sonic waves are their frequency and pathway through the rock. One of the challenges in performing geopressure prediction with seismic velocities is that the seismic and wellbore velocities often do not calibrate properly against each other. When this occurs, an obvious question arises regarding which data type provides the best base calibration for geopressure. The sonic velocities often do not provide the best calibration for seismic-based prediction because of differences in the frequency of measurement compared to seismic data and due to invasion and other deleterious wellbore effects. Velocity, quality factor, energy absorption and instantaneous velocity are seismic factors that are sensitive to pressure changes. The studies about pressure effect on seismic attribute are low, except velocity and quality factor. In fact, velocity and effective stress changes are the same concept for elasticity modulus of the rock. In geomechanical prospective, pressure changes cause change in the rock modulus. Two practical methods which are used for pore pressure estimation from seismic velocity are Eaton (1972) and Bowers (1995) methods.

Eaton relation
The relation was presented in 1969 by Eaton, and he modified it in 1972 (Eq. 6): where P p . is estimated pore pressure, P ob . is overburden pressure, P hyd . is hydrostatic pressure at the depth of investigation, and k is the constant of Eaton relation which almost is assumed 3 in all sandstone reservoirs. V n is calculated from the sonic log in shall zone with normal consolidation. V i also is obtained from seismic data and Dix formula. In fact this relation used the comparison between petrophysical properties of normal compression region and investigated area for pore pressure. The major disadvantage of this method is that only the non-consolidated shally layer is the overpressure agent, and also normal consolidated shale profile (depth-velocity curve) in the region is necessary which may be not available everywhere, especially in carbonate reservoir.

Bowers equation
Bowers relation is also called an effective pressure method. In 1994 Bowers presented the relationship between pressure and velocity. He considered normal consolidation, abnormal consolidation and unloading affects. In addition unlike the Eaton relation, it does not need normally consolidated shall profile (Bowers 2002).
where V 0 is velocity in unconsolidated sediments saturated with fluid, and A, B are controlling coefficients for velocity variation with effective pressure increasing (Bowers 2002).
In unloading region Eq. 7 should be corrected. For this purpose, maximum stress in which the unloading process is started should be known and the unloading path should be determined by an exponential function from adjacent wells. Incorporating unloading constants into Eq. 6 gives Eq. 7.
where u is unloading constant and P efc max is maximum unloading starting stress: where V max is maximum observed velocity. Given the Bowers characteristics of bowers' relationship, the probability that this relationship creates an appropriate response in the carbonate reservoir is high. If in some areas, especially carbonate reservoir any above-mentioned relation does not answer, a pressure-velocity cross-plot could help.

Case study
In this study 2D seismic data of the Sefid-Zakhor gas field are studied. The base of this study is obtaining accurate velocity data from seismic and converting it to pressure. Using the Hampson-Russell software (2006), the region's acoustic impedance model was obtained and by removing the density effect of the acoustic impedance model, the region velocity model was created. Finally, the velocity model is converted to pore pressure using Bowers' relation. Programming was carried out on MATLAB software to convert the velocity model to pressure.

Methodology
Using check shot data from Sefid-Zakhor well No. 1 that is located 100 m away from the 2D seismic line, sonic and density logs were corrected in terms of depth. Then, synthetic traces were calculated using the sonic and density logs and a Ricker wavelet. Initially, we tried to find the best Ricker wavelet to find the phase type of the wavelet. After trying several linear and minimum phase wavelets with different wavelet parameters such as wavelet length, sample rate, dominant frequency and the best Ricker wavelet was minimum phase wavelet with 25HZ and 125 ms of dominant frequency and wavelength, respectively. In this study, because the quality of the seismic section is poor as shown in Fig. 2, correlation greater than 50% was considered acceptable. Figure 3 shows the wavelet in frequency and time domain, cross-correlation between synthetic trace and seismic traces and corrected P-wave. Is Use Well Full Wavelet Extraction a module in Hampson-Russell.
One of the most important steps for acoustic impedance inversion is to build an initial model. This model is built using sonic logs. The initial acoustic impedance model is created using the data of well No. 1 and horizons of Dashtak, Kangan, upper Dalan and Nar (Fig. 4). The kriging method is used to control the manner of interpolation. This interpolation technique uses a simple linear variogram and then calculates interpolated points by kriging between the wells. Before inverting for P-impedance, an analysis is carried out to find the best parameters which provide minimum error for the model base inversion method.
Final inversion for the post-stack 2D seismic section is computed by the mentioned inversion parameters. This inversion is also done on the initial model which is made by only Sefid-Zakhor well No. 1 to compare the analysis result by an initial model made by manufactured logs and Sefid-Zakhor well No. 1. The result of the final inversion is given in Table 1. In well No. 1-Sefid-Zakhor well No. 1-the RMS error of full well (real and manufactured) initial model is 3051, and in an initial model created by only well No. 1, the RMS error is 5364.
Synthetic seismogram and seismic data correlation with the errors are shown in Fig. 5. The total correlation between the actual trace and synthetic trace used for inversion is 99%, and the RMS error is 0.05%.
The parameters of the inversion analysis were applied to the whole seismic section, and the final acoustic impedance model was obtained. The density effect was removed from the impedance to generate a velocity model. Effective pressure generated by the velocity model is shown in Fig. 6. As it is obvious the effective pressure layers are different from lithology layers. Also, effective pressure is increasing with depth; however, an effective pressure anomaly of overpressure zone is observed on the left edge.
Based on Terzaghi's rule the effective pressure is the difference between overburden pressure and pore pressure (Terzaghi 1943). Therefore, the overburden pressure and effective pressure are required to obtain pore pressure. The overburden pressure is calculated by Traugott (1997) equation and applied to the whole area (Eqs. 3 and 4). The pore pressure of the reservoir region is measured by MDT (Modular Formation Dynamic Tester) tools in 18  Fig. 7. The correlation coefficient is 78%. This relation was applied to the whole of the area. Finally, formation pressure for these points is obtained using Eq. 10:

Results and discussion
The relation between effective pressure and interval velocity using Bowers' equation for Sefid-Zakhor gas field is (Eq. 11): This procedure applied to the inverted velocity model using the MATLAB package. The converted effective pressure model is similar to the velocity model. The converted overburden pressure model increases with depth as expected (Fig. 7). The overburden pressure of the reservoir from the "Kangan" formation to total depth TD is in the range of 1.4 × 10 4 -1.8 × 10 4 psi. Effective pressure generated by the velocity model is shown in Fig. 6. As it is obvious, the effective pressure layers are different from lithology layers.
(10) P p = P ob − P eff P eff = 143.86V 0.461 int (11) R 2 = 0.7822  Synthetic seismogram and seismic data correlation with the errors. The red traces are synthetic, and the black one is the seismic trace. The red lines on the right-hand side of the seismic trace show the error between the synthetic used for inversion and real seismic data. The red line is the Zp across the well obtained from inversion, the black line is Zp from the initial guess model, and the blue line is Zp from the logs Also, effective pressure is increasing with depth; however, there is an effective pressure anomaly in the left edge of the anticline.
The pore pressure model is generated by the difference between overburden pressure and effective pressure (Fig. 8). An overpressure zone is observed on the left side of the anticline edge. Figure 9b shows the comparison between the estimated pore pressure and well test pore pressure data and its verification at the well location. The maximum drift is − 52 psi at depth of 4091 m. The red points are estimated pore pressure, the blue points are MDT test data, and greens are error-corrected estimated pore pressure. As it is described in the previous chapter the drift pore pressure between the estimated pore pressure and MDT pressure is generated as a function of velocity.

Conclusions
The important conclusions of this research are: 1. The accuracy of the initial velocity model affects the accuracy of the final refraction velocity model using the SIRT method. 2. Interpret and division of time-location data of refraction layers in detail with more layers cause lateral velocity variation to be seen better in the final refraction tomography velocity model. 3. Since pore pressure causes lateral velocity variations, more sonic logs for creating the initial acoustic model cause more accuracy of the final velocity model. 4. An accurate velocity model increases accuracy in pore pressure estimation using the seismic velocity model. 5. Reliable manufactured sonic logs can be used instead of real sonic logs to see the variation of velocity and pres-sure changes in the area where the appropriate number of real sonic logs is absent. 6. To assess whether manufactured wells have positive or negative results on the final model, post-stack inversion analysis should be performed. If the correlation between the synthetic seismogram is increased, RMSE is implying a positive effect and vice versa. 7. For the studied area general Bowers equation could give an appropriate conversion between velocity and effective pressure; because of the estimated pore pressure and MDT tool, pore pressure data had a good match in the reservoir region. 8. The effective pressure and pore pressure of the area are in the normal trend except on the left edge of the anticline.