Pore pressure prediction based on rock stress applied to Marmousi seismic data

Abstract

This paper focuses on the modeling of a sedimentary basin where the exploration of oil and gas takes place. In the modeling, we calculate the vertical variation of pore and rock pressure, which serve as natural pumps for the accumulation of fluids, where we use post-migration surface seismic data and information as necessary. We compare two methods of pore pressure calculation. In the first case, we calculate stress due to vertical loading created by gravitational geological overburden. In the second case, we propose a more complex method to calculate the stress distribution based on the mechanics of solids under gravity loading, using the concept of the first tensor invariant and the linear behavior of Hooke’s law. We prove the use of the P and S velocities and density information to calculate rock, pore, and effective pressure distribution, useful for characterizing potential zones for oil and gas accumulation. The proposed method allows formulation of rock pressure from different calculations instead of as a simple overburden pressure value. The joint analysis of the calculated sections can be used to identify patterns and correlations, outline and characterize the target zones, and make practical geological conclusions.

KEYWORDS:
pore pressure; low-pressure zone; high-pressure zone; rock pressure; solid mechanics; seismic data

INTRODUCTION

The main purpose of this work is to calculate pore pressure in a sedimentary basin, where the exploration of oil and gas takes places. Besides, we also aim at calculating low- and high-pressure zones that serve as natural pumps that accumulate fluids in the subsurface. To meet this goal, we calculate changes in stress (and, therefore, pressure types) in the subsurface based on surface seismic data.

We compare two methods of pore pressure calculation. In the first case, we calculate stress due to the vertical loading created by the gravitational geological overburden. In the second case, we propose a more complex method to calculate the stress distribution based on the mechanics of solids under gravity loading, using the concept of the first tensor invariant and the linear behavior of Hooke’s law.

We admit that the mechanics of solids is the proper method to describe the stress state in sedimentary basins and predict low- and high-pressure zones, where a reservoir is positioned in a low-pressure zone and capped by a high-pressure zone. Also, it is proper to consider a normal compaction trend (NCT) in the modeling and that the NCT curve is specific to a given lithology.

This work is more than a numerical experiment dedicated to pore pressure prediction, where the aim is to map pore pressure variations in the subsurface; this work is also to intended to show low- and high-pressure zones and correlate them to the mapped pore pressure profile. The method was applied to the Marmousi data described by Verrier and Branco (1972), Martin (2004), and Martin et al. (2006), among other authors. Their data supplies the information to the construction of stress prediction profiles using V_p and V_s velocities and density ρ distributions as the necessary data, as described by Sibiryakov and Sibiryakov (2015). Therefore, it is natural to classify the method as post-migration processing.

The method starts with the lithological identification for the construction of a porosity model, followed by least-squares regression to obtain the NCT. Then, we calculate the hydrostatic pressure, the alternative overburden weight, and pore and effective pressures.

The applied specific method is described from exploration company problems and the research literature classify pore pressure prediction as a relevant subject in oil and gas exploration. Besides, it consists of porosity and gravity loading phenomena, as described by Zhang (2013). Carcione and Helle (2002) and Carcione et al. (2003) address the subject of pore pressure based on a physical model that correlates porosity and elastic parameters that behave in the form of a depth-dependent exponential function, where the systematic use of the Hashin-Strickman model is made to describe upper and lower bounds for the rock parameters, as summarized by Mavko et al. (1999).

A method for calculating pore pressure is based on the NCT exponential fitting as described by Athy (1930) and Zhang (2011), where they call attention to the lithologies: clay, sandstone, and limestone. Several efforts are made for predicting pore pressure using porosity data, and Heppard et al. (1998) used an empirical porosity equation that is like Eaton’s method using shale porosity results to predict pore pressure. We can also add the importance of this subject to engineering aspects, as given by Holbrook et al. (2005), where they describe the use of petrophysical data and stress-strain relationships for pore pressure prediction in real-time drilling. Schneider et al. (2009) also describe porosity-stress relationships to predict over-pressure in claystone.

The prediction of low- and high-pressure zones and pore pressure variation in sedimentary basins for gas and oil exploration has many theoretical aspects related to engineering, geology, and geophysics, aiming to characterize potential reservoirs in the subsurface. Furthermore, the prediction and monitoring of abnormal pore pressure is a concern to avoid serious drilling hazards. Abnormally high pore pressures can cause blowouts, besides inducing geologic disasters, such as mud volcano eruptions.

METHODOLOGY

Solid mechanics and seismic modeling

Sibiryakov et al. (2020) describe applications of the theory of solid mechanics for pressure prediction in sedimentary basins, where the purpose is to show a mechanism for driving fluids from high- to low-pressure zones. Besides, it is necessary to analyze the participation of the so-called “effective pressure” in this fluid mechanism. As a first step in this direction, we start by defining the term “pressure”, which is a phenomenon present in fluids (oil and gas) and solid porous rock (matrix, skeleton), being a concept directly related to stress and gravity loading. Thus, pressure stands for a scalar quantity, and stress for a vectorial field.

The size of the hydrostatic pressure,P_h, is controlled by the fluid column (h = z − z₀) from the surface (z₀) to the depth of interest (z), the fluid density [ρ_f (z) ≈ ρ_f] and the gravity acceleration [g(z) ≈ g], expressed by the overload formula given by Equation 1.

(1) $$P_h=g{\int }_{z_0}^z{\rho }_f\left(z\right)dz\approx {\rho }_{f}gh.$$

The dimensions (and units) of pressure and stress are always [FL²], but the physical concepts can be different.

The overburden weight plays the role of vertical normal stress (P_z = σ_z) and it is defined by the total weight of the rock formation and fluids down to the reference point, usually taken around the reservoir. This model does not integrate the lateral variation in physical parameters, like density. Therefore, considering that the bulk density varies only with the depth, the overburden stress can be calculated by integrating the density function to the depth z, i. e. (Eq. 2):

(2) $$P_z\left(z\right)=g{\rho }_{\mathrm{W}}{h}_{\mathrm{W}}+g{\int }_{z_{\mathrm{w}}}^z{\rho }_{\mathrm{b}}\left(z\right)dz$$

