PORE PRESSURE PREDICTION METHODS USING NORMAL COMPACTIONTREND BASED ON SEISMIC INVERSION „

Pore pressure prediction plays an important role in hydraulic fracturing. The pore pressure in shale cannot be directly measured but to be inferred by the normal velocity trend, so methods based on the effective stress theory are dedicated to establishing a function between seismic interval velocity and pore pressure. Among them the mostly used is Eaton’s equation. However, how to precisely quantify the state of compaction remains unsolved. In this study, the AVO/AVA simultaneous inversion was introduced to estimate P-velocity. According to the exponential relationship between the pore pressure and the ratio of normal compaction acoustic transit time to the measured, three methods based on Eaton’s equation were proposed for shale gas pore pressure, respectively. The fitting method avoids the estimation of the normal compaction trend (NCT) by the fitted nonlinear relationship between borehole-side pore pressure (inferred by the Eaton’s equation) and P-velocity. The direct calculation method directly estimates the NCT with the linear trend line (LTL). Whereas the model-based calculation method constructs the NCT model using the NCT transit time well logs obtained by the LTL. The result shows that the approach of estimating the NCT impacts the pore pressure both on the accuracy and the horizontal continuity significantly, which implicates that the constraint between traces must be taken into account when computing the NCT, as well as lithology. theoretical velocities is used for the pore pressure prediction. The theoretical velocity is generally called as normal velocity trend, which indicates the normal compaction trend (NCT). In pore pressure prediction, NCT is used to compute the difference between the measures and the theoretical values, and the pore pressure can be predicted by this difference. Zhang [2011] has reviewed fracture gradient prediction methods and the commonly used empirical methods. Methods based on the effective stress theory [Hubbert and Rubby, 1959] are dedicated to establishing a relationship between seismic interval velocity and pore pressure [Fillippone, 1982]. In fact, other physical parameters are used as well, like resistivity [Eaton, 1972]. So accurate seismic interval velocity is needed for the estimate of pore pressure. On the point of this, there have been some researches for high precision velocity. Hong [2008] used neural network inversion to gain wave impedance, and extracted interval velocity by Gardner formula. Despite this, there remains a need for an efficient method that can precisely quantify NCT. Han et al. [2017] simulated the Pand S-velocities using the clay-plus-silt (CPS) model to gain normal velocity trend, and then completed the pore pressure prediction with Eaton’s equation. The result was good. Getting the proper NCT is a hard work, not only the knowledge to a field is needed for the determination of the coefficient, but also repeated attempts are performed to gain satisfied results. In this case, we tried to find a way to avoid the direct computation of NCT. In this study, AVO/AVA simultaneous inversion was performed to gain accurate seismic interval velocity. Then, the fitting method which avoids the NCT by the fitting formula, the direct calculation method which directly computes the NCT, and the modelbased calculation method which estimates the NCT by constructing the NCT model, were proposed to predict shale gas pore pressure. Lastly, the effect of the NCT generated by different methods on the result can be seen from comparative analysis.


