Pore pressure prediction in orthotropic medium based on rock physics modeling of shale gas

Journal Pre-proof Pore pressure prediction in orthotropic medium based on rock physics modeling of shale gas Ting Lei, Xingyao Yin, Zhaoyun Zong PII:

S1875-5100(19)30343-9

DOI:

https://doi.org/10.1016/j.jngse.2019.103091

Reference:

JNGSE 103091

To appear in:

Journal of Natural Gas Science and Engineering

Received Date: 28 August 2019 Revised Date:

21 November 2019

Accepted Date: 22 November 2019

Please cite this article as: Lei, T., Yin, X., Zong, Z., Pore pressure prediction in orthotropic medium based on rock physics modeling of shale gas, Journal of Natural Gas Science & Engineering, https:// doi.org/10.1016/j.jngse.2019.103091. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier B.V.

Ting Lei: Conceptualization, Methodology, Software, Writing - Original Draft, Visualization. Xingyao Yin: Validation, Investigation, Resources, Supervision. Zhaoyun Zong: Validation, Data Curation, Writing - Review & Editing.

Pore pressure prediction in orthotropic medium based on rock physics modeling of shale gas

Ting Lei1,2, Xingyao Yin1,2, Zhaoyun Zong1,2 China University of Petroleum, Qingdao, Shandong, China. Email: [email protected] 2 Laboratory for Marine Mineral Resources, Qingdao National Laboratory for Marine Science and Technology, Qingdao, 266071, China 1

Abstract: Pore pressure prediction is critical to reducing drilling costs and exploring hydrocarbons. Current seismic approaches for pore pressure prediction mainly initially establish the relation between vertical effective stress and physical parameters (velocity, resistivity or porosity), and then predict pore pressures with the known effective stress theory. The commonly used Eaton’s equation requires the estimation of normal compaction trend (NCT), which cannot be measured but only be inferred. However, the conventional NCT is a curve which is able to describe the normal trend but is relatively rough, limiting the accuracy of pore pressure prediction. In this study, the shale rock physics modeling in orthotropic (OA) medium is designed to equate the NCT. The modeled and the measured P-wave velocity match well in normal pore pressure segments, but they differ in abnormal pore pressure segments. The modeled NCT is utilized to predict pore pressures, and the comparison with the measured pore pressures proves that the accuracy can be improved. Then P-wave velocity is inverted from azimuthal seismic data based on the reflectivity approximation of OA medium, and the pore pressure of the research area is predicted using the modeled P-wave velocity (the target play is a shale gas field in Sichuan basin of China). The prediction result is consistent with the characteristics of overpressure mechanism of hydrocarbon generation, indicating the effectiveness of rock physics modeling in pore pressure prediction. Keywords: Pore pressure, rock physics modeling, orthotropic medium, seismic inversion. 1 Introduction The development of shale gas in North America led to a global upsurge in literatures for shale gas (Zou et al., 2013). Overpressures are common in sedimentary basins, and disequilibrium compaction overpressures are often identified by abnormally high porosity (Serebryakov et al., 2002, Swarbrick et al. 2002). Whereas, overpressures due to fluid expansion are distinguished by unloading curve (Bowers,

1995; Tingay et al., 2009; Zhang, 2011). Accurate pore pressure prediction helps to improve drilling safety and reduce drilling costs including prevention of lost circulation, borehole instability, stuck pipes, well kicks, and well blowouts et al. (Gutierrez et al., 2006; Cao et al., 2017). These accidents are usually caused by improper mud weight allocation due to the erroneous pore pressure predictions. According to the statistics of wells in the Gulf of Mexico, about 30% of deepwater drilling costs and 24% of non-productive time are associated with geopressure (Dutta, 2002; Zhang, 2011). Besides, accurate prediction of geopressure contributes to exploration and trap analysis, and helps to control the hydrofracture propagation during hydraulic fracturing (Song et al., 2016), as well as provides calibration for basin modeling (Gutierrez et al., 2006). In particular, areas with large identifiable discontinuities that are aligned for shear failure, have higher in-situ differential stress slope, and in-situ discontinuities and rock brittleness have significant impact on microseismicity distributions induced by hydraulic fracturing (Maity and Ciezobka, 2019). Therefore, accurate prediction of pore pressure is important to reduce drilling costs,

explore

hydrocarbon,

hydraulic

fracturing,

and

improve

reservoir

characterization and development. A number of pore pressure prediction methods have been proposed in literatures, which initially establish the relation between vertical effective stress and velocity (or resistivity, porosity) (Sayers et al. 2002; Swarbrick 2012; Wang and Wang, 2015; Das and Chatterjee, 2018), and then predict pore pressures by effective stress theory (Terzaghi et al., 1996). The effective stress theory holds that the overburden stress is shared by rock skeleton and pore fluid of the underlying strata,

σv = σe + p

(1)

where, p is the pore pressure, MPa . σ v is the overburden stress, MPa . σ e is the vertical effective stress, MPa . According to the relation between vertical effective stress and physical parameters, pore pressure prediction approaches can be divided into normal compaction trend (NCT)-based and empirical approaches (Lei et al., 2018). The methods based on the NCT calculate the vertical effective stress by the ratio of the infered normal-compacted formation properties to the observed properties including seismic interval velocity, resistivity or porosity, for instance, the commonly used Eaton’s equation (Hottmann and Johnson, 1965; Eaton, 1972; Eaton, 1975). While empirical methods directly establish the empirical relations to quantify the

vertical effective stress (Bowers, 1995; Lopez ea al., 2004). Unfortunately, the NCT is hard to infer, and the existing way is to fit the NCT with acoustic transit time in adjacent sandstone wells (van Ruth et al., 2004; Tingay et al., 2009), or using the normal compaction equation derived from the combination of Wyllie equation and Athy’s normal compaction porosity (Zhang, 2011). The resulting NCT is a curve which is able to describe the normal trend but is relatively rough. As for empirical methods, there are inherent errors in the empirical relations. Han et al., (2017) utilized the clay-plus-silt (CPS) (Pervukhina et al., 2015) model to an overpressure well, and found that the modeled P- and S-wave velocities were systematically higher than those from measured ones. The CPS model calculates the velocity through the total porosity of shale and it does not take into account the effects of pore pressures on velocity. However, the pore pressure affects velocity by enhancing compliant low-aspect-ratio pores, rendering the lower elastic moduli of the rock and the enhanced or unchanged total porosity, depending on the overpressure mechanism. Therefore, the simulation results of CPS model can be regarded as the NCT for pore pressure prediction. Rock physics theory helps to link the reservoir and elastic parameters. And rock physics modeling is able to simulate the elastic properties of subsurface rocks by combing mineral composition, pores and fluids (Mavko et al., 1998; Yin and Liu, 2016). Most of current shale equivalent models involve the self-consistent approximation (SCA) model (Budiansky, 1965; Berryman, 1995) or differential effective medium (DEM) model (Norris, 1985; Zimmerman, 1990). For SCA model, it places each mineral composition into an infinitely unknown background medium and adjusts the equivalent medium to approximate elastic interactions between inclusions. Whereas, DEM model needs to set the host mineral and approximates the equivalent medium by gradually adding inclusions into the host mineral. The hypothesis of idealized ellipsoid inclusions, and the approximation of aspect ratio and elastic interaction of inclusions make them suitable for simulating the coupled clay and kerogen (Wu et al., 2012; Jiang and Spikes, 2013). However, they are both high frequency models and cannot be used for low frequency seismic conditions. The layered shale with a suite of aligned or nearly aligned fractures is usually assumed as OA medium (Guo et al., 2012). Therefore, anisotropic theory is required to implement shale modeling. Hornby (1994) gave the anisotropic forms of SCA (ani-SCA) and DEM (ani-DEM), and simulated the shale thorough the combination of ani-SCA and ani-DEM.