where we included a liquid layer (over the rock basin formations) with water density ρ_w, height of water column h_w, seabed depth z_w (for onshore drilling z_w = 0 and h_w = 0), and depth of interest z. The gravity effect g = g(z) is considered constant (very small variation) in the interval (z_w, z).

It is relevant to note that for a horizontal plane stratified earth, the horizontal components of stress are null, and the loading action is due only to σ_zz. This means that the structure is constrained laterally and that the normal are symmetrical components. Therefore, as an approximation, only gentle dips should be admitted in the model.

The formation rock bulk density function, ρ_b (z), is admitted to vary only with depth (as a first approximation) and given by Equation 3:

(3) $${\rho }_b\left(z\right)=[1-\varphi (z\left)\right]{\rho }_m\left(z\right)+\varphi \left(z\right){\rho }_f\left(z\right).$$

Where:

ϕ = the porosity of the rock;

ρ_m = the matrix (or mineral) density;

ρ_f = the fluid (water, oil, or gas) density.

If the fluid is water, ρ_f is the density of formation water, which usually varies from 1.0 to 1.08 g/cm³ in sedimentary basins. The bulk density, ρ_b, can also be obtained from the density log, but these measures are not always available, or only for a borehole representing a limited area of the basin; therefore, we can make use of empirical prediction models.

There is confusion about the use of the word “pressure” in the geological and geophysical literature; therefore, it is convenient to add a second word to complement the concept behind the word “pressure” to make the physical description clearer. This terminology issue may derive from the fact that stress and pressure have the same physical dimension (force per unit area). The term “pressure” is clear for the hydrostatic (a scalar quantity) case. However, for rocks, we have nine elastic stresses (a second-order tensor) in the mechanics of solids.

Some authors use the term rock pressure, P_r, (instead of pressure), and we did the same in this paper, defined as the simple average of the normal stresses, and this corresponds to one-third of the first stress tensor invariant (I₁ = σ_xx + σ_yy + σ_zz) given by Equation 4:

(4) $$P_r=\frac{1}{3}\left({\sigma }_{xx}+{\sigma }_{yy}+{\sigma }_{zz}\right).$$

For the hydrostatic case, σ_xx = σ_yy = σ_zz, then, Pr has a clear physical meaning as a scalar. For the 3D geological medium, admitting the case of a flat-stratified earth model, the stress is non-hydrostatic and submitted to the symmetry condition σ_xx = σ_yy. In general, we can consider that σ_xx, σ_yy, and σ_zz can all have different values. We call the relevant vertical normal component stress σ_zz = P_z (Persen 1965).

The admitted state condition is of perfect elasticity where the Hooke’s law correlates linearly stress to strain, and isotropy for the stress tensor, σ_ij, given by Equation 5 (Fung 1965):

(5) $${\sigma }_{ij}=\lambda \theta {\delta }_{ij}+2\mu {\epsilon }_{ij}$$

for one specified material, and where μ and λ are the Lame’s elastic parameters, θ = ∇. u the dilation, δ_ij the Kronecker delta, ${\epsilon }_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$ the deformation tensor, and u = (u_i, i = 1,2,3; x,y,z) the displacement vector. The main stress components are σ_xx = σ_yy = λε_zz, σ_zz = (λ + 2μ)ε_zz, and the tangential components of the stress tensor are zero, σ_ij = 0, (i ≠ j) (Sibiryakov et al. 2004, Vieira et al. 2017).

Horizontal stress

Another problem with the concept of gravitational pressure on rock formations is that rock pressure (4), P_r, is simply taken as the loading σ_zz; that is, P_r ≡ σ_zz, without taking into account horizontal stress components σ_xx and σ_yy. Since we include this effect in our definition of rock pressure, we can show that the horizontal normal stress components (σ_xx and σ_yy) are sufficiently lower than σ_zz and given by Equation 6:

(6) $${\sigma }_{xx}={\sigma }_{yy}=P_s\left(1-2{\gamma }^2\right)=P_x$$

where we consider the lateral symmetry σ_xx = σ_yy, and the velocity ratio $\gamma =\frac{V_s}{V_p}$ (Sibiryakov et al. 2020).

Thus, using an algebra method with Equations 4 and 6 and using Hooke’s law (for the 3D case), we obtain an expression to represent the rock pressure as given by Equation 7:

(7) $$P_{rv}=P_z\left(1-\frac{4}{3}{\gamma }^2\right)$$

where the seismic information enters through the γ ratio. Equation 7 shows that P_rv is significantly reduced concerning the simple overlying load P_z = σ_zz, since γ < 1. For example, for the Poisson solid, μ = λ, γ = 0.57735; $\frac{P_{rv}}{P_r}\approx \frac{1}{2}$, i. e., P_rv corresponds approximately to half the value of P_r, which is a big difference. This means the rock pressure, P_rv, depends strongly on the velocity γ ratio (Sibiryakov E.P. et al. 2015). It should be clear that pressure on solids (rocks) does not include the effect of fluid (liquid and gas) in the different definitions.

From the rationale around Equations 4 and 7, the concept of rock pressure (or simply pressure in solids under gravity effect) does not have a clear physical meaning; nevertheless, it can be used to represent a loading force for driving fluids. A similar line of thought was honored by Cibin et al. (2011). A relevant observation is the limitations of the models since the layers are considered horizontal (or with a very smooth dip); also, fracturing and chemical transformations are not admitted in the modeling.

As a partial conclusion, loads P_z (x, z), P_r (x, z), and P_rv (x, z) should be compared as models for the driving fluid mechanism in the effective pore pressure concept, where P_z (x, z) is dependent on the density distribution, and P_rv (x, z) needs the addition of information about V_p (x, z) and V_s (x, z) seismic velocities. But, these loads should be calculated without considering the presence and effect of fluids.

Effective and pore pressures

Xu et al. (2006) and Li et al. (2009) showed in laboratory how transversal, V_s, and longitudinal, V_p, wave velocities in micro-inhomogeneous medium containing fluids depend on the difference between the rock matrix pressure [the three possibilities are P_z, P_r, or P_rv, all denominated by P_sca (and “sca” from scalar)] and the fluid pressure (or pore pressure, P_p).

The effective (or differential) pressure, P_eff, is the net pressure that acts on the solid rock matrix (skeleton) formed by the grains. According to Terzaghi (1943)’s principle, P_eff is defined as the difference between a general scalar overburden pressure, P_sca, and the pore pressure P_p, namely (Eq. 8):

(8) $$P_{\mathrm{eff}}=P_{\mathrm{sca}}-P_{p.}$$

We must assert that P_p is a natural scalar and that P_sca is also a scalar under the definition of the mean of the first stress invariant. The experimental results from Xu et al. (2006) and Li et al. (2009) show that sometimes it is necessary to introduce a dimensionless variable factor, κ, of unknown physical nature, based only on empirical evidence to reduce the participation of P_p for better data fitting, so that Equation 8 is rewritten as Equation 9:

(9) $$P_{\mathrm{eff}}=P_{\mathrm{sca}}-\kappa P_p.$$

It should also be clear that the model does not consider fluid flow.

Our proposed alternative is to change Equations 9 and 8 and use the rock pressures P_rv (and P_r) in Equations 8 and 9. We can obviously rewrite them as P_eff = P_rv − P_p and P_eff = P_rv − κP_p (we do not estimate κ in this work). The meaning of effective pressure is related to the part of the stress tensor which produces measurable effects such as compaction of sedimentary rocks, or an increase in shearing resistance. If P_eff is reduced, the compaction rate and pressure generation are also reduced; we can argue that P_eff can be considered a driving fluid mechanism.

As a basic principle, there are bounds for the variation of pore pressure in a sedimentary basin. As a reference, pore pressure can vary from the hydrostatic pressure to more than two times the hydrostatic pressure in a rock formation. When the pore pressure is equal or approximately equal to the hydrostatic pressure, the sediments are considered normally pressured, and the pore pressure is called normal pore pressure.

Following this principle, a usual qualitative description of the geological model for the sedimentation phenomenon and pressure variation with depth is as follows, and depicted in Figure 1. Then, the question is how to transform this qualitative model into a quantitative description. The initial marine sediments are formed by unconsolidated material and have high permeability and porosity. This means that the water is not confined in the pore spaces, and has a pressure communication with the seafloor surface. While compaction occurs, the sediments are buried at greater depths, and the increase of overburden load acting in the grain contacts results in lower porosity, lower permeability, and fluid ejection. If the sedimentation rate is higher than the rate at which fluids are expelled from the pores, or seals avoid dewatering of rock formations, the fluids in the pores become over pressured and bear part of the load of the overburden sediments. This compaction disequilibrium generates abnormally high pore pressure zones, with higher porosity than the normally compacted formations. Figure 1 illustrates this simple description; but, here we look at details that can contribute to more realistic pressure and stress distributions in a sedimentary basin based on seismic information (Carcione and Gangi 2000).

Figure 1.
(A) A schematic representation of a porosity profile in a sedimentary basin. The blue dotted curve represents the NCT, ϕ_n(z), in a formation. The compaction disequilibrium, caused by undercompaction (black curve), generates an abnormal pore pressure profile, P_p(z), as displayed in (B), which corresponds to overpressure. The details are redrawn from Zhang (2013).

The adopted model from solid mechanics for the description of pressure variation with depth is related to the parameter variations that directly define jumps in the $\gamma =\frac{{\overline{V}}_s}{V_p}$ ratio profile (V_s and V_p are the velocities of S and P waves, respectively). Therefore, this method will take us to a more consistent physical model to explain geological phenomena and, consequently, the engineering potential applications as described by Sibiryakov B.P. et al. (2015), Sibiryakov E.P. et al. (2015), and Vieira et al. (2017).

As further details, Sibiryakov (2014) and Sibiryakov and Sibiryakov (2018) addressed the problem of pore pressure and conditions for the disjoining pressure as the solution of a boundary value problem of the mixed type, with new computer technology for solving three-dimensional problems for the elastodynamic equation under stationary oscillations. The use of this new approach allowed the analysis of a few problems with the driving pressure in the granular medium, and they showed that the variation of properties of the static stress state under the influence of pore pressure depends essentially on the contact geometry and only a little on the type of boundary condition at the grain contact.

Estimation of porosity from density

The classical bulk density (or rock density) is defined as the total mass contained in a given unit volume of the rock. In a porous rock, this measure depends on the grain and fluid densities forming the rock, as expressed by the Equation 3. The bulk density is a relevant parameter for obtaining the overburden (or lithostatic) pressure.

In the present case, we are interested in estimating porosity from the density profile since it is a more direct measurement. Therefore, isolating for the porosity in the Equation 3, we obtain the Equation 10:

(10) $${\left.\varphi \left(z\right)\right|}_{\left(x\text{ fixed }\right)}=\frac{{\rho }_m\left(z\right)-{\rho }_b\left(z\right)}{{\rho }_m\left(z\right)-{\rho }_f\left(z\right)}.$$

The porosity estimate from the Equation 10 depends on the knowledge of the fluid density in the formations, ρ_f, and the density of the matrix (or main mineral), ρ_m, for a given rock type. Table 1 is an example of a list of matrix materials formed by common minerals. For the most solid granitic type of rocks, the matrix average density is around 2.65 g/cm³, also approximately for quartz. Clay minerals predominate in shale composition.

Table 1.
List of the average density of matrix materials formed by common minerals as presented by Zhang (2019).

Mineral	Density ρ (g/cm3)
Quartz	2.65
Calcite	2.71
Kaolinite	2.59
Dolomite	2.87
K-Feldspar	2.56
Sylvite	1.99
Halite	2.165

NCT and pore pressure from porosity and loading

To predict abnormal pore pressure generated by the compaction disequilibrium model, we need to obtain the porosity under normal compaction conditions. Athy (1930) proposed a relation for the decrease of porosity with depth under normal pressure conditions (P_z = σ_zz) with an exponential behavior given by Equation 11.

(11) $${\left.{\varphi }_n\left(z\right)\right|}_{\left(x\text{ fixed }\right)}={\varphi }_0{e}^{-cz},$$

Where:

ϕ_n = the porosity in normal compacted formation;

ϕ₀ = the porosity at the mudline (seafloor);

c = a constant that represents the compaction factor;

z = the true vertical depth below the sedimentary depositional mudline (x fixed).

The law expressed by Equation 11 is widely applied for shales because their normal porosity profiles generally show a concave downward curvature, as shown in Figure 1. Therefore, the main goal when analyzing compaction trends is to identify and evaluate the parameters that control the mathematical expression for normally compacted sediments, as in the case of shales. In the data used in this work, we have a variety of limestones, sandstones, and shales.

There are several methods for pore pressure prediction from formation porosity, and here we followed the method proposed by Zhang (2011) with some modification, which derived from the NCT of porosity proposed by Athy (1930) using Equation 11. Following this method, the pore pressure, the overburden pressure, and the porosity relate to one another by the optional expressions given by Equations 12 and 13:

(12) $${\left.P_p^{\left(z\right)}\left(z\right)\right|}_x=P_z\left(z\right)-\left[P_z\left(z\right)-P_h\left(z\right)\right]\frac{\ln {\varphi }_0-\ln \varphi \left(z\right)}{cz},$$

(13) $${\left.P_p^{\left(rv\right)}\left(z\right)\right|}_x=P_{rv}\left(z\right)-\left[P_{rv}\left(z\right)-P_{\mathrm{h}}\left(z\right)\right]\frac{\ln {\varphi }_0-\ln \varphi \left(z\right)}{cz},$$

Where:

P_p(^rv) and P_p^(z) = the pore pressure models;

P_rv = our present proposed overburden pressure model;

P_h = the hydrostatic pressure (or the normal pore pressure);

ϕ₀ = the initial porosity (at seabed level);

ϕ = the porosity at the depth of interest z and c is compaction factor.

The main point in the interpretation for the Equations 12 and 13 is that the calculated pore pressure P_p(z) has porosity as depth-dependent. If the porosity ϕ(z) (at a target depth z) is greater than the normal porosity ϕ_n for the same depth, this zone corresponds to an over-pressure; otherwise, to an under-pressure.

From the concept of effective pressure, in this work, we compare the Equations 12 and 13 calculated with P_sca = P_rv and with P_sca = P_z and present as results. The vast literature uses only P_z for the overburden loading in the geological description for the mechanics of sedimentation. However, from the concepts of mechanics of solids, this description is incomplete due to the horizontally developed stresses.

Flowchart for pressure distribution calculation

The Marmousi2 data (an update of Marmousi data) is a collection of seismic distribution of the P and S wave velocities [V_p(x, z) and V_s(x, z), respectively], and the distribution of density [ρ(x, z)] built from a region of the Cuanza basin (Angola), as described in the sequel. The basin is characterized by intense tectonics, a clear sequence of lithological units with well-defined layer geometries. In such a way, with the distribution of physical parameters and geometry, forward modeling allows different studies in a complex and important geological environment.

The following flowchart and block diagram in Figure 2 describe the main steps for the modeling of pore and effective pressures applied to the Marmousi2 seismic data, and it represents a summary of the methodology. We use a 3D grid (Δx, Δy, Δz) for all calculations, but the transversal Y direction is formed by only one cell (Δy), while along with the (x, z) coordinates we have (M x N) grid cells (brick wall type of profile).

Figure 2.
Block diagram showing the main steps for pore pressure prediction and the relationships between the input data, calculated parameters, and pressure sections.

Flowchart

• Field and seismic parameters of the migrated seismic section;

• The input of V_p(x, z), V_s(x, z), ρ(x, z) distributions in an equal mathematical matrix format;

• Separation of the main lithologies in the seismic section;

• Calculation of the porosity distribution;

• Least-square fitting of the NCT based on the chosen geological formation: definition of the parameters ϕ₀ (porosity at the mudline level) and c (compaction factor for the geological formation);

• Calculation of different pressure models: hydrostatic P_h(x, z), overburden P_z(x, z), and rock P_rv(x, z) as a function of the velocity γ(x, z) ratio;

• Calculation in loop of the pore pressure P_p^(z)(x, z) for the profiles using Equation 12;

• Calculation in loop of the pore pressure using rock pressure P_p^(rv)(x, z) for the profiles using Equation 13;

• Calculation of effective pressure P_eff^(z)(x, z) using Terzaghi’s Equation 8;

• Calculation of effective pressure P_eff^(rv)(x, z) using Terzaghi’s Equation 8;

• A Comparison of the different results aiming at the target zones.

RESULTS

We present our findings in a sequence of color figures that also show the modeling steps taken in mapping the low- and high-pressure zones. In general, the target is to point to new potential oil and gas accumulations or extend production reservoirs.

Geological and seismic information

The Marmousi model was created at the Institut Français du Pétrole (IFP) in 1988, based on a profile through the North Quenguela Trough in the Cuanza basin (Angola), as described in the work by Verrier and Branco (1972). The Marmousi model was developed to meet several research demands formulated at the time, to which both the model and the data would need to adhere. First, the model should be geologically plausible, preferably based on a real situation. Second, the model should be complex, containing several reflectors, steep dips, and strong velocity gradients in both vertical and lateral directions. Thus, the generated data should be so complex that it would not satisfy the assumptions of conventional processing, requiring depth migration methods to get good images.

The model used in this work is an upgrade of the original Marmousi model, now called Marmousi2, as described in detail by Martin (2004) and Martin et al. (2006). This update consisted of the reconstruction and spatial expansion of the original Marmousi model, preserving its structural complexity. Besides, they added new channels and hydrocarbon traps and kept the original model reservoirs.

The lithological composition of the Marmousi model is predominantly of shale units, with occasional sand layers. The central area has large faults, with an anticline composed of marl (carbonate-rich shale). In addition, this model has an unconformity overlapped by large salt layers separating the marls from the deeper anticlinal units, which are also mostly shales with some interweaving of sand layers (Fig. 3).

Figure 3.
Structural elements, formations, and lithologies of Marmousi2 according to Martin et al. (2006). In general, the vertical exaggeration is of the order of 4:1.

The Marmousi model has a series of hydrocarbon reservoirs contained in its structural model, displayed as gas (red) and oil (green) in Figure 3. In this model, all layers containing hydrocarbons are sand. The reservoirs are distributed within a complex fault zone at varying depths, also in the layers with simple structural settings and more structurally complex locations, such as folds.

All figures relate to the complex Marmousi2 used to apply the developed methodology for pore pressure prediction based on porosity. The general vertical exaggeration of the following figures is of the order 4:1, and the interpretation principle is that all computed sections should resemble the original data in Figure 3.

The available data has seismic velocities [V_p(x, z) and V_s(x, z)] and density [ρ(x, z)] distributions, as shown in Figures 4, 5 and 6, respectively, which are the necessary starting information for calculating stress distribution. Also, another relevant information is the γ(x, z) ratio distribution shown in Figure 7, which participates in the modeling and represents the contribution of the seismic information for the prediction of low- and high-pressure zones and pore pressure.

Figure 4.
P-wave velocity distribution, V_p(x, z).

Figure 5.
S-wave velocity distribution, V_s(x, z).

Figure 6.
Density distribution, ρ(x, z).

Figure 7.
Velocity ratio distribution, γ(x, z).

In the absence of borehole data, we adopted a strategy to create a model with the lithological distinction between the main lithologies as described by Martin et al. (2006) and depicted in Figure 3. Using the information for the velocities V_p, V_s, and density ρ, we constructed a model giving a constant value for each type of lithology as described in the caption of Figure 8.

Figure 8.
A lithological separation where each integer in the caption corresponds to a rock type in the following order: 1, water (dark blue); 2, shale (blue); 3, sandstone (light blue); 4, sandstone with hydrocarbons (yellow); 5, limestone (orange); 6, salt (red wine color).

The porosity model used Equation 10 and the matrix of the model with the lithological separation as shown in Figure 8. Using the density model as an input, as displayed in Figure 6, the porosity estimate was applied at all points of density distribution, applying the matrix density values given in Table 1 for each type of lithology. Figure 9 shows a wide range of porosity values, which highlights the large contrast between the sandstone reservoirs and the shales in the region between 500 and 2,000 meters of depth.

Figure 9.
Porosity section, ϕ(x, z), calculated using Equation 10.

Hydrostatic and overburden pressure sections

Figure 10 shows the hydrostatic pressure P_h(x, z) corresponding to a linear trend according to Equation 1 where the constant value for water density used was ρ_f = 1.01 g/cm³. Figure 11 shows the variation of pressure P_z according to Equation 2. The outline of map 11 is similar to 10, but it shows some details, as can be seen in the central region, where the complex fault zone is located. We should keep in mind the limitations due to the color scale of the figures.

Figure 10.
Hydrostatic pressure section, P_h(x, z), according to Equation 1. The map is characterized by its smoothness.

Figure 11.
Overburden pressure section, P_z(x, z), according to Equation 2. We can see some very light variations with respect to Figure 10.

In the other hand, Figure 12 represents a partial conclusion of this work, where shows the rock pressure distribution P_rv(x, z) according to our Equation 7. We notice in this figure a very different outline with respect to the previous described pressure sections 10 and 11.

Figure 12.
Rock pressure section, P_rv(x, z), according to Equation 7, where we can trace the characteristics of the original geological section, and some strong variations with respect to Figures 10 and 11. The rectangle around 3,000 meters indicates the main reservoir of the basin, while the arrows indicate other smaller reservoirs, and they correspond to low-pressure zones.

The outline shown in Figure 12 is due to the contribution of lateral variations of horizontal stress, σ_xx = σ_yy, which depend on $\gamma =\frac{V_s}{V_p}$ ratio for calculating P_rv. A more complete modeling would consider the condition of non-symmetry in σ_xx, σ_yy, and σ_zz stresses.

Porosity profiles and NCT estimation

The next step was to select some ordered porosity profiles for specific analysis, taken from the calculated model displayed in Figure 9, and shown in Figure 13. These profiles are positioned starting at 10,500 m from the left border of the model, and they cross two reservoirs in different depths. The first reservoir is located in the complex central fault region, around 1,000 m depth, and the second is located at the top of the lower anticlinal and around 3,000 m depth. We selected for exemplification only the region with the presence of shale to find the best fitting exponential curve to the porosity values. That is, the purpose of this step is to find the NCT involving the shale formations.

Figure 13.
The selected profiles (vertical lines) over the porosity distribution as displayed in Figure 9. The profiles serve for further analysis of the anticline region that includes the main reservoir (see Fig. 14). The orange lines are placed over the reservoir and the white lines outside and near the reservoir.

Figure 14 shows details of the selected profiles across the main reservoir, where the porosity was calculated using Equation 10. The green line represents the fit in the least-squares sense between Athy (1930)’s exponential model and the red dots (shale formations). The obtained values for the model were ϕ₀ = 0.43367 (fixed porosity for the mudline), and c (compaction factor) is variable (example, c = 0.0006773 m^-1 for one line). Furthermore, we can see that Athy’s exponential model (NCT) is suitable for a general representation of the variation of porosity with depth. Besides, based on this model, the line segments positioned to the right side of the green reference line indicate the condition of higher pore pressure and, to the left, of lower pore pressure. Also, the reservoir is located at a depth of around 3,000 meters, and we can see that the condition is systematic of higher pressure, but it does not stand out with high values as in other places of the section.

Figure 14.
The porosity profiles across the main reservoir as shown in Figure 13. Blue line: complete profiles of porosity based on the lithological description and Athy’s law. Green line: the NCT for the porosity profile based on the shale formations distributed along with the profile. Red dots: highlights the presence of shale formations over the blue lines. Dashed lines: delimits in profiles P2 to P6 the reservoir zone in the anticline structure.

Pore pressure prediction sections and profiles

In the following figures, we are interested in looking for patterns, combinations, correlations, coherences, values, and caption scales, to compare, understand, and interpret the obtained sections and profiles. The combination of all information is the method for interpreting and characterizing the target zones (reservoirs, collectors).

Figures 15 and 16 represent the sections of pore pressure calculated, respectively, by Equations 12 and 13 where we can see the good coherence between them. Also, that Figure 16 more clearly displays the details of the target reservoirs, which we consider an improvement and a partial conclusion of this work. Observe that the caption of Figure 15 indicates the presence of low pore pressure and even a strange negative value (see Equation 12) as a result of high contrast for the physical parameters of the salt formation. However, this effect does not appear for the model of the Equation 13.

Figure 15.
Pore pressure section, P_p^(z)(x, z), based on Equation 12.

Figure 16.
Pore pressure section, P_p^(rv)(x, z), based on Equation 13.

The effective pressure is considered as a potential driving mechanism for fluid percolation, and a lower pressure zone, capped by a high-pressure zone, is the proper zone for fluid accumulation (Sibiryakov et al. 2020). Figures 17 and 18 represent sections of effective pressure calculated by Equation 8 with its variants 12 and 13, respectively. As in the previous case of Figures 15 and 16, Figure 18 more effectively reveals details with respect to the target reservoirs, and a better color gradient for the figure; besides, Figure 18 represents a partial conclusion in this work by more effectively highlighting the reservoirs.

Figure 17.
Effective pressure section, P_eff^(z)(x, z), based on Equation 8.

Figure 18.
Effective pressure section, P_eff^(rv)(x, z), based on Equation 8.

Figures 19 and 20 represent the seven profiles selected for detail analysis, as shown in Figure 13 across a target zone. The curves shown are calculated with their respective models where the hydrostatic, rock, pore, and effective pressures are included.

Figure 19.
Pressure profiles for Figure 13. Blue: hydrostatic linear P_h(z). Red: overburden linear P_z(z). Black: pore variable P_p^(z)(z). Green: effective variable P_eff^(z)(z).

Figure 20.
Pressure profiles for Figure 13. Blue: hydrostatic linear P_h(z). Red: overburden non-linear P_rv(z). Black: pore variable P_p^(rv)(z). Green: effective variable P_eff^(rv)(z) . The dashed lines in profiles P2 to P6 indicates the reservoir zone contained in the anticline structure.

For the set of profiles in Figure 19, the reference bounds are linear and shown in blue [hydrostatic P_h(z)] and red [overburden P_z(z)] colors. Furthermore, the pore pressure P_p^(z)(z) and the effective pressure P_eff^(z)(z) vary around these reference bounds, but mainly around the hydrostatic trend and in the opposite sense (negative correlation: P_eff = P_sca −P_p). Also, for the target reservoir located at 3,000 meters depth, the pore pressure is more indicative.

For the set profiles in Figure 20, the linear reference bound is shown in blue [hydrostatic P_h(z)]. On the other hand, the red reference bound [overburden P_rv(z)] is rather variable, but with different mathematical meaning concerning its similar in Figure 19. Furthermore, the pore pressure P_p^(rv)(z) follows the tendencies of the rock pressure P_rv(z) with a good positive correlation. Also, the effective pressure shown in green [P_eff^(rv)(z)] is systematically below the reference hydrostatic bound shown in blue, and with a negative correlation concerning the pore pressure curve shown in black [P_p^(rv)(z)]. Also, for the target reservoir located at 3,000 meters depth, the combination of the rock pressure P_rv(z) with the pore pressure P_p^(rv)(z) is more indicative for the target.

The complementary Figures 21 and 22 depict the deviation between the pore pressure models, (12) and (13), and the hydrostatic model, (1). The first observation is that these findings are very coherent, as can be seen by the color distribution. The second is that the scales of the captions show the double magnitude interval, which is explained by the important presence of the γ ratio that represents the seismic data. Therefore, we can wind up that the modeling is consistent, with these figures showing only a relative value since the difference could be the opposite.

Disclaimer: This is an incomplete paper and, therefore, represents research in progress. It is not meant to represent professional results, but rather an academic opinion of the authors. The authors are not responsible for any practical business applications.

Figure 21.
Pore pressure deviation section calculated by the difference between the pore pressure from overburden pressure and the hydrostatic model: Δ_p^(z)(x, z) = P_p^(z) (x, z) − P_h (x, z).

Figure 22.
Pore pressure deviation section calculated by the difference between the pore pressure from rock pressure and the hydrostatic model: ΔP_p^(rv)(x, z) = P_p^(rv) (x, z) − P_h (x, z).

GENERAL CONCLUSIONS

We outlined throughout the text some partial conclusions for the different pressure sections and profiles. In this section, we provide general comments only.

We demonstrated in this work the possibility of using the proposed method that uses seismic information, V_p(x, z) and V_s(x, z) velocities and density ρ(x, z), to calculate rock, pore, and effective pressure distribution. These are useful for characterizing potential zones for oil and gas accumulation. However, it is necessary to have a lithological description of the basin, and possibly use empirical and petrophysical models. The proposed macro model establishes the use of the rock pressure P_rv(x, z) in the models for different pressure calculations, instead of the simple overburden pressure P_z(x, z). We can state that the present method has adaptative characteristics in the sense that it can be reviewed as new information is incorporated in the processing.

It is relevant to consider a combination of different modeling for the analysis of the target zones. Therefore, as a macro model, we should consider driving rock pressure P_rv(x, z), pore pressure P_p^(rv)(x, z), and effective pressure P_eff^(rv)(x, z) to understand the fluid percolation mechanism between high- and low-pressure zones in the basin, together with the fundamentals of petroleum geology.

For the application of the method, it is necessary to know the porosity distribution values. Therefore, in the absence of borehole logs, an estimation of this property is made using empirical relationships to produce an optimum porosity model for the basin.

The calculation of the NCT is used for the prediction of pore pressure. Therefore, lithological separation is necessary to identify the geological formations and adjust the NCT for the same type of rock. Namely, an ill-fit curve can create a false over-pressure or under-pressure zone.

A usual geological descriptive mechanism for pore pressure generation is the compaction unbalance, as shown in Figure 1, but this description is physically limited. Furthermore, our physical modeling for the present state can be guided by the principles of petroleum geology concerning the generation and accumulation of oil and gas. In addition, the pore pressure-porosity model employed in this work is closely related to detecting changes in compaction rates (ϕ₀ and c parameters), which are reflected in porosity values.

Our main conclusions in this work are directly related to Figures 12, 16, 18, and the companion profile figures, that are based on the P_rv(x, z) pressure model. We can summarize our conclusions based on the calculated images by laying out the following sections:

• The geological information highlighting the potential reservoirs and specifying the targets;

• The rock pressure to show the low- and high-pressure zones;

• The pore pressure P_p^(rv)(x, z);

• The effective pressure P_eff^(rv)(x, z);

• Then, display the profiles over the target zones, and focus on the low-pressure zones with relative high pore pressure [P_p^(rv)(x, z)], when these conditions correspond to low effective pressure [P_eff^(rv)(x, z)]; We analyze these trends taking the hydrostatic linear trend [P_h(x, z)] as reference;

• To extend the reservoir target zone, it can be done laterally and vertically.

The combined analysis of these three figures will identify patterns and correlations, outline and characterize the target zones, and make conclusions.

The modeling of 3D data requires more computation power. However, a complete application of this methodology needs 3D data where we can calculate the stress distributions in full, and no symmetry is adopted.

APPENDIX

We used least squares regression to calculate the NCT, shown in Figure 14, to fit the porosity values corresponding to the shale formations. The observed data consist of N pairs of points (ϕ_i^(sh), z_i^(sh)), where z_i^(sh) is the depth of the shale formation, and ϕ_i^(sh) is the corresponding porosity value. The chosen regression model is the classical one due to Athy (1930)’s compaction that depends only on the vertical loading and is defined by Equation 14.

(14) $${\varphi }_n\left(z\right)={\varphi }_0{e}^{-cz},$$

Where:

ϕ_n^(z) = the NCT curve;

z = the depth below the mudline level;

ϕ₀ = the porosity at the mudline;

c = the compaction factor, both estimated from the curve fitting.

Then, applying the natural logarithm operator to the observed and model (Equation 14) data, we transform the expressions into a linear relationship, namely Equation 15 and Equation 16:

(15) $$y^{\left(obs\right)}=\ln {\varphi }^{\left(sh\right)},$$

(16) $$y^{\left(m\right)}=a+bz,$$

where y^(m) = lnϕ_n, a = lnϕ₀, b = −c. Figure 14 shows y^(obs) by red dots and y^(m) by a green line.

The misfit between the observed and calculated data is defined by Equation 17.

(17) $$r_i=y_i^{\left(obs\right)}-y_i^{\left(m\right)},$$

and the optimizing functional by Equation 18.

(18) $$S(a,b)=\sum _{i=1}^N{\left[y_i^{\left(obs\right)}-a-b{z}_i\right]}^2,$$

which is a concave positive and from where the parameters a and b are estimated. The minimum of Equation 18 is obtained by deriving S(a, b) with respect to a and b and equalizing to zero, i. e. (Eq. 19).

(19) $$\begin{gathered} \frac{\mathrm{\partial }S(a,b)}{\mathrm{\partial }a}=-2\sum _{i=1}^N \left[y_i^{\left(abs\right)}-a-b{z}_i\right]=0 \\ \frac{{\displaystyle \partial S(a,b)}}{{\displaystyle \partial b}}=-2\sum _{i=1}^N{z}_i\left[y_i^{\left(abs\right)}-a-b{z}_i\right]=0 \end{gathered}$$

which is a simultaneous system of Equation 19 with the solution given by Equation 20.

(20) $$\begin{gathered} a={\overline{y}}^{\left(obs\right)}-b\overline{z}, b=\frac{{\sum }_{i=1}^N{z}_i\left[y_i^{\left(obs\right)}-{\overline{y}}^{\left(obs\right)}\right]}{{\sum }_{i=1}^N{z}_i\left(z_i-{\hat{z}}^{\text{'}}\right)} \\ \overline{z}=\frac{1}{N}\sum _{i=1}^N{z}_i, {\overline{y}}^{\left(obs\right)}=\frac{1}{N}\sum _{i=1}^N{y}_l^{\left(obs\right)}, \end{gathered}$$

from where by inversion, ϕ₀ = e^a, and c = −b.

Depending on the description and classification of the given data, the optimizing functional can have a different expression. In this case, the estimated parameters can be different, but not that different from the least-square fit.

ACKNOWLEDGEMENTS

The authors would like to thank Universidade Federal do Pará (UFPA), the National Institute of Science and Technology (INCT-GP, MCT/CNPq/FINEP), and PETROBRAS/CENPES for the research support (GEOMEC project) aiming at oil and gas exploration. The authors would like thank the reviewers for their patience and positive contributions.

Article Information

• Manuscript ID: 20200131.
  How to cite: Leite L.W.B., Andrade F.T.B. Pore pressure prediction based on rock stress applied to the Marmousi seismic data. Brazilian Journal of Geology, 52(3):e20200131, 2022. https://doi.org/10.1590/2317-4889202220200131

REFERENCES

• Athy L. 1930. Density, porosity and compaction of sedimentary rocks. Bulletin of the American Association of Petroleum Geologists, 14(1):1-24. https://doi.org/10.1306/3D93289E-16B1-11D7-8645000102C1865D
  » https://doi.org/10.1306/3D93289E-16B1-11D7-8645000102C1865D
• Carcione J.M., Gangi A.F. 2000. Non-equilibrium compaction and abnormal pore-fluid pressures: effects on rock properties. Geophysical Prospecting, 48(3):521-537. https://doi.org/10.1046/j.1365-2478.2000.00197.x
  » https://doi.org/10.1046/j.1365-2478.2000.00197.x
• Carcione J.M., Helle H.B. 2002. Rock physics of geopressure and prediction of abnormal pore fluid pressures using seismic data. CSEG Recorder, 27(7):8-32.
• Carcione J.M., Helle H.B., Pham N.H., Toverud T. 2003. Pore pressure estimation in reservoir rocks from seismic reflection data. Geophysics, 68(5):1569-1579. https://doi.org/10.1190/1.1620631
  » https://doi.org/10.1190/1.1620631
• Cibin P.P., Mantovani M., Moglia A. 2011. Pore pressure prediction from vp/vs ratio volumes. In: EAGE Conference and Exhibition Incorporating SPE EUROPEC, 73., 2011. Annals … Vienna: EAGE.
• Fung Y.C. 1965. Foundations of solid mechanics London: Prentice-Hall.
• Heppard P., Cander H., Eggertson E. 1998. Abnormal pressure and the occurrence of hydrocarbons in offshore eastern trinidad, west indies. Abnormal Pressures in Hydrocarbon Environments: AAPG Memoir, 45:215-246. https://doi.org/10.1306/M70615C13
  » https://doi.org/10.1306/M70615C13
• Holbrook P., Maggiori D., Hensley R. 2005. Real-time pore pressure and fracture gradient evaluation in all sedimentary lithologies. Society of Petroleum Engineers, 10(4):215-222. https://doi.org/10.2118/26791-PA
  » https://doi.org/10.2118/26791-PA
• Li M., Xiao W-L., Guo X., Zhang L., Zheng L-L. 2009. Laboratory study of the effective pressure law for permeability of the low-permeability sandstones from the Tabamiao area, inner Mongolia. Chinese Journal of Geophysics, 52(6):1402-1413. https://doi.org/10.1002/cjg2.1465
  » https://doi.org/10.1002/cjg2.1465
• Martin G.S. 2004. The Marmousi2 model, elastic synthetic data, and an analysis of imaging and AVO in a structurally complex environment MS Dissertation, University of Houston, Houston, 213 p.
• Martin G.S., Wiley R., Marfurt K.J. 2006. Marmousi2: An elastic upgrade for marmousi. The Leading Edge, 25(2):113-224. https://doi.org/10.1190/1.2172306
  » https://doi.org/10.1190/1.2172306
• Mavko G., Mukerji T., Dvorkin J. 1999. The rock physics handbook London: Cambridge University Press.
• Persen L.N. 1965. Rock dynamics and geophysical exploration Amsterdam: Elsevier.
• Schneider J., Flemings P.B., Dugan B., Long H., Germaine J.T. 2009. Overpressure and consolidation near the seafloor of brazos-trinity basin iv, northwest deepwater gulf of Mexico. Journal of Geophysical Research, 114(B5):1-13. https://doi.org/10.1029/2008JB005922
  » https://doi.org/10.1029/2008JB005922
• Sibiryakov B.P., Leite L.W.B., Sibiryakov E.P., Vieira W.W.S. 2015. Local low pressure areas in anticline structures. Brazilian Journal of Geophysics, 33(2):225-236. https://doi.org/10.22564/rbgf.v33i2.716
  » https://doi.org/10.22564/rbgf.v33i2.716
• Sibiryakov B.P., Sibiryakov E.B., Glebov A.F., Nesterov V.N., Sokolov E.P. 2004. Prediction of stress and hydrodynamic behavior of fluids from multicomponent seismic data. Russian Geology and Geophysics, 6:709-715.
• Sibiryakov B.P., Sibiryakov E.B., Leite L.W.B. 2020. Dynamics of underground rocks containing fluids: application to exploration geophysics with emphasis on oil and gas. Belém: UFPA, IG, Geosciences Library, 360 p.
• Sibiryakov E.B. 2014. Modification of the method of boundary integral equations for the solution of mixed problems of elasticity theory. (Manuscript), p. 1-20.
• Sibiryakov E.B., Sibiryakov B.P. 2015. Local pressure lows as possible sinks of fluids in geologic structures. Russian Geology and Geophysics, 56(7):1091-1095. https://doi.org/10.1016/j.rgg.2015.06.010
  » https://doi.org/10.1016/j.rgg.2015.06.010
• Sibiryakov E.B., Sibiryakov B.P. 2018. Destruction of granular media under pore pressure. IOP Conference Series: Earth and Environmental Science, 134:012061. https://doi.org/10.1088/1755-1315/134/1/012061
  » https://doi.org/10.1088/1755-1315/134/1/012061
• Sibiryakov E.P., Leite L.W.B., Vieira W.W.S. 2015. Behavior of stresses in structures and the effect on hydrodynamics analyzed from multicomponent seismic data. Brazilian Journal of Geophysics, 33(1):57-70. https://doi.org/10.22564/rbgf.v33i1.601
  » https://doi.org/10.22564/rbgf.v33i1.601
• Terzaghi K. 1943. Theoretical soil mechanics New York: John Wiley and Sons, 510 p.
• Verrier G., Branco F.C. 1972. La fosse tertiaire et le gisement de quenguela-nord. Revue de l’Institut Francais du Petrole, 27:1-72.
• Vieira W.W.S., Leite L.W.B., Sibiryakov B.P. 2017. Modeling and pressure prediction of a block of the Jequitinhonha basin. Revista Brasileira de Geofísica, 35(3):173-186.
• Xu X., Hofmann R., Batzle M., Tshering T. 2006. Influence of pore pressure on velocity in low-porosity sandstone: Implications for time-lapse feasibility and pore-pressure study. Geophysical Prospecting, 54(5):565-573. https://doi.org/10.1111/j.1365-2478.2006.00569.x
  » https://doi.org/10.1111/j.1365-2478.2006.00569.x
• Zhang J. 2011. Pore pressure prediction from well logs: methods, modifications, and new aproaches. Earth Science Reviews, 108(1-2):50-63. https://doi.org/10.1016/j.earscirev.2011.06.001
  » https://doi.org/10.1016/j.earscirev.2011.06.001
• Zhang J. 2013. Effective stress, porosity, velocity and abnormal pore pressure prediction accounting for compaction disequilibrium and unloading. Marine and Petroleum Geology, 45:2-11. https://doi.org/10.1016/j.marpetgeo.2013.04.007
  » https://doi.org/10.1016/j.marpetgeo.2013.04.007
• Zhang J. 2019. Applied petroleum geomechanics Gulf Professional Publishing, 1, 534 p.

Publication Dates

• Publication in this collection
  24 June 2022
• Date of issue
  2022

History

• Received
  18 Dec 2020
• Accepted
  21 Mar 2022