INTRODUCTION
Shale gas usually has the characteristic of high pore pressure [Serebryakov et al., 2002].Pressure prediction plays an important role in hydraulic fracturing, as well as reducing the risks of drilling hazards.It has been studied that the main mechanisms of abnormal high pore pressure are compaction disequilibrium (under-compaction) and hydrocarbon generation [Gutierrez, Braunsdor, and Couzens, 2006].The primary depositional setting for under-compaction is dominantly delta and the lithology is mainly shale [Law and Spencer, 1998].Under-compaction is caused by the rapidly subside sediments in basins, in this case the fluids in pores could only be partially expelled when the formation has extremely low permeability.The remained fluids in pores have to support part of the weight of overlying sediments and so produce abnormal high pressure.Besides, porosity in this case is higher than the normal.But it has been confirmed [Meissner, 1981;Tissot, 1984] that the increase of porosity can be caused by hydrocarbon generation as well, which is more likely to be in gas-bearing shale.Sonic velocities observed by acoustic log is lower than it is thought to be with higher porosity and this change is mainly controlled by the amount of compliant porosity [Sviridov, Mayr and Shapiro, 2017].The difference between the sonic measures and the theoretical velocities is used for the pore pressure prediction.The theoretical velocity is generally called as normal velocity trend, which indicates the normal compaction trend (NCT).
In pore pressure prediction, NCT is used to compute the difference between the measures and the theoretical values, and the pore pressure can be predicted by this difference.Zhang [2011] has reviewed fracture gradient prediction methods and the commonly used empirical methods.Methods based on the effective stress theory [Hubbert and Rubby, 1959] are dedicated to establishing a relationship between seismic interval velocity and pore pressure [Fillippone, 1982].In fact, other physical parameters are used as well, like resistivity [Eaton, 1972].So accurate seismic interval velocity is needed for the estimate of pore pressure.On the point of this, there have been some researches for high precision velocity.Hong [2008] used neural network inversion to gain wave impedance, and extracted interval velocity by Gardner formula.Despite this, there remains a need for an efficient method that can precisely quantify NCT.Han et al. [2017] simulated the P-and S-velocities using the clay-plus-silt (CPS) model to gain normal velocity trend, and then completed the pore pressure prediction with Eaton's equation.The result was good.
Getting the proper NCT is a hard work, not only the knowledge to a field is needed for the determination of the coefficient, but also repeated attempts are performed to gain satisfied results.In this case, we tried to find a way to avoid the direct computation of NCT.In this study, AVO/AVA simultaneous inversion was performed to gain accurate seismic interval velocity.Then, the fitting method which avoids the NCT by the fitting formula, the direct calculation method which directly computes the NCT, and the modelbased calculation method which estimates the NCT by constructing the NCT model, were proposed to predict shale gas pore pressure.Lastly, the effect of the NCT generated by different methods on the result can be seen from comparative analysis.

EATON'S EQUATION
Pore pressure was first predicted by the equivalent depth method [Hottman and Johnson, 1965], which was based on NCT only considered under-compaction and the error would be very large if there existed other mechanisms.Forster [1966] analyzed the feasibility and the existent problems of pore pressure prediction.Eaton [1975] summed up the work of former researchers and proposed Eaton's equation to settle the problem that the predicted pore pressure didn't match the measures well when only under-compaction was considered in geologically complicated basins.The stratigraphic factors, diagenetic processes, stratigraphic assemblage conditions and other genetic mechanisms were took into account.It was actually an empirical method and the empirical coefficient required a lot of logging data to determine.Adapted Eaton's methods depend on the NCT, and pore pressure can either be predicted by electrical resistivity or sonic velocity.When using velocity, NCT can be estimated by the linear trend line (LTL) [Peter, Richard, and Peter, 2004;Tingay et al., 2009]: Where Δt n is normal acoustic transit time at the depth of H, Δt 0 is surface acoustic transit time, and C 1 is compaction coefficient which is determined by several attempts referring to transit time well logs.Pore pressure prediction equation is: (2) Where: P p is pore pressure, Mpa; P 0 is overburden pressure, Mpa; P n is normal compaction pore pressure, Mpa; C 2 is the region exponent relates to overpressure formation mechanisms.In this studied area, a value of 2.0 is taken based on the calibration work with actual pressure data; Δt n is the normal compaction acoustic transit time, μs/m; Δt s is measured acoustic transit time, μs/m.
Borehole-side pore pressure was calculated by Eaton's equation using well logs.The NCT contributes to the pore pressure significantly and the impact can be evidently seen in this process (Figure 1).In Figure 1 (a), the NCT was estimated on Well A and B disregarding the lithology, so the gradient C 1 of the LTL was taken as 0.0001 for Well A and 0.00007 for Well B and were kept the same from the start to end.The pore pressures were presented next to the NCTs.In Figure 1 (b), different lithology was considered.On the acoustic transit time logs, the sudden change implicates the change of lithology.The lithology underneath the sudden change in the depth range of 3284 to 3580m on Well A is dominantly mudstone and shale, while the lithology in the depth range of 3163 to 3283m on Well B is dominantly mudstone.The gradients C 1 were 0.00005 and 0.00016 on Well A for the two parts, and C 1 were 0.00008 and 0.00006 on Well B. The result shows that NCT impacts both value and trend.Compare the pore pressure on Well A in Figure 1 (a) and (b), the maximum value of (b) is lower, and the decrease reaches to 7Mpa.The pressure in (b) shows more details due to the larger difference at the adjacent depth.Besides, the pressure gradient in (b) is smaller than that of (a).It is more reasonable to consider lithology that different lithology gets different background Δt s , so the gradient C 1 needs to be carefully determined to identify the overpressure zone.Overpressure zone is the area where Δt s > Δt n .

THEORY
Prestack inversion techniques use seismic angle stacks which contain rich amplitude information and can effectively improve the ability of lithology and fluid characterization [Li et al., 2007].AVA simultaneous inversion can be divided into two steps: reflectivity with sparse spike inversion and P-, S-wave impedance and density with simultaneous inversion [Yang, 2010;Yuan et al., 2015;Zong et al., 2017].This is performed as follows: Reflectivity is inverted by sparse spike inversion.Sparse spike inversion [Debeye and RIEL, 1990] assumes that large reflectivities are sparse, and they are superimposed on a series of small reflectivities that corresponds to Gaussian distribution.Initial reflectivity is estimated by maximum-likelihood deconvolution.Then the optimal reflectivity can be found by minimizing the objective function 1 [Wang et al., 2006]: (3) Where r is reflectivity, d is recorded seismic data, s is synthetic seismic data, λ is residual weight factor, and p,q are norm factors.The lower λ is, the more sparse r is.A low λ emphasizes the sum of the reflectivity is small and will result in little details.A high λ emphasizes the residual of seismic data is small and will result in more details.
P-, S-wave impedance and density are inverted by simultaneous inversion.Refer to the model-based poststack inversion [Russell and Hampson, 1991], the reflectivity is: (4) Written in matrix form is: ( Where N is the number of the layers, R p is the Pwave reflectivity, L lpi = ln Z pi = ln(V pi ρ i ) is the logarithmic P-wave impedance, and D is the difference matrix.
The low frequency content of P-, S-wave impedance and density well logs are taken the logarithm as the initial solution [Hampson, Russell, and Bankhead, 2005], and using the conjugate gradient method to gain P-, S-wave impedance and density.V p , V s , ρ can be then acquired after taking antilog.The inversion process is presented in Figure 2.

APPLICATION
The studied area is located in Sichuan Basin and the structure spreads to the northeast.Three angle stacks are 7-13°, 17-23°, 27-33°.The survey in Figure 3 displays the horizon.
Incorrect values of the two well-logs were removed first.Well-seismic calibration aims to estimate wavelet which is very important for seismic inversion, the quality of the wavelet directly effects the inversion result [Wang et al., 2015].This had been done in a repeated process.A wavelet was extracted on each well and an average wavelet was calculated by the two wavelets.Three angle stacks were used in AVA simultaneous inversion so three average wavelets were computed.One of the average wavelets is shown in Figure 4 which is represented by the green line.
Well logs' band ranges from 2 to thousand Hz and are usually taken to establish low-frequency model.Low-frequency model is to constrain the inversion and supplement low frequency for seismic data [Zong et al., 2017;Zong et al., 2018].In this inversion, P-, S-velocity and density models were established.The P-velocity model is shown in Figure 5.
Parameters including V p , V s , ρ contrast misfit uncertainties, merge cutoff frequency, and wavelet scale factor were tested to gain the optimal estimation.Then AVA simultaneous inversion could be started.The inversion result is shown in Figure 6.The borehole-side inversion results were extracted and compared with logging curves to see how match they were, this is shown in Figure 7.
According to the comparison, P-velocity inversion result is good, but the density is not accurate.It's hard to gain satisfied density because density is insensitive to the reflectivity.

METHODOLOGY 4.1 THE FITTING METHOD
The fitting method is actually a rock physics method.It first uses the Eaton's equation to calculate the boreholeside pore pressure.Then borehole-side pore pressure and the well log P-velocity are fitted to gain a formula, and pore pressure of the adjacent area around the well can be estimated by it.Since only two parameters that pore pressure and P-velocity are contained in the formula, the fitting method can avoid the calculation of regional NCT.

THE DIRECT CALCULATION METHOD
The direct calculation method directly uses the Eaton's equation to calculate regional pore pressure in the target layer.The LTL Formula (1) is directly used to compute the normal compaction trend.Before calculating, time-to-depth domain conversion needs to be performed to gain each point's depth.The applicable way is to use the average velocity.

THE MODEL BASED CALCULATION METHOD
As the aforementioned conclusion says, the NCT contributes to the pore pressure significantly.But the fitting method only computes the borehole-side NCT and omits the region's, which is a bit of careless.Besides, the horizontal continuity of the direct calculation method is poor because the NCT is computed trace by trace.There must exist constraint between traces for a set of sediments of the stratum.Ignoring this kind of constraint incurs the discontinuity, this is presented in Figure 11(a).The model-based calculation method considers the constraint.Well logs are used to construct the NCT model with an interpolation method, and then the pore pressure will be calculated by Eaton's equation.This method is applicable for the area with several wells.

THE FITTING METHOD
Borehole-side pore pressure was calculated with the consideration of lithology, as is shown in Figure 1 (b).The pore pressure versus P-velocity crossplot and its fit are shown in Figure 8, and the prediction result is shown in Figure 9.
The fitted formula is: The fitting method only needs to estimate the borehole-side NCT, and the pore pressure is calculated by the simple fitted formula, which contributes to the high efficiency.

THE DIRECT CALCULATION METHOD
The regional NCT was calculated by Formula (1), and the compaction coefficient C 1 was the average of the two wells' in the depth of the target layer, which is 0.00011.Then the pore pressure was calculated by Formula (2).The result is shown in Figure 10 (a), and the pore pressure coefficient is shown in Figure 10 (b).Pore pressure coefficient is the ratio of pore pressure to hydrostatic pressure: (10) Where P c is pore pressure coefficient, P p is pore pressure, and P n is hydrostatic pressure (i.e.normal compaction pore pressure).P c > 1 indicates the overpressure, whereas P c < 1 indicates the abnormal low LEI ET AL.  pressure.It is obvious that horizontal continuity is very poor.The direct calculation method calculates each trace's NCT independently and without constraint between traces, which causes horizontal discontinuity.

THE MODEL BASED CALCULATION METHOD
The NCT model was established with the two wells' NCT (shown in Figure 1 (b)).Traces between the two wells were interpolated by the Inverse Distance method.Figure 11 shows the comparison of the NCT in the direct calculation method and the established model.It can be seen that the established model has a very good horizontal continuity.The pore pressure and pore pressure coefficient obtained by the model-based calculation method are shown in Figure 12.The result shows an excellent horizon-tal continuity.
Predicted pore pressures beside the borehole of the fitting method, the direct calculation method and the model-based calculation method were exacted as pseudo logs shown in Figure 13.The lest figure is the comparison of pore pressure, and the right figure is the comparison of each method's relative error.The blue medium line is the inferred borehole-side pore pressure as is shown in Figure 1 (b) (pore pressure can't be directly measured but must be inferred by quantifying the state of compaction), the black medium line is the result of the fitting method, the green small dashed line is of the direct calculation method, and the red large line is of the model-based calculation method.
The pore pressure of the model-based calculation method is the most accurate at the Well A, but is the lowest at the Well B. The pore pressure of the fitting method is neither good nor bad.However, the pore pressure of the direct calculation method is not good enough to be put into product.Although there is a better match than the other two methods on the Well B at the time range of 1.15 to 1.2s.The direct calculation method directly uses the LTL to calculate the NCT, it is the easiest access to pore pressure but the precision is too low to be used.This implicates that the direct calculation of the NCT is not desirable.The fitting method is the fastest but the result is not good enough.The model-based calculation method is at both the good and the bad extremes.So in production, both the fitting method and the model-based calculation method should be performed and guide the hydraulic fracturing based on the integrated consideration.

CONCLUSION
Pore pressure in shale cannot be directly measured but to be inferred by the NCT.There are two elements of predicting pore pressure using Eaton's equation.One is the accurate seismic interval velocity, another is to quantify the NCT.Gain accurate seismic interval velocity has been settled, but how to precisely quantify the NCT remains unsolved.
The NCT contributes to the pore pressure a lot.When estimating the NCT, lithology must be considered to gain more reasonable and precise result.In this study, three different methods are proposed to predict pore pressure, the difference between them is how to quantify the NCT.The fitting method only bothers to calculate the borehole-side NCT and it can be regarded that the borehole-side NCT is took to rep- resent the region's.This method has the highest efficiency, and the accuracy is at a medium level.The direct calculation method is the easiest.This method directly calculates the regional NCT by the LTL, but the precision of the pore pressure is too low to be used, which implicates that the direct calculation of the NCT is not desirable.The model-based calculation method estimates the regional NCT by constructing the NCT model, and the accuracy of the result is at two extremes.It is suggested that both the fitting method and the model-based calculation method be integratedly used in practice.
Besides, the result also shows that the different approach of quantifying the NCT impacts the horizontal continuity significantly.The direct calculation method calculates each trace's NCT independently, and the result shows awful horizontal continuity.Whereas, the model-based calculation method considers the constraint with the interpolation method and gets a better estimation.When computing the NCT, the constraint between traces must be taken into account.

FIGURE 3 .
FIGURE 3. Survey of the studied area.The color shows the horizon.

FIGURE 4 .
FIGURE 4. One of the average wavelets coded by green.

FIGURE 5 .
FIGURE 5. P-velocity model.Low-frequency model is to constrain inversion and supplement low frequency for seismic data.

FIGURE 7 .FIGURE 9 .
FIGURE 7. (a) Well A. (b) Well B. Comparison of well logs and the inversion result.Red dashed lines are the inversion results, and blue lines are the well logs.P-velocity matches the log well, but density is not so good.It's hard to gain satisfied density because density is insensitive to the reflectivity.

FIGURE 8 .
FIGURE 8.The fit of the borehole-side pore pressure and P-velocity.

8
FIGURE 10.(a) The direct calculation method.Pore pressure in the target layer.Horizontal continuity is very poor, which is caused by independent computes of each trace's NCT.(b) The direct calculation method.Pore pressure coefficient in the target layer.Pressure prediction results of the direct calculation method.The Horizontal continuity is poor.

FIGURE 11 .
FIGURE 11.(a) The NCT of the direct calculation method.The horizontal continuity is poor, since it ignores the constraint between traces.(b) The Established NCT model.Traces between the two wells were interpolated by the Inverse Distance method.The horizontal continuity is good.Comparison of the NCT.Horizontal continuity of the direct calculation method is poor.Whereas the NCT model enhances the continuity.

FIGURE 12 .
FIGURE 12. (a) The model-based calculation method.Pore pressure in the target layer.The horizontal continuity has a great improvement.(b) The model-based calculation method.Pore pressure coefficient in the target layer.Pressure prediction results of the model-based calculation method.The horizontal continuity has a great improvement.

FIGURE 13 .
FIGURE 13.(a) Comparison of well A's result.The left is the comparison of pore pressure, and the right is the comparison of each method's relative error.The relative error of the model-based calculation method is the smallest, which is about 11%.(b) Comparison of well B's result.The left is the comparison of pore pressure, and the right is the comparison of each method's relative error.The relative error of the fitting method is the smallest, which is about 10%.Comparison of the borehole-side result.Blue medium line is the inferred borehole-side pore pressure, the black medium line is the result of the fitting method, the green small dashed line is of the direct calculation method, and the red large line is of the model-based calculation method.