The shale reservoir is usually characterized by diverse minerals, complex pore types including nano-scale pores, occurrence of organic matter, and plenty of adsorbed gas (Zou et al., 2013), making the rock physics modeling tough. A representative organic-rich shale model (Zhu et al., 2012) initially calculated the moduli of isotropic mixed mineral matrix by Voigt-Reuss-Hill average, then added pores to the matrix with ani-DEM, and then added fluids by Gassmann equation. Finally, kerogen was added to the pores by the anisotropic solid replacement equation (Ciz and Shapiro, 2007). In general, how to add kerogen and multi-type pores is a challenge for the rock physics modeling of organic-rich shale. According to maturity, Zhao et al. (2016) classified kerogen into three categories: immature, mature and over-mature, and added kerogen with different models to study its effect on shale. As nanopores are developed in shale, researches on nanopores are increasing (Miller and Shenoy, 2000; Gor and Gurevich, 2018). However, shale rock physics models have not been widely applied to pore pressure prediction. Because the commonly used NCT calculation method is simple and efficient, which is derived from the combination of Wyllie equation and Athy’s normal compaction porosity (Zhang, 2011). In this study, a pragmatic shale rock physics modeling approach for OA medium is introduced. Firstly, Clay and kerogen are mixed up by ani-SCA and ani-DEM to form a clay block. Secondly, the elastic tensor of a rotating clay block in the observation coordinate system is calculated by using Bond transformation (Wu, 2006), and the layered shale formed by superposition of multiple rotating clay blocks is approximated by Backus average (Backus, 1962). Then, brittle minerals are mixed up by Hashin-Shtrikman-Hill (HSH) average (Berryman, 1995), and the brittle compound and pores are added to shale matrix with ani-DEM. In addition, cracks are added by Hudson model (Hudson, 1980) to form the rock skeleton. Finally, fluids are added by Brown-Korringa equation (Brown and Korringa, 1975). We find that the modeled and the measured P-wave velocity match well in normal pore pressure segments, but they differ in abnormal segments, which may be caused by abnormal pore pressures. And the modeled P-wave velocity is considered to be the NCT. In order to obtain the interval velocity of the research area, the normalized elastic impedance equation is derived on the basis of the reflectivity approximation for OA medium (Psencik and Martins, 2001; Bachrach, 2014), and then P-wave velocity is inverted by seismic inversion. Thereafter, Eaton’s equation is utilized to predict pore

pressures by using the modeled NCT and the NCT calculated by the conventional method, respectively. The pore pressure prediction results are compared with the measured ones and it is found that the pore pressure from the modeled NCT is more accurate, which verifies the effectiveness of rock physics modeling in pore pressure prediction. 2 Rock physics modeling for orthotropic shale The main composition of shale is clay, and clay is coupled with kerogen (Vernik and Landis, 1996; Wu, 2012; Algazlan et al., 2019). Shale usually owns strong vertical transverse isotropic (VTI) properties due to stratification (Vernik and Nur, 1992; Sayers, 2005). Scanning electron microscopies have displayed laterally distributed cracks in shale (Sondergeld et al., 2010), making shale into OA medium (Guo et al., 2012). Therefore, anisotropic rock physics theory is necessary. The introduced rock physics modeling for orthotropic shale is displayed in Fig. 1.

Fig. 1. The introduced rock physics modeling for orthotropic shale is divided into four steps: (1) clay block simulation; (2) shale stratification simulation; (3) establishment of dry rock skeleton; (4) shale saturation.

2.1 Simulation of elastic stiffness tensor of a clay block The anisotropic SCA model approximates interactions between closely spaced inclusions and can be equivalent to mixed minerals. The anisotropic SCA model is (Hornby, 1994), C %

(

ˆ (C − C = ∑ vn C I + G % % % % % n =1 N

SCA

n

n

SCA

))

−1

−1 N ˆ ( C p − CSCA )    ∑ v p  I + G  % % % %  p =1 

ˆ = 1 (G + G ) G ikjl jkil % 8π

−1

(2)

(3)

where C is the elastic sitffness tensor. v is the volume fraction of the mineral. I % %

ˆ describes the effect of inclusion shape on the equivalent is the identity matrix. G % medium. N is the total number of minerals, and subscripts n, p indicate a certain

mineral, respectively. Lin and Mura (1973) gave the specific expression of Gijkl , see Appendix A. If the aspect ratio is α = 1 , i.e. spherical, the results of ani-SCA would be the same as the isotropic SCA. N

The initial value of CSCA can be given by Voigt upper bound, i.e. K 0 = ∑ f i K i , % i =1 where K0 is the Voigt upper bound, and fi is the volume fraction of a certain mineral. Solving it iteratively, given an accuracy ε , and CSCA is the result of the n % SCA − CSCA would be the final result. nth iteration, if CSCA n n −1 < ε , C n % % % Fig. 2 displays the results of ani-SCA mixing clay and kerogen, and Hashin-Shtrikman (HS) bounds. If the elastic moduli and volume fraction of each component are known, and without the knowledge of geometric details, HS would be the narrowest allowable upper and lower bounds of average mineral moduli (Mavko et al., 1998). The bulk and shear moduli of clay are 22.9GPa and 10.6GPa, respectively. And those of kerogen are 2.9GPa and 2.7GPa , respectively. The aspect ratio of inclusion is α = 1 . When kerogen fraction is about within the range of 0-15%, the equivalent modulus is close to the upper bound of HS, indicating that the duplex minerals are disconnected and the simulation is not good. The ani-SCA simulation results are identical to those from the isotropic SCA models when inclusions are spherical. Generally, the kerogen fraction of organic-rich shale would not exceed 10%, however the ani-SCA cannot accurately simulate shale with less than 15% kerogen fraction. For better simulation, the combination of ani-SCA and ani-DEM is utilized. The anisotropic DEM model is (Hornby, 1994), −1 d 1 ˆ ( C i − C DEM ( v ) )  C DEM ( v ) ) = Ci − C DEM ( v ) )  I + G ( (  % % % dv % % (1 − v ) % %

(4)

where CDEM is the elastic stiffness tensor of host mineral. Ci is the elastic stiffness % % ˆ describes the shape of tensor of inclusion. v is the volume fraction of inclusion. G % inclusion and the expression is the same as that of the anisotropic SCA. Firstly, the ani-SCA model is utilized to mix 50% clay and 50% kerogen, and

then clay is gradually added by the ani-DEM model to reach the measured fraction of clay (Hornby, 1994). As shown in Fig. 3, the equivalent modulus is not close to the HS upper bound when kerogen fraction is below 15%, ensuring the connectivity of clay. And when inclusions are spherical, the results of ani-SCA and ani-DEM are the same as those from the isotropic ones. 10

ani-SCA iso-SCA Hashin Shtrikman bounds

9

K, GPa

8 7 6 5 4 3 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Kerogen fraction

(a)

(b)

Fig. 2. The equivalent results of ani-SCA mixing clay and kerogen. (a) Equivalent bulk modulus; (b) Equivalent shear modulus. 10

ani-SCA+ani-DEM iso-SCA+iso-DEM Hashin Shtrikman bounds

20

ani-SCA+ani-DEM iso-SCA+iso-DEM Hashin Shtrikman bounds

9 8

15 7 6 10 5 4 5 3 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Kerogen fraction

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Kerogen fraction

(a) (b) Fig. 3. 50% clay and 50% kerogen are mixed up by ani-SCA first, and then clay or kerogen is gradually added by anisotropic DEM. (a) Equivalent bulk modulus; (b) Equivalent shear modulus.

2.2 Simulation of shale stratification The layered shale can be assumed to be the superposition of a plurality of rotated clay blocks, and the elastic stiffness tensor of each rotated clay block can be determined by Bond transformation. It is assumed that the polarization angle θ and the azimuth ϕ of a clay block obey a two-dimensional Gaussian distribution,

f (θ , ϕ ) = 1 2πσ 1σ 2

2  −1  θ − µ 2 θ − µ1 )(ϕ − µ2 )  ϕ − µ2    (5) ( 1  exp  +  − 2χ   2 σ σ σ σ χ 2 1 −  ( ) 1− χ 2 1 2  2      1 

where f (θ , ϕ ) is a two-dimensional Gaussian distribution. σ 1 is the standard deviation of polarization angle, controlling the extent to which shale is stratified. µ1

is the mean of polarization angle, controlling the direction of shale stratification. σ 2 is the standard deviation of azimuth which controls the strength of azimuth rotation.

µ2 is the mean of azimuth which controls the direction of the azimuth. χ is the combination tightness parameter, and θ , ϕ are independent of each other when

χ = 0 . In practice, the direction of shale formation is determined by scanning electron microscopy images. In the research area, σ 1 = 30°, µ1 = 0 , σ 2 = 60°, µ2 = 0 , χ = 0 , the joint probability density of polarization angle and azimuth of a clay block is displayed in Fig. 4.

Fig. 4. Joint probability density of polarization angle and azimuth of a clay block.

When observing with seismic waves, the observation coordinate system may be inconsistent with the constitutive coordinate system, and coordinate transformation of the elastic matrix is required. Wu (2006) listed the coordinate rotation transformation of elastic matrix, namely Bond transformation, see Appendix B. The volume fraction of each clay block is,

vk = f (θi , ϕ j ) ∆θ∆ϕ

(6)

∑ ∑ f (θ , ϕ ) ∆θ∆ϕ = 1

(7)

m

n

j =1 i =1

i

j

where f (θ , ϕ ) is a two-dimensional Gaussian distribution; m is the number of azimuth angles, n is the number of polarization angles, k is the total number of clay blocks, k = mn . Then according to Backus average theory (Backus, 1962) (see

Appendix C), the elastic stiffness tensor of stratified shale formed by the superposition of clay blocks can be equated. Backus demonstrated that under the long-wave limit (wavelength was at least 10 times the layer thickness), and without the internal energy consumption caused by friction or liquid viscosity, the fine layered

medium could be approximated as a horizontal transverse isotropic media. The elastic coefficients of the effective medium were averages of algebraic combinations of the original medium’s elastic coefficients. 2.3 Addition of brittle minerals, pores and cracks Brittle minerals like quartz, calcite, dolomite, feldspar and pyrite can be mixed up by Hashin-Shtrikman-Hill (HSH) average (Berryman, 1995), which is an isotropic theory, as shown in Appendix D. The brittle composite are then added to the layered shale by ani-DEM to form the anisotropic shale matrix. The porosity of shale is generally less than 6% (Zou et al., 2013), and oil and gas can be effectively stored in organic, intragranular and intergranular pores (Loucks et al., 2009). In high-or-over-mature shale gas reservoirs, organic pores and micro-cracks are developed, providing sufficient storage space for gas. To simulate rock properties under normal compaction, normal compaction porosity (Athy, 1930) is used in the modeling,

φn = φ0 e− cZ

(8)

where, φn is the normal compaction porosity. φ0 is the porosity of ground surface or sea floor. Z is the depth, m , and c is the constant. In this research area,

φn = 0.15e−0.0007 Z

(9)

Since cracks are the cause of shale’s HTI property, the normal compaction porosity is divided into pores and cracks,

φn = φ pores + φcracks

(10)

φ pores = φclay + φbrittle + φker ogen

(11)

1 φclay = Vshφ ， φbrittle = Vbrittleφ ， φker ogen = Vker ogenφ 2

(12)

Among which,

where Vsh is the clay content. Vker ogen is the kerogen content. Vbrittle is the brittle mineral content. Isolated pores are added by the ani-DEM model, and cracks are then added by Hudson model (Hudson, 1980), see Appendix E. Hudson model is based on thin coin-shaped ellipsoid cracks (crack radius and crack spacing must be much smaller than wavelength, and crack density is assumed in advance and is less than 0.1), and is derived by the scattering theory analysis of average wave field. Hudson model

performs first- and second-order correction on isotropic background. Second-order expansion is not a single convergent sequence, and exceeds the formal scope of Hudson model, its predicted modulus increases with the increase of crack density. Using only first-order correction will give better results than improper use of second-order correction. The isotropic background can be replaced by a VTI background for first- and second-order correction, so as to obtain an OA medium. Since Hudson model assumes isolated cracks, it is suitable for high frequencies. 2.4 Addition of fluids Ani-SCA and ani-DEM both assume isolated inclusions and fluids cannot flow, so the increase in pore pressure caused by wave propagation cannot be balanced, making the models only suitable for high-frequency laboratory conditions. Under low frequencies, the equivalent modulus of dry rock skeleton should be first sought, then fluids can be added by low frequency fluid substitution theory. Brown and Korringa (1975) derived the anisotropic fluid substitution formula, which is an extension of Gassmann’s, S

( dry ) ijkl

−S

( sat ) ijkl

(S = (S

( dry ) ijαα ( dry )

ααββ

) 0 − Sij0αα )( S kl( dry αα − S klαα ) 0 − Sααββ ) + ( β fl − β0 )φ

(13)

( dry ) ( sat ) is the elastic compliance component of rock skeleton. Sijkl is the where, Sijkl 0 elastic compliance component of fluid-saturated rock. Sijkl is the elastic compliance

component of rock matrix. β fl is the compressibility of pore fluid. β 0 is the 0 compressibility of rock matrix, β 0 = Sααββ . φ is the porosity. Repeated subscripts

indicate traversal range summation, i.e. Sij(αα ) = Sij( 11 ) + Sij( 22 ) + Sij( 33 ) , dry

dry

dry

dry

( i, j = 1, 2, 3) ,

( dry ) ( dry ) ( dry ) ( dry ) ( dry ) ( dry ) ( dry ) ( dry ) ( dry ) ( dry ) Sααββ = S1111 + S1122 + S1133 + S 2211 + S 2222 + S 2233 + S3311 + S3322 + S3333

The model assumes the wave-induced increase in pore pressure has sufficient time to be balanced by fluid flow, so it is suitable for low-frequency seismic conditions. According to gas-bearing interpretation, gas and brine are added to gas-bearing segments and gas-free segments, respectively. P-wave velocity of OA medium propagating along Z-axis in both XZ plane and YZ plane is,

VP ( 0° ) =

c33

ρ

(14)

where ρ is the density, kg m3 . 3 Modeling example application and pore pressure prediction for single well 3.1 Aspect ratio inversion of pores and cracks The research area is in Sichuan Basin, and a shale gas well is selected. The structure of this area has experienced multi-stage tectonic movement. Longmaxi Formation shale has a high degree of thermal evolution, with Ro (vitrinite reflectance) greater than 2.2%, and the lower part develops high-quality shale strata with 35~45m thickness and TOC greater than 2%. The reservoir is overpressured and the pore pressure coefficient can reach 1.55. Shale gas production and pressure are stable. The pore aspect ratio is changed from 0.01 to 1, the interval is 0.01, and the crack aspect ratio is changed from 1 to 10, the interval is 0.1. Then for each combination of aspect ratio of pores and cracks (α pores , α cracks ) , the shale rock physics model is established in accordance with the introduced modeling, and P-wave velocity in the vertical direction is calculated by Equation (14). The inversion objective function is:

F = VP ( 0° ) − VPW → min

(15)

where VPW is the logging P-wave velocity, m s . Solving it iteratively, and (α pores , α cracks ) which minimums F is the inversion result. The inverted pore aspect ratio is shown in Fig. 5 as “Pore aspect ratio”, and the inverted crack aspect ratio is as “Crack aspect ratio”. Pore aspect ratio 0.2

0.4

0.6

Crack aspect ratio

0.8

2

1600

1600

1800

1800

2000

2000

2200

2200

2400

2400

2600

2600

2800

2800

3000

3000

3200

3200

3400

3400

4

6

8

10

Fig. 5. The inversion results of aspect ratios.

3.2 Rock physics modeling example application The component parameters used for shale modeling in this research area are

listed in Table 1. The inverted aspect ratio of pores and cracks (Fig. 5) are utilized. The modeled P-wave velocity of each step in the introduced modeling is displayed in Fig. 6. There are differences between the modeled and the measured P-wave velocity within the depth range of 2900-3600m, while the simulations above 2900m are accurate. Since the normal compaction porosity (Equation (9)) is utilized, and the bulk moduli of brine and gas are approximately the values at standard atmospheric pressure, i.e. K brine = 2.2GPa , K gas = 0.3GPa . Therefore, in addition to the inherent error of rock physics models, the pore pressure should also be taken into account, especially in hydrocarbon-bearing regions. In the “Gas-bearing & TOC” panel, the background red indicates containing Class I gas, the yellow indicates containing Class II gas, the green indicates containing Class III gas, and the light brown indicates no gas-bearing. The three categories are mainly divided by TOC (total organic carbon), other criteria include brittleness and DHSR (Differential Horizontal Stress Ratio). TOC is 0-3% for Class III, 3-4% for Class II, and greater than 4% for Class I. In segments below 3400 m, the differences are mainly caused by overpressures due to hydrocarbon generation. Table 1. Bulk and shear moduli of each component in rock physics modeling. Clay

Quartz

Calcite

Dolomite

Feldspar

Pyrite

Kerogen

gas

Brine

K / ( GPa )

31.4

37.0

76.8

76.4

61.0

147.4

2.9

0.3

2.2

µ / ( GPa )

17.0

44.0

32.0

49.7

47.0

132.5

2.7

0

0

Fig. 6. The vertical P-wave velocity VP ( 0° ) obtained at each step of the modeling. The solid blue curves are well logs, and the red dashed curves are the modeled P-wave velocities. In “Gas-bearing & TOC”, the solid black curve is TOC.

3.3 Pore pressure prediction for single well The commonly used pore pressure prediction method Eaton’s equation is,  ∆t  p = σ v − (σ v − σ n )  n   ∆t s 

c

(16)

where p is the pore pressure, MPa . σ v is the overburden stress, MPa . σ n is the hydrostatic pressure, MPa . ∆tn is the acoustic transit time under normal compaction, µ s m . ∆ts is the measured acoustic transit time, µ s m , and the coefficient c is a constant. ∆tn can only be inferred. The commonly used calculation method of the normal compaction trend is (van Ruth et al., 2004; Tingay et al., 2009; Zhang, 2011), ∆tn = ∆tm + ( ∆tml − ∆tm ) e − aZ

(17)

where ∆tm is the acoustic transit time of shale matrix, µ s m . ∆tml is the acoustic transit time of sea floor or ground surface, µ s m . Z is the depth, m , and the constant a is related to regions. The overburden stress σ v is (Yin et al., 2018a),

σ v = ∫ ρ ( Z )gdZ z

0

(18)

where g is the acceleration of gravity, m s 2 . ρ ( Z ) is the denstiy, kg m3 , and Z is the depth, m . The hydrostatic pressure is,

σ n = ρ w gZ

(19)

where ρ w is the density of brine, kg m3 . The vertical P-wave velocity VP ( 0° ) simulated by the orthotropic shale rock physics modeling is regarded as the NCT (the red dashed curve in “Add fluids” of Fig.

6), and pore pressure prediction is performed by Eaton’s equation. The prediction results are displayed in Fig. 7. In the research area, the Eaton exponent is c = 5 . Compared with the measured pore pressures (converted from the pore pressure coefficient collected by repeat formation test (RFT), which is 1.323), the pore pressure predicted by Eaton’s equation using the modeled P-wave velocity is more accurate than that using the conventional method (Equation (17)).

Fig. 7. The “Clay” panel is the clay content; In the “Modeling Vp(0º)” panel, the blue curve is the well log, and the red dashed curve is the modeled P-wave velocity; In the “Pore pressure” panel, the solid blue curve is the pore pressure predicted by Eaton’s equation using the modeled P-wave velocity, while the green dotted curve is the pore pressure predicted by Eaton’s equation using the NCT calculated by the conventional method (Equation (17)), the solid black curve is the hydrostatic pressure, and the red circle is the measured pore pressure. In “Gas-bearing & TOC”, the solid black curve is TOC.

4 Inversion of P-wave velocity in orthotropic medium 4.1 Reflectivity approximation of orthotropic medium and its sensitivity analysis Psencik and Martins (2001) derived the approximation of P-wave incidence and P-wave reflection coefficient in OA medium, and Bachrach (2014) gave the linearized form of the approximation of P-wave reflection coefficient in OA medium, iso ani RPP (ϕ , θ ) = RPP + RPP

iso RPP =

1 ∆α β 2 2 ∆β 1  β 2 2  ∆ρ 2 θ θ 1 + tan − 4 sin + 1 − 4 sin θ  ( ) α α2  β 2 α2 2  ρ ani RPP = ( b1 cos 2 ϕ + b2 sin 2 ϕ ) sin 2 θ +

( b cos 3

4

ϕ + b4 sin 4 ϕ + b5 cos 2 ϕ sin 2 ϕ ) sin 2 θ tan 2 θ

where, b1 = ∆Γ x = ∆δ x − 8 b3 =

β2 β2 ， ∆ b = ∆Γ = ∆ − 8 ∆γ γ δ 2 y y α2 x α2 y

∆ε ∆ε x ∆δ ， b4 = y ， b5 = z 2 2 2

(20) (21)

(22)

A13 + 2 A55 − α 2 A23 + 2 A44 − α 2 A12 + 2 A66 − α 2 ，δ y = ，δz = δx = 2α 2 2α 2 2α 2

γx =

A55 − β 2 A44 − β 2 ， γ = y 2β 2 2β 2

εx =

A11 − α 2 A22 − α 2 ， = ε y 2α 2 2α 2

θ is the incident angle. ϕ is the azimuth. α is the P-wave velocity, m s . β is the S-wave velocity, m s , and Aij = Cij ρ is the density normalized elastic stiffness component. The contribution of each variable to the reflectivity is different, and the higher the contribution, the more accurate the inversion of the variable. Therefore, sensitivity analysis should be carried out for each variable. Let Rα = ∆α α , Rβ = ∆β β , and Rρ = ∆ρ ρ , and the model parameters are listed in Table 2. Table 2. Model parameters for OA medium.

α / ( m ⋅ s −1 )

β / ( m ⋅ s −1 )

ρ / ( g ⋅ cm−3 )

3000

2000

2.4

The upper The lower

2 Rα + 1 vP1 2 Rα − 1

2 Rβ + 1 2 Rβ − 1

vS1

2 Rρ + 1 2 Rρ − 1

b1

b2

b3

b4

b5

−0.4 : 0.08 : 0.4

−0.4 : 0.08 : 0.4

−0.4 : 0.08 : 0.4

−0.4 : 0.08 : 0.4

−0.4 : 0.08 : 0.4

ρ1

Rα is changed from -0.4 to 0.4 with an interval of 0.08, and the other four

variables are kept unchanged at 0.4. Analyze the sensitivity of the reflectivity to Rα . Then, Rβ is changed from -0.4 to 0.4 with an interval of 0.08, and the other four variables are kept unchanged at 0.4. Analyze the sensitivity of the reflectivity to Rβ . The remaining variables are analyzed in the same way, and the analysis results are displayed in Fig. 8. It can be seen that P-wave velocity contributes the most to the reflectivity, so theoretically the inversion result of P-wave velocity would be the most accurate. Fig. 9 is the two-dimensional plan view in terms of the incident angle of Fig. 8. When incident angle θ < 30° , the effect of b3 , b4 , b5 on the reflectivity is approximately zero, so it can be omitted, then,

1  ∆ρ ∆α   1 ∆α β 2 ∆β β 2 ∆ρ  2 + + − 4 − 2 RPP (ϕ ,θ ) =    sin θ α   2 α α2 β α2 ρ  2 ρ + ( b1 cos ϕ + b2 sin ϕ ) sin θ , 2

2

2

(23)

θ < 30 °

Prestack seismic inversion for P-wave velocity can be performed based on this equation.

Fig. 8. Sensitivity analysis of variables of the approximation of P-wave incidence and P-wave reflection coefficient in OA medium. Each graph is an analysis result in which one variable is changed from -0.4 to 0.4 with an interval of 0.08, and the other variables are kept constant at 0.4.

Fig. 9. Sensitivity analysis of variables of the approximation of P-wave incidence and P-wave reflection coefficient in OA medium, this figure is a two-dimensional plan view in terms of the incident angle of Fig. 8.

4.2 Extraction of anisotropic gradient parameters The differentials of anisotropic gradient parameters b1 , b2 in Equation (23) cannot be measured so they need to be estimated from borehole-side trace on profile. Equation (23) can be written as (Chen et al., 2014),

RPP (ϕ , θ ) = P + G sin 2 θ

θ < 30 °

(24)

where the intercept and gradient are: 1  ∆ρ ∆α  P=  + 2 ρ α   1 ∆α β 2 ∆β β 2 ∆ρ  iso ani 2 2 G= −4 2 −2 2  + ( ∆Γ x cos ϕ + ∆Γ y sin ϕ ) = G + G α β α ρ  2 α

If the gradients of three azimuthal angle stacks are known as,

G (ϕ1 ) = G iso + ( ∆Γ x cos 2 ϕ1 + ∆Γ y sin 2 ϕ1 ) G (ϕ 2 ) = G iso + ( ∆Γ x cos 2 ϕ 2 + ∆Γ y sin 2 ϕ2 ) G (ϕ3 ) = G iso + ( ∆Γ x cos 2 ϕ3 + ∆Γ y sin 2 ϕ3 ) two-two subtraction,

B1 = G (ϕ2 ) − G (ϕ1 ) = ∆Γ x ( cos 2 ϕ 2 − cos 2 ϕ1 ) + ∆Γ y ( sin 2 ϕ2 − sin 2 ϕ1 ) B2 = G (ϕ3 ) − G (ϕ 2 ) = ∆Γ x ( cos 2 ϕ3 − cos 2 ϕ2 ) + ∆Γ y ( sin 2 ϕ3 − sin 2 ϕ 2 ) then the differentials of anisotropic gradient parameters in the plane plane

( x1 , x3 )

∆Γ x =

∆Γ y =

( cos (sin

2

2

( x2 , x3 )

and

are, B1 ( sin 2 ϕ3 − sin 2 ϕ 2 ) − B2 ( sin 2 ϕ2 − sin 2 ϕ1 )

(25)

B2 ( cos 2 ϕ 2 − cos 2 ϕ1 ) − B1 ( cos 2 ϕ3 − cos 2 ϕ2 )

(26)

ϕ 2 − cos 2 ϕ1 )( sin 2 ϕ3 − sin 2 ϕ 2 ) − ( cos 2 ϕ3 − cos 2 ϕ2 )( sin 2 ϕ2 − sin 2 ϕ1 ) ϕ3 − sin 2 ϕ2 )( cos 2 ϕ2 − cos 2 ϕ1 ) − ( sin 2 ϕ2 − sin 2 ϕ1 )( cos 2 ϕ3 − cos 2 ϕ2 )

Fig. 10 displays an angle stack with the azimuth of 45° and the incident angle of 8° in the research area of Sichuan basin.

Fig. 10. An angle stack with the azimuth of 45° and the incident angle of 8°.

The AVO intercept and gradient of each angle stack obtained by linear fitting in R (θ ) − sin 2 θ coordinate system according to Equation (24) are displayed in Fig. 11.

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 11. AVO intercepts and gradients of angle stacks with azimuth of 75º, 45º, and 165º. Azimuth-x-P indicates the intercept of the angle stack with azimuth of x; Azimuth-x-G indicates the gradient of the angle stack with azimuth of x.

The differentials of anisotropic gradient parameters calculated by Equation (25) and (26) are displayed in Fig. 12. The differentials of anisotropic gradient parameters

∆Γ x , ∆Γ y and the recursively calculated anisotropic gradient parameters Γ x , Γ y on the borehole-side trace are displayed in Fig. 13.

(a)

(b)

Fig. 12. The differentials of anisotropic gradient parameters in the △Γx -0.1

0

Γx

△Γy 0.1

-0.1

0

( x2 , x3 ) and ( x1 , x3 ) planes.

0.1

-1 -0.5

Γy 0

-1 -0.5

1.3

1.3

1.3

1.3

1.35

1.35

1.35

1.35

1.4

1.4

1.4

1.4

1.45

1.45

1.45

1.45

1.5

1.5

1.5

1.5

1.55

1.55

1.55

1.55

1.6

1.6

1.6

1.6

1.65

1.65

1.65

1.65

0

Fig. 13. The differentials of anisotropic gradient parameters ∆Γ x , ∆Γ y and the recursively

calculated anisotropic gradient parameters Γ x , Γ y on the borehole-side trace.

4.3 Elastic impedance inversion for P-wave velocity The elastic impedance form of Equation (23) is, EI (θ , ϕ ) = α

−8

β2 2 sin θ α2

1− 4

β2 2 sin θ 2 Γ cos 2 ϕ +Γ sin 2 ϕ sin 2 θ x y α2

β ρ The normalized form of the above formula is, 1+ sin 2 θ

1+ sin 2 θ

α  EI (θ , ϕ ) = α 0 ρ 0    α0  P, S-wave velocities

−8

β2

sin 2 θ

e

1− 4

β2

(

)

, θ < 30°

(27)

sin 2 θ

 β  α2  ρ  α2 2( Γ cos 2 ϕ +Γ y sin 2 ϕ ) sin 2 θ e x , θ < 30° (28)     β ρ  0  0 and density well logs α , β , ρ and the extracted

borehole-side anisotropic gradient parameters Γ x , Γ y are utilized to calculate the elastic impedances with different azimuths according to Equation (28), and then azimuthal low frequency models (Yin et al., 2018b; Zong et al., 2018) are established. The elastic impedance sections obtained by poststack constrained sparse spike inversion (Zong et al., 2013; Yin et al. 2014) are displayed in Fig. 14.

(a) EI-45º-8º

(c) EI-45º-16º

(b) EI-165º-8º

(d) EI-165º-16º

(e) EI-75º-16º

Fig. 14. Azimuthal elastic impedances obtained by poststack constrained sparse spike inversion. EI-x-y, where x represents the azimuth and y represents the incident angle: (a) EI-45º-8º, (b) EI-165º-8º, (c) EI-45º-16º, (d) EI-165º-16º, (e) EI-75º-16º.

Equation (28) can be written as Ax = B ,  EI (θ1 , ϕ1 , t1 )  ln EI 0   EI (θ , ϕ , t ) 1 1 2  ln EI 0   M   EI (θ1 , ϕ1 , tn ) ln EI 0 

ln

EI (θ 2 , ϕ1 , t1 )

ln

EI (θ 2 , ϕ1 , t2 )

EI 0 EI 0 M

ln

ln

EI (θ1 , ϕ 2 , t1 )

ln

EI (θ1 ,ϕ 2 , t2 )

EI 0

EI 0

EI (θ 2 , ϕ2 , t1 )

ln

EI (θ 2 ,ϕ 2 , t2 )

M

EI (θ 2 ,ϕ1 , tn ) EI (θ1 , ϕ2 , tn ) ln EI 0 EI 0

let x1 = [ a b c d

ln

EI 0

EI 0 M

ln

EI (θ 2 ,ϕ 2 , tn ) EI 0

EI (θ 2 , ϕ3 , t1 )   α ( t1 )    ln  α0  EI 0  a    EI (θ 2 , ϕ3 , t2 )   b   α ( t2 )    ln  ln α0  EI 0  c =     d   M    M  EI (θ 2 , ϕ3 , t3 )   e   α ( tn )  ln  ln  α 0  EI 0   ln

(29)

e ] . The borehole-side elastic impedances on azimuthal elastic

impedance sections and the logging P, S-wave velocities, density and anisotropic gradient parameters are utilized to solve for x1 , x 2 , x3 , x 4 , x5 corresponding to

ln

α

α0

, ln

β

β0

, ln

ρ

ρ0

, Γ x , Γ y , respectively, by least-square fitting. Then five azimuthal

elastic impedance sections are used as inputs, and sections of P, S-wave velocities, density and anisotropic gradient parameters can be inverted. The inversion results are displayed in Fig. 15. The comparisons of borehole-side inversion results and well logs are displayed in Fig. 16. P-wave velocity inversion result is accurate, which confirms the conclusion that the reflectivity approximation of OA medium is the most sensitive to P-wave velocity.

(a) P-wave velocity

(b) S-wave velocity

(c) Density

(d) Γ

(e) Γ Fig. 15. P, S-wave velocities, density and anisotropic gradient parameters obtained by least-square fitting inversion. Vp,m/s 4500 5500

2500

Vs,m/s

Density,kg/m 3

3000

2550 2650 2750

Γx -1

-0.5

Γy 0

-1

1.4

1.4

1.4

1.4

1.4

1.45

1.45

1.45

1.45

1.45

1.5

1.5

1.5

1.5

1.5

1.55

1.55

1.55

1.55

1.55

1.6

1.6

1.6

1.6

1.6

1.65

1.65

1.65

1.65

1.65

1.7

1.7

1.7

1.7

1.7

-0.5

0

Fig. 16. Comparisons of borehole-side inversion results and well logs. The blue curves are well logs and the dashed red curves are inversion results.

5 Pore pressure prediction example on profile The modeled P-wave velocity (the dashed red curve in panel “Add fluids” in Fig. 6) is used to construct the model by interpolation method of “Inverse Distance” (The weight is inversely proportional to the distance from the sampling point to the known point. In the example section, there is only one well, so weights of sampling points are 1, i.e., the values of sampling points are the same as the known points), representing the normal compaction trend. The P-wave velocity inversion result of OA medium (Fig. 15 (a)) serves as the observed P-wave velocity. Pore pressures are predicted by Eaton’s equation with c = 5 , and the predicted results are displayed in Fig. 17. Single-channel contrast is displayed in Fig. 18.

Fig. 17. Pore pressure prediction results. The black curve is the pore pressure calculated by Eaton’s equation using well logs (the blue curve in the “Pore pressure” of Fig. 7), the white indicates no gas-bearing, the green indicates Class III gas, the yellow indicates Class II gas, and the red indicates Class I gas.

Fig. 18. Single-channel comparison of pore pressure prediction results. The blue curve in the “Pore pressures” panel is the pore pressure predicted by Eaton’s equation using well logs, the dashed red curve is the borehole-side pore pressure in the predicted pore pressure profile (Fig. 17), the green curve is the hydrostatic pressure, and the red circle is the measured pore pressure. The black curve in “Gas-bearing & TOC” is TOC.

6 Discussion The introduced organic-rich shale model (Fig. 1) is equivalent to adding immature kerogen, because the model is designed to estimate the NCT, and the porosity added in the third step is normal compaction porosity which does not include the pores generated on kerogen surface due to maturity. Alternatively, the measured total porosity can be added in the third step to study mature kerogen. However, a new model is required for the study of overmature kerogen. Because overmature kerogen resembles the suspended matter in pores and cannot be regarded as part of the matrix (Zhao et al., 2016), it can be added by solid substitution equation (Ciz and Shapiro, 2007). In addition, the pores are not classified detailedly, the porosity is only divided into pores and cracks, and the fraction of cracks are half of the brittle-related pores (Equation (12)). Furthermore, nanoscale pores are not taken into account either, which

are abundantly present in shale (Zou et al., 2013). Another difficulty in modeling is how to determine the bulk modulus of gas. Under reservoir conditions, the bulk modulus of gas is non-negligible (Mavko et al., 1998). In order to estimate the NCT, the brine and gas at standard atmospheric pressure are added, i.e. K brine = 2.2GPa , K gas = 0.3GPa , which needs further study. Since the bulk and shear moduli of the added clay ( K clay = 31.4GPa , µclay = 17.0GPa ) are larger than the general values ( K clay < 30.0GPa , µclay < 10.0GPa ), it can be inferred that the bulk moduli of the added brine and gas are too small, but this may also be caused by inaccurate logging of brittle mineral fraction. Besides, how to approximate plenty of adsorbed gas is not included in the model. The commonly used Eaton’s equation quantifies the vertical effective stress by the difference between the normal velocity trend (NVT) and the measured velocity. In this study, the rock physics modeling result is utilized as NVT. However, the aim of rock physics modeling is to accurately simulate the elastic properties of rocks, which seems to be contrary to the pore pressure prediction. Recall modeling results (Fig. 7), the modeled and the measured P-wave velocity match well in normal pore pressure segments, and they differ in abnormal segments, which can be explained to be caused by abnormal pore pressures. Conversely, if pore pressures are incorporated in the modeling, the modeled P-wave velocity should be consistent with the measured theoretically. In sediments where overpressures are caused by disequilibrium compaction, the total porosity is abnormally high, and the original Eaton exponent is 3. Eaton’s equation can be flexibly applied to different overpressure mechanism sediments, because the difference between NVT and the measured P-wave velocity can be amplified by increasing Eaton exponent (Bowers, 1995; Tingay et al., 2009). In fluid expansion overpressures, the total porosity increases little (Zhao et al., 2018). If P-wave velocity is simulated with the normal compaction porosity, the difference between the NVT and the measured would be small, and the Eaton exponent should be amplified. In this research area, the genetic mechanism of overpressures is hydrocarbon generation. We use the normal compaction porosity in the modeling, and an amplified Eaton exponent of 5 is calibrated, which is the same as the above conclusion. Therefore, it is necessary to understand the mechanism of overpressures before pore pressure prediction (Swarbrick et al., 2002; Tingay et al., 2009).

Conclusions This study mainly includes anisotropic shale rock physics modeling and the inversion of P-wave velocity from azimuthal elastic impedances in OA medium. The single-well test has been performed to verify that the rock physics modeling can improve the accuracy of pore pressures. And the practicability of predicting pore pressures using P-wave velocity inversion result and the modeled P-wave velocity in OA medium is illustrated through an example in the research area. This study advocates the use of rock physics models to quantify the normal compaction trend. It is found that the difference between the modeled and the measured velocity is due to the lack of consideration of pore pressures. Therefore, pore pressures can be utilized to correct conventional rock physics models to gain more reasonable simulations. Acknowledgements We would like to acknowledge the sponsorship of National Nature Science Foundation of China (U1762103, U1562215), National Grand Project for Science and Technology (2017ZX05036005; 2017ZX05009001; 2017ZX05032003), Science Foundation from SINOPEC Key Laboratory of Geophysics and Young Elite Scientists Sponsorship Program by CAST (2017QNRC001). References Algazlan, M., Pinetown, K., Grigore, M., Chen, Z., Sarmadivaleh, M., Roshan, H., 2019. Petrophysical Assessment of Australian Organic-rich Shales: Beetaloo, Cooper and Perth basins. Journal of Natural Gas Science and Engineering, 102952. https://doi.org/10.1016/j.jngse.2019.102952. Athy, L.F., 1930. Density, porosity, and compaction of sedimentary rocks. AAPG Bulletin, v. 5, p. 1–22. Bachrach, R., 2014. Linearized orthorhombic AVAZ inversion: Theoretical and practical consideration. In: SEG Technical Program Expanded Abstracts 2014 (pp. 528-532). Society of Exploration Geophysicists. https://doi.org/10.1190/segam2014-0767.1. Backus, G.E., 1962. Long-wave elastic anisotropy produced by horizontal layering. Journal of Geophysical Research, 67(11), 4427-4440. https://doi.org/10.1029/JZ067i011p04427. Berryman, J.G., 1995. Mixture theories for rock properties. Rock physics and phase relations: A handbook of physical constants, 3, 205-228. Bowers, G.L., 1995. Pore pressure estimation from velocity data: Accounting for overpressure mechanisms besides undercompaction. SPE Drilling & Completion, 10(02), 89-95. https://doi.org/10.2118/27488-PA. Brown, R.J., Korringa, J., 1975. On the dependence of the elastic properties of a porous rock on the compressibility of the pore fluid. Geophysics, 40(4), 608-616. https://doi.org/10.1190/1.1440551. Budiansky, B., 1965. On the elastic moduli of some heterogeneous materials. Journal of the

Mechanics and Physics of Solids, 13(4), 223-227. https://doi.org/10.1016/0022-5096(65)90011-6. Cao, W., Deng, J., Liu, W., Yu, B., Tan, Q., Yang, L., Li, Y., Gao, J., 2017. Pore pressure and stress distribution analysis around an inclined wellbore in a transversely isotropic formation based on the fully coupled chemo-thermo-poroelastic theory. Journal of Natural Gas Science and Engineering, 40, 24-37. https://doi.org/10.1016/j.jngse.2017.02.002. Chen, H., Zhang, G., Chen, J., Yin, X., 2014. Fracture filling fluids identification using azimuthally elastic impedance based on rock physics. Journal of Applied Geophysics, 110, 98-105. https://doi.org/10.1016/j.jappgeo.2014.09.006. Ciz, R., Shapiro, S.A., 2007. Generalization of Gassmann equations for porous media saturated with a solid material. Geophysics, 72(6), A75-A79. https://doi.org/10.1190/1.2772400. Das, B., Chatterjee, R., 2018. Mapping of pore pressure, in-situ stress and brittleness in unconventional shale reservoir of Krishna-Godavari basin. Journal of Natural Gas Science and Engineering, 50, 74-89. https://doi.org/10.1016/j.jngse.2017.10.021. Dutta, N.C., 2002. Geopressure prediction using seismic data: Current status and the road ahead. Geophysics, 67(6), 2012-2041. https://doi.org/10.1190/1.1527101. Eaton, B.A., 1972. The effect of overburden stress on geopressure prediction from well logs. Journal of Petroleum Technology, 24(08), 929-934. https://doi.org/10.2118/3719-PA. Eaton, B.A., 1975. The equation for geopressure prediction from well logs. In: Fall meeting of the Society of Petroleum Engineers of AIME. Society of Petroleum Engineers. https://doi.org/10.2118/5544-MS. Gor, G. Y., Gurevich, B., 2018. Gassmann theory applies to nanoporous media. Geophysical Research Letters, 45(1), 146-155. https://doi.org/10.1002/2017GL075321. Guo, Z., Chapman, M., Li, X., 2012. A shale rock physics model and its application in the prediction of brittleness index, mineralogy, and porosity of the Barnett Shale. In: SEG technical program expanded abstracts 2012 (pp. 1-5). Society of Exploration Geophysicists. https://doi.org/10.1190/segam2012-0777.1. Gutierrez, M.A., Braunsdor, N.R., Couzens, B.A., 2006. Calibration and ranking of pore-pressure prediction models. The Leading Edge, 25(12), 1516-1523. https://doi.org/10.1190/1.2405337. Han, T., Pervukhina, M., Ben Clennell, M., Neil Dewhurst, D., 2017. Model-based pore-pressure prediction in shales: An example from the Gulf of Mexico, North America. Geophysics, 82(3), M37-M42. https://doi.org/10.1190/geo2016-0504.1. Hornby, B.E., Schwartz, L.M., Hudson, J.A., 1994. Anisotropic effective-medium modeling of the elastic properties of shales. Geophysics, 59(10), 1570-1583. https://doi.org/10.1190/1.1443546. Hottmann, C.E., Johnson, R.K., 1965. Estimation of formation pressures from log-derived shale properties. Journal of Petroleum Technology, 17(06), 717-722. https://doi.org/10.2118/1110-PA. Hudson, J.A., 1980. Overall properties of a cracked solid. In: Mathematical Proceedings of the Cambridge Philosophical Society (Vol. 88, No. 2, pp. 371-384). Cambridge University Press. https://doi.org/10.1017/S0305004100057674. Jiang, M., Spikes, K.T., 2013. Estimation of reservoir properties of the Haynesville Shale by using rock-physics modelling and grid searching. Geophysical Journal International, 195(1), 315-329. https://doi.org/10.1093/gji/ggt250.

Kobchenko, M., Panahi, H., Renard, F., et al., 2011. 4D imaging of fracturing in organic-rich shales during heating. Journal of Geophysical Research: Solid Earth, 116(B12). https://doi.org/10.1029/2011JB008565. Lei, T., Yin, X.Y., Zong, Z.Y., 2018. Pore pressure prediction methods using normal compaction trend based on seismic inversion. Annals of Geophysics, 61(4), 442. https://doi.org/10.4401/ag-7650. Lin, S.C., Mura, T., 1973. Elastic fields of inclusions in anisotropic media (II). Physica status solidi (a), 15(1), 281-285. https://doi.org/10.1002/pssa.2210150131. López, J.L., Rappold, P.M., Ugueto, G.A., Wieseneck, J.B., Vu, C.K., 2004. Integrated shared earth model: 3D pore-pressure prediction and uncertainty analysis. The Leading Edge, 23(1), 52-59. https://doi.org/10.1190/1.1645455. Loucks, R.G., Reed, R.M., Ruppel, S.C., Jarvie, D.M., 2009. Morphology, genesis, and distribution of nanometer-scale pores in siliceous mudstones of the Mississippian Barnett Shale. Journal of sedimentary research, 79(12), 848-861. https://doi.org/10.2110/jsr.2009.092. Maity, D., Ciezobka, J, 2019. Using microseismic frequency-magnitude distributions from hydraulic fracturing as an incremental tool for fracture completion diagnostics. Journal of Petroleum Science and Engineering, 176, 1135-1151. https://doi.org/10.1016/j.petrol.2019.01.111. Mavko, G., Mukerji, T., Dvorkin, J., 1998. The rock physics handbook: Tools for seismic analysis in porous media. University of Cambridge. Miller, R.E., Shenoy, V.B., 2000. Size-dependent elastic properties of nanosized structural elements. Nanotechnology, 11(3), 139. Norris, A.N., 1985. A differential scheme for the effective moduli of composites. Mechanics of materials, 4(1), 1-16. https://doi.org/10.1016/0167-6636(85)90002-X. Pervukhina, M., P. Golodoniuc, B. Gurevich, M.B. Clennell, D.N. Dewhurst, H.M. Nordgård-Bolås, 2015. Prediction of sonic velocities in shale from porosity and clay fraction obtained from logs: A North Sea well case study. Geophysics, 80(1), 1JF-Z39. http://doi.org/10.1190/geo2014-0044.1. Pšenčík, I., Martins, J.L., 2001. Properties of weak contrast PP reflection/transmission coefficients for weakly anisotropic elastic media. Studia Geophysica et Geodaetica, 45(2), 176-199. https://doi.org/10.1023/A:102186832. Sayers, C.M., 2005. Seismic anisotropy of shales. Geophysical prospecting, 53(5), 667-676. https://doi.org/10.1111/j.1365-2478.2005.00495.x. Sayers, C.M., Johnson, G.M., Denyer, G., 2002. Predrill pore-pressure prediction using seismic data. Geophysics 67(4), 1286-1292. https://doi.org/10.1190/1.1500391. Serebryakov, V.A., Robertson Jr, J.O., Chilingarian, G.V. (Eds.), 2002. Origin and prediction of abnormal formation pressures (Vol. 50). Gulf Professional Publishing. Sondergeld, C.H., Ambrose, R.J., Rai, C.S., Moncrieff, J., 2010. Micro-structural studies of gas shales. In: SPE unconventional gas conference. Society of Petroleum Engineers. https://doi.org/10.2118/131771-MS. Song, C., Lu, Y., Tang, H., Jia, Y., 2016. A method for hydrofracture propagation control based on non-uniform pore pressure field. Journal of Natural Gas Science and Engineering, 33, 287-295. https://doi.org/10.1016/j.jngse.2016.05.029. Swarbrick, R., 2012. Review of pore-pressure prediction challenges in high-temperature areas.

The Leading Edge, 31(11), 1288-1294. https://doi.org/10.1190/tle31111288.1. Swarbrick, R.E., Osborne, M.J., Yardley, G.S., 2002. Comparison of overpressure magnitude resulting from the main generating mechanisms. In: Huffman, A.R., Bowers, G.L. (Eds.), Pressure Regimes in Sedimentary Basins and their Prediction: AAPG Memoir, 76, pp. 1–12. Terzaghi, K., Peck, R.B., Mesri, G., 1996. Soil mechanics in engineering practice. John Wiley & Sons. Thomsen, L., 1986. Weak elastic anisotropy. Geophysics, 51(10), 1954-1966. https://doi.org/10.1190/1.1442051. Tingay, M.R., Hillis, R.R., Swarbrick, R.E., Morley, C.K., Damit, A.R., 2009. Origin of overpressure and pore-pressure prediction in the Baram province, Brunei. Aapg Bulletin, 93(1), 51-74. https://doi.org/10.1306/08080808016. van Ruth, P., Hillis, R., Tingate, P., 2004. The origin of overpressure in the Carnarvon Basin, Western Australia: Implications for pore pressure prediction. Petroleum Geoscience, 10(3), 247-257. https://doi.org/10.1144/1354-079302-562. Vernik, L., Landis, C., 1996. Elastic anisotropy of source rocks: Implications for hydrocarbon generation and primary migration. AAPG bulletin, 80(4), 531-544. Vernik, L., Nur, A., 1992. Petrophysical analysis of the Cajon Pass scientific well: implications for fluid flow and seismic studies in the continental crust. Journal of Geophysical Research: Solid Earth, 97(B4), 5121-5134. Wang, Z., Wang, R., 2015. Pore pressure prediction using geophysical methods in carbonate reservoirs: Current status, challenges and way ahead. Journal of Natural Gas Science and Engineering, 27, 986-993. https://doi.org/10.1016/j.jngse.2015.09.032. Wu, G.C., 2006. Seismic wave propagation and imaging in anisotropy media. Dongying: China University of Petroleum Press. Wu, X., Chapman, M., Li, X. Y., Dai, H., 2012. Anisotropic elastic modelling for organic shales. In: 74th EAGE Conference and Exhibition incorporating EUROPEC 2012. http://doi.org/10.3997/2214-4609.20148651. Yin X.Y., Liu X.X., 2016. Research status and progress of the seismic rock-physics modeling methods. Geophysical Prospecting for Petroleum 55(3): 309-325. http://doi.org/10.3969/j.issn.1000-1441.2016.03.001. Yin, X.Y., Cao, D.P., Wang, B.L., Zong, Z.Y., 2014. Advances in fluid discrimination methods based on prestack seismic inversion. Oil Geophysical Prospecting 49(1), 22-34. http://doi.org/10.3969/j.issn.1000-0952.2010.15.020. Yin, X.Y., Ma, N., Ma, Z.Q., Zong, Z.Y., 2018a. Review of in-situ stress prediction technology. Geophysical Prospecting for Petroleum 57(4): 488-504. http://doi.org/10.3969/j.issn.1000-1441.2018.04.001. Yin, X.Y., Zhang, H.X., Zong, Z.Y., 2018b. Research status and progress of 5D seismic data interpretation in OVT domain. Geophysical Prospecting for Petroleum 57(2): 155-178. http://doi.org/10.3969/j.issn.1000-1441.2018.02.001. Zhang, J., 2011. Pore pressure prediction from well logs: Methods, modifications, and new approaches. Earth-Science Reviews, 108(1-2), 50-63. https://doi.org/10.1016/j.earscirev.2011.06.001. Zhao, J., Li, J., Xu, Z., 2018. Advances in the origin of overpressures in sedimentary basins. Petroleum Research, 3(1), 1-24. https://doi.org/10.1016/j.ptlrs.2018.03.007.

Zhao, L., Qin, X., Han, D. H., Geng, J., Yang, Z., Cao, H., 2016. Rock-physics modeling for the elastic properties of organic shale at different maturity stages. Geophysics, 81(5), D527-D541. https://doi.org/10.1190/geo2015-0713.1. Zhu, Y., Xu, S., Payne, M., Martinez, A., Liu, E., Harris, C., Bandyopadhyay, K., 2012. Improved rock-physics model for shale gas reservoirs. In: SEG Technical Program Expanded Abstracts 2012 (pp. 1-5). Society of Exploration Geophysicists. https://doi.org/10.1190/segam2012-0927.1. Zimmerman, R.W., 1990. Compressibility of sandstones (Vol. 29). Elsevier. Zong, Z., Wang, Y., Li K., Yin, X., 2018. Broadband Seismic Inversion for Low-Frequency Component of the Model Parameter. IEEE Transactions on Geoscience & Remote Sensing, 99:1-8. http://doi.org/10.1109/TGRS.2018.2810845. Zong, Z., Yin, X., Wu, G., 2013. Direct inversion for a fluid factor and its application in heterogeneous reservoirs. Geophysical Prospecting 61(5), 998-1005. https://doi.org/10.1111/1365-2478.12038. Zou, C. N., Tao, S. Z., Hou, L. H., 2013. Unconventional oil and gas geology, 2nd Edition. Beijing: Geological Publishing House.

Appendix A. The expression of tensor Gijkl Lin and Mura (1973) gave the expression of component M ijkl , M ijkl = Gikjl , defing,

c11 − c12 1 ， f = c44 ， g = c13 + c44 ， h = c33 ， ρ = 2 α where α is the aspect ratio of the inclusion, then the expression of Gijkl is,

d = c11 ， e =

{

}

∆ −1 = e (1 − x 2 ) + f ρ 2 x 2   d (1 − x 2 ) + f ρ 2 x 2   f (1 − x 2 ) + h ρ 2 x 2  − g 2 ρ 2 x 2 (1 − x 2 )

G1111 = G2222 = ∆ (1 − x ) { f (1 − x ) + h ρ 2∫

π

1

2

2

2

0

}

x 2  ( 3e + d ) (1 − x 2 ) + 4 f ρ 2 x 2  − g 2 ρ 2 x 2 (1 − x 2 ) dx

G3333 = 4π ∫ ∆ρ 2 x 2  d (1 − x 2 ) + f ρ 2 x 2   e (1 − x 2 ) + f ρ 2 x 2  dx 1

0

G1122 = G2211 = ∆ (1 − x ) { f (1 − x ) + h ρ 2∫

π

1

2

2

0

2

}

x 2  ( e + 3d ) (1 − x 2 ) + 4 f ρ 2 x 2  − 3 g 2 ρ 2 x 2 (1 − x 2 ) dx

G1133 = G2233 =

{

}

2π ∫ ∆ρ 2 x 2 ( d + e ) (1 − x 2 ) + 2 f ρ 2 x 2   f (1 − x 2 ) + hρ 2 x 2  − g 2 ρ 2 x 2 (1 − x 2 ) dx 0 1

G3311 = G3322 = 2π ∫ ∆ (1 − x 2 )  d (1 − x 2 ) + f ρ 2 x 2  e (1 − x 2 ) + f ρ 2 x 2  dx 1

0

G1212 =

π

∆ (1 − x ) 2∫ 1

0

2 2

{g ρ x 2

2

2

}

− ( d − e )  f (1 − x 2 ) + h ρ 2 x 2  dx

G1313 = G2323 = ( −2π ) ∫ ∆g ρ 2 x 2 (1 − x 2 )  e (1 − x 2 ) + f ρ 2 x 2  dx 1

0

Appendix B. Bond transformation The polarization angle of transverse isotropic medium is θ , and the azimuth angle is ϕ , then,

 cos 2 θ   0  sin 2 θ M θ =  0 1  sin 2θ 2  0

0 1 0 0

sin 2 θ 0 cos 2 θ 0 1 0 − sin 2θ 2 0 0

 cos 2 ϕ  2  sin ϕ  0 M φ =  0  0  1  2 sin 2ϕ

sin 2 ϕ cos 2 ϕ 0 0 0 1 − sin 2ϕ 2

0 0 0 cos θ

− sin 2θ 0 sin 2θ 0

0

cos 2θ

− sin θ

0

0 0 0 0 0 0 1 0 0 0 cos ϕ sin ϕ 0 − sin ϕ cos ϕ 0

0

0

0   0  0   sin θ   0   cos θ  − sin 2ϕ   sin 2ϕ  0   0  0    cos 2ϕ  

Let the elastic stiffness tensor of the transverse isotropic medium in the constitutive coordinate system be C0 , then the elastic stiffness matrix of the polarized azimuthal anisotropic medium in the observation coordinate system has the following expression (Wu, 2006):

 c11 c  12 c Cθφ = M φ M θC0 M θT M φT =  13 c14 c15  c16

c12 c22 c23 c24 c25

c13 c23 c33 c34 c35

c14 c24 c34 c44 c45

c15 c25 c35 c45 c55

c26

c36

c46

c56

c16  c26  c36   c46  c56   c66 

Appendix C. Backus average Each clay block can be equivalent to a VTI medium, and its elastic stiffness tensor can be written as,

a b  f     

b a f

f f c d d

    ,    m

m=

1 ( a − b) 2

Backus (1962) demonstrated that, in the long-wave limit, the equivalent stiffness of a layered medium composed of multi-layer transverse isotropic materials is, A B F B A F  F F C     

    , D   D  M

M=

1 ( A − B) 2

The volume fraction of each clay block is vi , see equation (6), then, −1

2

−1

2

N  N   N  A = ∑ vi ( a − f 2 c −1 ) +  ∑ vi c −1   ∑ vi fc −1  i =1  i =1   i =1 

 N   N  B = ∑ vi ( b − f c ) +  ∑ vi c −1   ∑ vi fc −1  i =1  i =1   i =1  N

2 −1

−1

−1

N  N   N  C =  ∑ vi c −1  , D =  ∑ vi d −1  , M = ∑ vi m i =1  i =1   i =1 

 N  F =  ∑ vi c −1   i =1 

−1 N

∑ v fc i =1

−1

i

Appendix D. Hashin-Shtrikman-Hill (HSH) average If the elastic moduli and volume fraction of each component are known, and without the knowledge of geometric details of how components are combined with each other, Hashin-Shtrikman (HS) bounds would be the narrowest allowable upper and lower bounds of average mineral moduli (Mavko et al., 1998). A more general form of boundaries for more than two compositions is, −1

K

HS +

 N  1 4 =  ∑ vi  − µmax 3  i =1 K i + 4 3 µmax  −1

K

HS −

 N  1 4 =  ∑ vi  − µ min 3  i =1 K i + 4 3 µmin 

−1

µ HS +

  N 1  = ∑ vi  i =1 µ  9 K + 8µmax µi + max  max  6  K max + 2 µ max 

   − µ max  9 K max + 8µmax    6  K max + 2 µ max     

µ HS −

  N 1 =  ∑ vi  i =1 µ  9 K + 8µ min µi + min  min  6  K min + 2 µ min 

   − µmin  9 K min + 8µ min    6  K min + 2 µ min     

−1

where K HS + is the upper bound of the bulk modulus. K HS − is the lower bound.

µ HS + is the upper bound of the shear modulus. µ HS − is the lower bound. K max , µmax are the largest bulk and shear moduli among all mineral components, and K min , µmin are the smallest bulk and shear moduli among all mineral components. N is the total number of components. i indicates a certain mineral component, and vi is the volume fraction of a certain component. Hashin-Shtrikman-Hill (HSH) average is the averaging of HS upper and lower bounds,

M HSH =

M HS + + M HS − 2

where M can be any modulus.

Appendix E. Hudson model for fractured medium Hudson introduced Hudson model based on the scattering theory analysis of the average wavefield in an elastic solid containing thin coin-shaped ellipsoidal cracks or inclusions (Hudson, 1980), the equivalent modulus is, cijeff = cij0 + cij1 + cij2

where cij0 is the modulus of isotropic background. cij1 and cij2 are the first- and second-order corrections, respectively. For a single crack group with a vertical direction along the normal of ground, the fractured medium exhibits transverse isotropic symmetry, and the corrections are,

λ ( λ + 2µ ) ( λ + 2µ ) εU λ2 1 1 c = − ε U 3 ， c13 =− ε U 3 ， c33 =− 3 µ µ µ 2

1 11

1 c144 = − µε U1 ， c66 =0

and, c112 =

q q 2 2 q λ2 2 ( εU 3 ) ， c132 = λ ( εU 3 ) ， c332 = ( λ + 2µ )( εU 3 ) 15 15 15 λ + 2 µ 2 = c44

q = 15

2 µ ( 3λ + 8µ ) 2 ( εU1 ) ， c662 = 0 15 λ + 2 µ

N 3 3φ λ2 λ + 28 + 28 ， ε = a = 2 V 4πα µ µ

where λ , µ are the elastic moduli of the isotropic background. a is the radius of crack. α is the aspect ratio of crack. cij1 , cij2 follows the usual symmetrical properties of transverse isotropy. ε is the crack density. φ is the porosity. U1 ,U 3 depends on the state of cracks, and for dry cracks,

U1 =

16 ( λ + 2µ ) 4 ( λ + 2µ ) ， U3 = 3 ( 3λ + 4µ ) 3(λ + µ )

1. The rock physics modeling of shale gas in orthotropic (OA) medium is designed to equate the normal compaction trend (NCT) to predict pore pressures. 2. P-velocity is inverted from azimuthal seismic data based on the reflectivity approximation of OA medium, and the pore pressure of the research area is predicted using the modeled P-velocity. 3. The prediction result is consistent with the characteristics of overpressure mechanism of hydrocarbon generation.

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: