Fluid Discrimination Based on Inclusion-Based Method for Tight Sandstone Reservoirs

Abstract

Fluid discrimination is challenging for reservoir prediction, especially for tight sandstones with special petrophysical properties. In this paper, we first review the effective medium models that are widely used in seismic exploration and a variety of inversion methods and reservoir prediction strategies in reservoir prediction. Rock physics modeling takes an important role in reservoir prediction by linking petrophysical properties and elastic parameters. We also review the theoretical implications for different rock physics models that are based on the inclusion-based method, focusing specifically on the modeling workflow for conventional sand-shale reservoirs and two models for tight sandstones. The applicability of the conventional fluid substitution equations is analyzed in detail. Then, a new inclusion-based rock physics model for tight sandstones is proposed by considering the fluid pressure ratio between cracks and stiff pores. The proposed model helps to highlight the difference between different pores and present reasonable fluid information. In the application, a detailed prediction process for fluid discrimination is given, in which the Bayes posterior prediction framework is adopted to provide the maximum posterior probability solution and its posterior probability. Field data applications demonstrate the effectiveness of the proposed method.

FormalPara Article Highlights

• The theoretical implications for the rock physics models that are based on the inclusion-based method are reviewed

• A new inclusion-based rock physics model for tight sandstones is proposed by considering the fluid pressure ratio between cracks and stiff pores

• A detailed prediction process of fluid discrimination for tight sandstones is presented

1 Introduction

Rock physics modeling, elastic parameter inversion, and reservoir prediction strategy are essential for reservoir fluid discrimination. Studies of these aspects can improve the reservoir prediction accuracy, including elastic parameter inversion, lithology discrimination, and fluid identification, which can serve for the fine description of reservoirs and finally identify the hydrocarbon-bearing zone.

1.1 Inclusion-Based Rock Physics Models

Fluid discrimination is the focus of reservoir prediction for hydrocarbon exploration. However, in addition to the pore fluid, other components, such as pore structure and solid grains, can also affect the wave propagation, which limits the accuracy of fluid identification (Song et al. 2016a; Guo et al. 2020). Thus, the contribution of different components should be clarified. Rock physics modeling provides a tool to relate petrophysical properties (e.g., water saturation and porosity) to elastic parameters (e.g., P- and S-wave velocities) (Mavko et al. 2020). Without considering the microstructure of inclusions, the upper and lower bounds of the elastic modulus can be predicted, where the volume fractions and the constituent moduli should be specified. The Hashin–Shtrikman bounds and Voigt-Reuss bounds are commonly used to estimate the range of mineral modulus for a mixture of mineral grains or mineral and pore fluid (Mavko et al. 2020). In the rock physics modeling, the Voigt–Reuss–Hill average (Kobayashi and Mavko 2016; Wang et al. 2020a) is preferred to calculate the modulus of rock matrix, and the Wood’s formula (Mavko et al. 2020) is used for the fluid. If the microstructure is considered, the detailed variation in rock modulus can be obtained. Inclusion-based effective medium models (EMMs), such as Kuster–Toksöz (KT) model (Kuster and Toksöz 1974), differential effective medium (DEM) model (Berryman 1992), self-consistent approximation (SCA) model (Wu 1966), and Mori–Tanaka (MT) model (Mori and Tanaka 1973), are commonly applied. In these models, the complex pore structure is described using different pore aspect ratios. Large pore aspect ratios refer to stiff pores or sand pores. In contrast, small pore aspect ratios indicate cracks or shale pores (Xu and White 1996; Ruiz and Dvorkin 2010; Wang et al. 2020c). For sandstones, two types of pores including cracks and stiff pores are considered, while for carbonate rocks, more than three types of pores including interparticle pores, microcracks, and stiff pores are considered (Xu and Payne 2009). All EMMs have a presupposition that the inclusions are assumed to be isolated with respect to fluid flow, which is appropriate for high-frequency laboratory conditions (Mavko et al. 2020). For low-frequency situations, such as seismic data applications, EMMs are used to provide dry inclusions, which are then filled with fluid using Gassmann’s equations (Biot 1956; Song et al. 2016b). At low frequency limit, the pore space is considered fully connected and there is no fluid pressure gradient (Mavko et al. 2020; Wang et al. 2021). The modulus variation falling between high and low frequency limits can be depicted by considering the fluid flow (Borgomano et al. 2017; Pang et al. 2019; Wang et al. 2021). Many researchers have studied the attenuation and dispersion caused by wave-induced fluid flow including squirt flow (Dvorkin and Nur 1993; Chapman et al. 2002, 2003; Tang 2011; Song et al. 2016a), patchy saturation (White 1975; Dutta and Odé 1979; Kobayashi and Mavko 2016) and double porosity (Ba et al. 2016, 2017; Sun et al. 2018; Zhang et al. 2020) theories. These models introduce more parameters, which can improve the accuracy of rock physics modeling, but are not conducive to inverse problems such as reservoir prediction. In seismic inversion, the combination of EMMs and Gassmann’s equations is commonly adopted for lithology discrimination and fluid identification (Russell et al. 2011; Zong et al. 2015; Zhou et al. 2020). For conventional sand-shale reservoirs, this approach is feasible because the reservoir pores are usually connected. Considering the special petrophysical properties of tight sandstone reservoirs, such as low porosity, low permeability, and complex pore structure, the conventional rock physics models show limitations (Wang et al. 2020a). Some researchers have investigated the rock physics modeling of tight reservoirs (Smith et al. 2009; Ruiz and Cheng 2010; Avseth et al. 2014; Ba et al. 2019; Wang et al. 2020c), however, relevant papers are few. In this paper, we review and analyze two rock physics models for tight sandstones and propose a new model by considering the fluid pressure ratio.

1.2 Seismic Elastic Parameter Inversion

In addition to rock physics modeling, seismic inversion is required to link elastic parameters with seismic data for fluid discrimination (Chen and Zhang 2017; Pan et al. 2017). Compared with post-stack seismic inversion which can only give one parameter (e.g., P-wave impedance), pre-stack seismic inversion can provide abundant information including P- and S-wave velocities and density (Zhi et al. 2016; Pan and Zhang 2018; Wang et al. 2020a). Amplitude variation with offset (AVO) inversion is an effective way to obtain elastic parameters from seismic angle gathers (Zong et al. 2012; Pan et al. 2018; Chen 2020). The Zoeppritz equations and their linear approximations are widely applied in seismic exploration (Shuey 1985; Smith and Gidlow 1987; Fatti et al. 1994). These linear equations provide different combinations of elastic parameters with clear physical meaning and can reduce the computational time. However, they are only suitable for near offsets and can lead to biases in the calculation of converted S-wave reflection coefficients. Thus, the exact Zoeppritz equations are used directly to predict elastic parameters (Zhou et al. 2020). Although Zoeppritz equations-based methods can provide reflected P- and S-waves as well as transmitted P- and S-waves, they are limited to a single interface and cannot describe transmission loss and multiples in multiple layers (Mavko et al. 2020). Wang et al. (2020a) pointed out that due to the special petrophysical properties, weak reflection events are more likely to appear in tight sandstones and they can be weakened by the transmission and obscured by multiples. The reflectivity method is a useful way to consider the cumulative effect of multiple layers, which has been studied by many researchers (Fryer 1980; Mallick and Adhikari 2015; Liu et al. 2016) and proven to be valid for tight sandstone reservoirs (Wang et al. 2020a). Thus, the reflectivity method is adopted to predict the P- and S-wave velocities and density in this paper. Prior constraint is an important part of the inversion algorithm, which can be used to enhance inversion stability and reduce multiple solutions. For deterministic inversion, a number of regularization norms are used to improve the inversion accuracy, including l₁ norm (Hu et al. 2012), l₂ norm (Tihonov 1963), hybrid l₁/l₂ norm (Bube and Langan 1997), l_1–2 norm (Esser et al. 2013), and l_p norm (Woodworth and Chartrand 2016). For probabilistic prediction, the prior probability is usually assumed to follow a specific distribution, such as Gaussian distribution (Fjeldstad and Grana 2018), Cauchy distribution (Alemie and Sacchi 2011), and Laplace distribution (Theune et al. 2010). In the reflectivity method-based AVO inversion, the framework of Bayesian inversion is adopted here, where the prior distributions of elastic parameters are assumed to follow the Gaussian distribution.

1.3 Reservoir Prediction Strategies

The process from seismic data to elastic parameters and then to petrophysical properties is indirect reservoir prediction, while the process from seismic data directly to petrophysical properties is named direct reservoir prediction. For direct reservoir prediction, a simplified rock physics model is needed to obtain a direct expression from seismic angle gathers to petrophysical properties (Grana et al. 2012; Chen and Zhang 2017). However, the complex petrophysical properties of tight sandstones bring challenges to linearization. In addition, considering the difficulty of combing the reflectivity method and the rock physics model of tight sandstones, indirect reservoir prediction is preferred in this paper. Bayes posterior prediction framework is adopted as the inversion method, where the Markov chain Monte Carlo (McMC) and Metropolis–Hastings (MH) algorithm (Maschio and Schiozer 2014; Grana 2020) is used to obtain the maximum posterior probability solution and its posterior probability. If fluid-related petrophysical properties are difficult to obtain, sensitivity analysis of elastic parameters can be performed to obtain a fluid factor that is sensitive to the pore fluid. Then, the fluid factor is predicted in an inversion framework to indirectly identify the fluid. This strategy is widely used in seismic inversion (Fawad et al. 2020; Wang et al. 2022). Another approach for the reservoir prediction is to leverage machine learning algorithms including random forests (Silva et al. 2020), hidden Markov models (Lindberg and Omre 2015; Feng et al. 2018; Wang et al. 2020b), deep neural networks (Yu and Ma 2021), and convolutional neural networks (Wu et al. 2020), which avoid the complex rock physics modeling. Training data are critical to machine learning algorithms. Insufficient training data and measurement bias can affect the prediction results significantly, and the transfer ability of machine learning algorithms also limits their application. Machine learning algorithms have both advantages and disadvantages. In order to improve the general applicability of prediction methods, the support of rock physics theory is still needed. Wang et al. (2020c) proposed a joint posterior probability approach for fluid discrimination by combining the contributions of deterministic rock physics and statistical rock physics, which is driven by both data and rock physics model.

With the continuous exploration and development of hydrocarbon resources, as well as the continuous depletion of oil and gas from conventional reservoirs, the study on unconventional reservoirs, such as tight sandstones, has received increasing attention. Fluid discrimination is a key technology for the fine description of reservoirs, of which rock physics modeling is an essential part. In this paper, we first review the modeling workflow of conventional sand-shale reservoirs that is widely used in seismic exploration. Then, a series of rock physics models for tight sandstones are presented, including a newly proposed rock physics model considering fluid pressure ratio. In the application, a detailed prediction process for fluid discrimination is provided, which helps to improve the prediction accuracy. Both well-logging data and seismic data applications are shown to demonstrate the effectiveness of the proposed method.

2 Methods

In this section, we review the modeling workflow of conventional sand-shale reservoirs and analyze the applicability of Gassmann’s fluid substitution equations. Then, two rock physics models for tight sandstones are shown and investigated. Finally, a new inclusion-based model considering the fluid pressure ratio is proposed.

2.1 Conventional Sand-Shale Reservoir Modeling

For hydrocarbon-bearing sandstone reservoirs, pore structure and fluid are critical, which affect wave propagation in rocks. Brittle minerals, such as quartz, are prone to develop cracks during the geologic deposition process, which causes the complexity of the pore structure. To identify the pore fluid with high accuracy, a fine delineation of the pore structure is needed. Some researchers divided the pore space into two types, in which one is considered as the sand pore with a large aspect ratio and the other is the shale pore with a small aspect ratio (Xu and White 1996; Keys and Xu 2002). Others assumed that the pore structures contain cracks and stiff pores, which are two extreme pore structures (Ruiz and Dvorkin 2010). The stiff pore is considered to be spherical, and its aspect ratio is set to 1. To better solve the inverse problem in the reservoir prediction, some researchers simplified the rock physics model by treating the pore structure as a type of aspect ratio equivalently (Wang et al. 2022).

For conventional sand-shale reservoirs (Fig. 1a), induced pore pressures are assumed equilibrated throughout the pore space, especially in the seismic inversion. Thus, Fig. 1a can be equivalently depicted by Fig. 1b. For Fig. 1b, EMMs and Gassmann’ equations are commonly used to establish the connection between elastic parameters and petrophysical properties. Among EMMs, the KT model is commonly used for the conventional sand-shale and carbonate reservoirs (Xu and Payne 2009). However, this does not mean that the KT model is the best and it is not suitable for dilute inclusions. Besides, other EMMs also have limitations. For the DEM model, the adding path of different components should be considered, and for the SCA model, an iterative process is needed. In contrast, the MT scheme is independent of the adding path, takes more effects of inclusions into account, and has an explicit expression (Song et al. 2016b). Guéguen and Kachanov (2011) pointed out that the MT scheme is more suitable for cracked rocks.

Fig. 1

Diagrammatic sketch of real rock for conventional sand-shale reservoirs and its equivalent sketch for rock physics modeling

Here, the KT model is used as an example to introduce the modeling process for conventional sand-shale reservoirs. Figure 2 displays the rock physics modeling workflow. First, the minerals present in the rock are mixed to get the moduli of rock matrix by using the Reuss–Voigt–Hill average (Mavko et al. 2020). Then, different pores are added to form the dry rock frame by using the KT model, where the pore space is filled with no fluid by assuming that the bulk and shear moduli of the fluid equal to 0. The expressions are given by

$$\left( {K_{{\rm m}} - K_{{\rm d}} } \right)\frac{{K_{{\rm m}} + \frac{4}{3}\mu_{{\rm m}} }}{{K_{{\rm d}} + \frac{4}{3}\mu_{{\rm m}} }} = \sum\limits_{i = 1}^{N} {\phi_{i} K_{{\rm m}} P^{{\rm mi}} } ,$$

(1)

$$\left( {\mu_{{\rm m}} - \mu_{{\rm d}} } \right)\frac{{\mu_{{\rm m}} + \frac{{\mu_{{\rm m}} \left( {9K_{{\rm m}} + 8\mu_{{\rm m}} } \right)}}{{6\left( {K_{{\rm m}} + 2\mu_{{\rm m}} } \right)}}}}{{\mu_{{\rm d}} + \frac{{\mu_{{\rm m}} \left( {9K_{{\rm m}} + 8\mu_{{\rm m}} } \right)}}{{6\left( {K_{{\rm m}} + 2\mu_{{\rm m}} } \right)}}}} = \sum\limits_{i = 1}^{N} {\phi_{i} \mu_{{\rm m}} Q^{{\rm mi}} } ,$$

(2)

where K_d and μ_d are the bulk and shear moduli of the dry rock frame. K_m and μ_m are the bulk and shear moduli of the rock matrix. $\phi_{i}$ is the volume content of ith inclusion. N is the total number of inclusions. Coefficients $Q^{{\rm mi}}$ and $P^{{\rm mi}}$ are given in Appendix 1, which are determined by the pore structure and fluid.

Fig. 2

Rock physics modeling workflow for conventional sand-shale reservoir

In Eqs. 1 and 2, the bulk and shear moduli are coupled, which can be solved by iteration. The KT model is limited to dilute inclusions. Thus, different types of pores should be added gradually. This is similar to the DEM model.

The fluid-saturated rocks are finally obtained by using Gassmann’s equations:

$$K_{{\rm Ga}} = K_{{\rm d}} + \frac{{\left( {1 - {{K_{{\rm d}} } \mathord{\left/ {\vphantom {{K_{{\rm d}} } {K_{{\rm m}} }}} \right. \kern-\nulldelimiterspace} {K_{{\rm m}} }}} \right)^{2} }}{{{{\phi^{*} } \mathord{\left/ {\vphantom {{\phi^{*} } {K_{{\rm f}} + {{(1 - \phi^{*} )} \mathord{\left/ {\vphantom {{(1 - \phi^{*} )} {K_{{\rm m}} - {{K_{{\rm d}} } \mathord{\left/ {\vphantom {{K_{{\rm d}} } {K_{{\rm m}}^{2} }}} \right. \kern-\nulldelimiterspace} {K_{{\rm m}}^{2} }}}}} \right. \kern-\nulldelimiterspace} {K_{{\rm m}} - {{K_{{\rm d}} } \mathord{\left/ {\vphantom {{K_{{\rm d}} } {K_{{\rm m}}^{2} }}} \right. \kern-\nulldelimiterspace} {K_{{\rm m}}^{2} }}}}}}} \right. \kern-\nulldelimiterspace} {K_{{\rm f}} + {{(1 - \phi^{*} )} \mathord{\left/ {\vphantom {{(1 - \phi^{*} )} {K_{{\rm m}} - {{K_{{\rm d}} } \mathord{\left/ {\vphantom {{K_{{\rm d}} } {K_{{\rm m}}^{2} }}} \right. \kern-\nulldelimiterspace} {K_{{\rm m}}^{2} }}}}} \right. \kern-\nulldelimiterspace} {K_{{\rm m}} - {{K_{{\rm d}} } \mathord{\left/ {\vphantom {{K_{{\rm d}} } {K_{{\rm m}}^{2} }}} \right. \kern-\nulldelimiterspace} {K_{{\rm m}}^{2} }}}}}}}},$$

(3)

$$\mu_{{\rm Ga}} = \mu_{{\rm d}} ,$$

(4)

where $\phi^{*} = \sum\limits_{i = 1}^{N} {\phi_{i} }$. K_Ga and $\mu_{{\rm Ga}}$ are the bulk and shear moduli of the fluid-saturated rock.

Equations 1 and 2 give the dry rock moduli in the high-frequency limit, where the inclusions are isolated. In contrast, Eqs. 3 and 4 provide the fluid-saturated rock moduli in the low-frequency limit, where the inclusions are fully connected. Thus, the workflow shown in Fig. 2 can be used to depict the conventional sand-shale reservoirs, especially for the seismic data with low frequencies.

According to Keys and Xu (2002), Eqs. 1 and 2 can be reduced to a simple form by reasonable assumptions about dry rock elastic properties:

$$K_{{\rm d}} = K_{{\rm m}} \left( {1 - \phi } \right)^{{P^{*} }} ,$$

(5)

$$\mu_{{\rm d}} = \mu_{{\rm m}} \left( {1 - \phi } \right)^{{Q^{*} }} ,$$

(6)

where $P^{*} = \sum\limits_{i = 1}^{N} {\phi_{i} P^{{\rm mi}} }$, $Q^{*} = \sum\limits_{i = 1}^{N} {\phi_{i} Q^{{\rm mi}} }$. Thus, the workflow shown in Fig. 2 can be implemented with the dry rock approximation instead of repeatedly solving Eqs. 1 and 2 to estimate the dry rock moduli.

2.2 Applicability Analysis of Gassmann’s Fluid Substitution Equations

Compared with conventional sand-shale reservoirs (Fig. 1a), tight sandstone reservoirs have special petrophysical properties, such as poor connectivity and low porosity (Fig. 3a). Thus, the conventional rock physics modeling may show limitations. In this part, the applicability of Gassmann’s fluid substitution equations for the tight sandstone is analyzed. Due to the low permeability of tight sandstone reservoirs, the assumption that the pore space is fully connected may cause bias. Considering the special petrophysical properties, the pore space in the tight sandstone is divided into four types: connected cracks, connected stiff pores, unconnected cracks, and unconnected stiff pores (Fig. 3b). Its equivalent model is defined here as the partially connected model (Fig. 4a), which can be described by weighting the contributions of the unconnected model (Fig. 4b) and the connected model (Fig. 4c). The building of connected and unconnected models follows the workflows in Figs. 5 and 6, respectively. For the connected model, Gassmann’s equations are used for fluid substitution to ensure that the pores are fully connected (Fig. 5). Instead, the unconnected model is built by adding pore fluid directly into the EMM model (Fig. 6), where the inclusions are isolated. Therefore, by analyzing the differences between these two extreme models, we can understand the limitation of Gassmann’s equations in describing tight sandstone reservoirs.

Fig. 3

Diagrammatic sketch of real rock for tight sandstone reservoirs and its equivalent sketch for rock physics modeling

Fig. 4

Partially connected model (a), unconnected model (b), and connected model (c)

Fig. 5

Rock physics modeling workflow of connected model

Fig. 6

Rock physics modeling workflow of unconnected model

The mineral and fluid parameters used are shown in Table 1. Figure 7 shows the variations in P-wave velocity in different cases. If the rock contains only stiff pores and the pore space is filled with gas, the P-wave velocities are almost the same for the unconnected and connected models. In addition, the pore space of stiff pores filled with water causes a difference in velocity, although this difference is not significant. Thus, for the tight sandstone with stiff pores, the use of Gassmann’s equations can cause a bias, but this bias is small, especially for gas-bearing reservoirs. Notice that the velocities obtained using the parameters in Table 1 are higher for the gas-filled rocks compared with water-filled rocks. This is due to the difference in density. The moduli are actually higher for water-filled rocks (highlighted with black arrows in Fig. 7a). For rocks containing only cracks, the difference in velocity is also small if the pore space is filled with gas. However, the difference becomes large if the pore contains water. Thus, the use of Gassmann’s equations is limited for crack-developed reservoirs. Considering that tight sandstones are prone to cracks, other rock physics modeling methods need to be investigated for tight sandstones.

Table 1 The moduli of mineral/fluid for rock physics modeling

Fig. 7

Variation in P-wave velocity (Vp) with porosity for different petrophysical properties. a refers to the rocks containing only stiff pores; b refers to the rocks containing only cracks. M indicates the P-wave modulus

2.3 Weighted Model for Tight Sandstone

One approach to depict the equivalent model in Fig. 3b is to consider the contributions of both the unconnected and connected models shown in Fig. 4, which was adopted by Wang et al. (2020c) for the reservoir prediction of tight sandstones. The velocities of tight sandstones are given by using the weighted model:

$$V_{{\rm P}} = \sqrt {\frac{{K_{{{\text{tight}}}} + \frac{4}{3}\mu_{{{\text{tight}}}} }}{\rho }} ,$$

(7)

$$V_{{\rm S}} = \sqrt {\frac{{\mu_{{{\text{tight}}}} }}{\rho }} ,$$

(8)

where

$$K_{{\rm tight}} = \lambda K_{{\rm Ga}} + \left( {1 - \lambda } \right)K_{{\rm EMM}} ,$$

(9)

$$\mu_{{\rm tight}} = \lambda \mu_{{\rm Ga}} + \left( {1 - \lambda } \right)\mu_{{\rm EMM}} ,$$

(10)

in which λ indicates the content of connected pores, which is set to 1 if pores are connected, otherwise it is set to 0. The intermediate values imply partial connectivity of the pore space. $K_{{\rm Ga}}$ and $\mu_{{\rm Ga}}$ are the bulk and shear moduli using Gassmann’s equations. $K_{{\rm EMM}}$ and $\mu_{{\rm EMM}}$ are the bulk and shear moduli using the workflow shown in Fig. 6.

Figure 8 shows the variation in P-wave velocity versus porosity. Compared to the rocks with connected pores, the P-wave velocities of unconnected model are greater. These two extreme models are limited to describe the velocity variation in rocks with partially connected pores. The weighted model provides a feasible approach, where λ is a relative parameter which may be affected by pore permeability, gas saturation, fluid viscosity, etc. Following Wang et al. (2020c), the parameter λ can be estimated with well-logging data by using McMC MH algorithm. In their study, both well-logging and seismic data are applied to test this approach. From Fig. 8, the velocities between the upper and lower bounds given by the unconnected and connected models can be described effectively using the weighted model by adjusting the parameter λ. Besides, the velocities also change when the crack content varies, where the crack content refers to the volume content of cracks in the pore space. As the crack content increases, the velocities become smaller. Thus, the fine description of the reservoir requires the investigation of crack content.

Fig. 8

Comparison among unconnected, partially connected and connected models. Vp indicates the P-wave velocity; Sw indicates the water saturation; α indicates the crack content in the pore space

2.4 Same Aspect Ratio Model by Inclusion-Based Method

The theory of inclusion-based method is given in Appendix 2, which was discussed by Song et al. (2016a) to describe the porous rock. By using averaging method and Eshelby’s single-inclusion solution, the strain and stress of the porous rocks can be derived.

According to Wang et al. (2022), the theory of inclusion-based method can be simplified if the pore structure is considered to be the same. The equivalent sketch of the tight sandstone is displayed in Fig. 9, where two types of pore connectivity and one pore aspect ratio are considered. Then, Eq. 52 in Appendix 2 reduces to

$${\mathbf{L}}_{0}^{ - 1} = \left[ {\left( {{\mathbf{L}} - {\mathbf{L}}_{0} } \right)^{ - 1} + \left( {{\mathbf{L}}_{0} - {\mathbf{L}}_{{\rm d}} } \right)^{ - 1} } \right]\left( {{\mathbf{I}} - {\mathbf{L}}_{0} {\mathbf{L}}_{ai}^{ - 1} } \right)\phi^{*} ,$$

(11)

and the dry rock modulus ${\mathbf{L}}_{{\rm d}}$ changes to

$${\mathbf{L}}_{{\rm d}} = {\mathbf{L}}_{0} - {\mathbf{L}}_{0} \left\{ {{\mathbf{S}}_{0} + \left[ {\phi^{*} \left( {{\mathbf{I}} - {\mathbf{S}}_{sr} } \right)^{ - 1} } \right]^{ - 1} } \right\}^{ - 1} ,$$

(12)

where ${\mathbf{L}}_{0}$, ${\mathbf{L}}_{ai}$, and ${\mathbf{L}}$ are the effective elastic moduli of solid matrix, total inclusions and fluid-saturated rock, respectively. ${\mathbf{S}}_{sr}$ is the Eshelby tensor for the inclusion with the same aspect ratio. ${\mathbf{S}}_{0}$ is the Eshelby tensor for the spherical inclusion. Here, this model is called the same aspect ratio model (SARM).

Fig. 9

The equivalent sketch of tight sandstones with the same aspect ratio

The first element of ${\mathbf{L}}_{ai}$ given by

$$M_{{\rm f}}^{\rm IEFM} = {\mathbf{L}}_{ai} \left( {1,1} \right),$$

(13)

is defined as the inclusion-based effective fluid modulus (IEFM).

The introduction of IEFM has no assumption about the connectivity between inclusions, which is different from Gassmann’s equations. Here, we refer to the rock physics model using Gassmann’s equations as the Gassmann model. Figure 10 shows that for the conventional sand-shale rock with connected pores, the Gassmann model is feasible and its predicted crack content is hardly affected by the water saturation. However, for the tight sandstone with unconnected pores, when the pores contain a small amount of gas, the predicted results will cause deviation (highlighted by the black dashed line in Fig. 10a). The predicted crack content and water saturation are also incorrect (blue lines in Fig. 10b). In contrast, when the SARM is adopted, the predictions match well with both unconnected-pore and connected-pore rocks (Fig. 11a). Figure 11b shows that the fluid can be highlighted by IEFM and almost independent of the aspect ratio. The prediction method adopted here is the simulated annealing algorithm. Compared with the Gassmann model, the SARM is superior in the reservoir prediction for tight sandstones (Wang et al. 2022).

Fig. 10

Reservoir prediction using Gassmann model for unconnected-pore and connected-pore rocks. Sw indicates the water saturation; M indicates the P-wave modulus

Fig. 11

Reservoir prediction using SARM for unconnected-pore and connected-pore rocks. Sw indicates the water saturation; r* indicates the aspect ratio of the equivalent pore; M indicates the P-wave modulus

2.5 Proposing a Fluid Pressure Ratio Model

To highlight the difference between cracks and stiff pores, a new inclusion-based model considering the fluid pressure ratio is proposed here, which is named as the fluid pressure ratio model (FPRM). The proposed model is an extension of the SARM. According to Wang et al. (2022), the SARM is based on the KT scheme. However, the KT model shows limitations for the non-dilute inclusions. Instead, the MT scheme is adopted in this part.

Instead of using Eq. 51 in Appendix 2, assuming that the average strain is given by

$${\overline{\mathbf{e}}} = \sum\limits_{i = 0}^{N} {\phi_{i} {\mathbf{e}}_{i} } ,$$

(14)

Equation 52 in Appendix 2 changes to

$$\left( {{\mathbf{L}} - {\mathbf{L}}_{{\rm d}} } \right)^{ - 1} \left\{ {{\mathbf{L}}_{{\rm d}} \left[ {\sum\limits_{i = 1}^{N} {\phi_{i} \left( {{\mathbf{I}} - {\mathbf{S}}_{i} } \right)^{ - 1} {\mathbf{S}}_{i} {\mathbf{L}}_{0}^{ - 1} {{\varvec{\upsigma}}}_{i} } } \right] + \phi^{*} {{\varvec{\upsigma}}}_{ai} } \right\} = \phi^{*} \left( {{\mathbf{L}} - {\mathbf{L}}_{0} } \right)^{ - 1} \left( {{\mathbf{I}} - {\mathbf{L}}_{0} {\mathbf{L}}_{ai}^{ - 1} } \right){{\varvec{\upsigma}}}_{ai} ,$$

(15)

where the elastic modulus tensor of dry rock skeleton ${\mathbf{L}}_{{\rm d}}$ becomes

$${\mathbf{L}}_{{\rm d}} = \left( {1 - \phi^{*} } \right){\mathbf{L}}_{0} \left[ {\left( {1 - \phi^{*} } \right){\mathbf{I}} + \sum\limits_{i = 1}^{N} {\phi_{i} \left( {{\mathbf{I}} - {\mathbf{S}}_{i} } \right)^{ - 1} } } \right]^{ - 1} .$$

(16)

Assuming that the pore space contains two types of pores: cracks and stiff pores, ${{\varvec{\upsigma}}}_{i}$ can be divided into ${{\varvec{\upsigma}}}_{{\rm c}}$ and ${{\varvec{\upsigma}}}_{{\rm p}}$ correspondingly. Then, ${{\varvec{\upsigma}}}_{ai}$ changes to

$${{\varvec{\upsigma}}}_{ai} = \frac{{\phi_{{\rm c}} {{\varvec{\upsigma}}}_{{\rm c}} + \phi_{{\rm p}} {{\varvec{\upsigma}}}_{{\rm p}} }}{{\phi^{*} }} = \alpha {{\varvec{\upsigma}}}_{{\rm c}} + \left( {1 - \alpha } \right){{\varvec{\upsigma}}}_{{\rm p}} ,$$

(17)

where $\phi_{{\rm c}}$ and $\phi_{{\rm p}}$ are the contents of cracks and stiff pores in the rock, respectively. α is the crack content in the pore space, defined as $\alpha = {{\phi_{{\rm c}} } \mathord{\left/ {\vphantom {{\phi_{{\rm c}} } {\phi^{*} }}} \right. \kern-\nulldelimiterspace} {\phi^{*} }}$.

Then, Eq. 15 reduces to

$$\begin{aligned} & \left( {{\mathbf{L}} - {\mathbf{L}}_{{\rm d}} } \right)^{ - 1} \left\{ {{\mathbf{L}}_{{\rm d}} \left[ {\alpha \left( {{\mathbf{I}} - {\mathbf{S}}_{{\rm c}} } \right)^{ - 1} {\mathbf{S}}_{{\rm c}} {\mathbf{L}}_{0}^{ - 1} {{\varvec{\upsigma}}}_{{\rm c}} + \left( {1 - \alpha } \right)\left( {{\mathbf{I}} - {\mathbf{S}}_{{\rm p}} } \right)^{ - 1} {\mathbf{S}}_{{\rm p}} {\mathbf{L}}_{0}^{ - 1} {{\varvec{\upsigma}}}_{{\rm p}} } \right] + \alpha {{\varvec{\upsigma}}}_{{\rm c}} + \left( {1 - \alpha } \right){{\varvec{\upsigma}}}_{{\rm p}} } \right\} \\ & \quad = \left( {{\mathbf{L}} - {\mathbf{L}}_{0} } \right)^{ - 1} \left( {{\mathbf{I}} - {\mathbf{L}}_{0} {\mathbf{L}}_{ai}^{ - 1} } \right)\left[ {\alpha {{\varvec{\upsigma}}}_{{\rm c}} + \left( {1 - \alpha } \right){{\varvec{\upsigma}}}_{{\rm p}} } \right]. \\ \end{aligned}$$

(18)

Letting $\gamma = {{{{\varvec{\upsigma}}}_{{\rm p}} } \mathord{\left/ {\vphantom {{{{\varvec{\upsigma}}}_{{\rm p}} } {{{\varvec{\upsigma}}}_{{\rm c}} }}} \right. \kern-\nulldelimiterspace} {{{\varvec{\upsigma}}}_{{\rm c}} }}$ leads to

$$\begin{aligned} & \left( {{\mathbf{L}} - {\mathbf{L}}_{{\rm d}} } \right)^{ - 1} \left\{ {{\mathbf{L}}_{{\rm d}} \left[ {\alpha \left( {{\mathbf{I}} - {\mathbf{S}}_{{\rm c}} } \right)^{ - 1} {\mathbf{S}}_{{\rm c}} {\mathbf{L}}_{0}^{ - 1} + \gamma \left( {1 - \alpha } \right)\left( {{\mathbf{I}} - {\mathbf{S}}_{{\rm p}} } \right)^{ - 1} {\mathbf{S}}_{{\rm p}} {\mathbf{L}}_{0}^{ - 1} } \right] + \left( {\alpha + \gamma \left( {1 - \alpha } \right)} \right){\mathbf{I}}} \right\} \\ & \quad = \left( {{\mathbf{L}} - {\mathbf{L}}_{0} } \right)^{ - 1} \left( {{\mathbf{I}} - {\mathbf{L}}_{0} {\mathbf{L}}_{ai}^{ - 1} } \right)\left[ {\alpha + \gamma \left( {1 - \alpha } \right)} \right], \\ \end{aligned}$$

(19)

and the elastic modulus tensor ${\mathbf{L}}_{{\rm d}}$ changes to

$${\mathbf{L}}_{{\rm d}} = \left( {1 - \phi^{*} } \right){\mathbf{L}}_{0} \left\{ {\left( {1 - \phi^{*} } \right){\mathbf{I}} + \phi^{*} \left[ {\alpha \left( {{\mathbf{I}} - {\mathbf{S}}_{{\rm c}} } \right)^{ - 1} { + }\left( {1 - \alpha } \right)\left( {{\mathbf{I}} - {\mathbf{S}}_{{\rm p}} } \right)^{ - 1} } \right]} \right\},$$

(20)

where

$${\mathbf{L}}_{ai} = \left( {\begin{array}{*{20}c} {M_{{\rm f}} } & {M_{{\rm f}} - 2\mu_{{\rm f}} } & {M_{{\rm f}} - 2\mu_{{\rm f}} } & 0 & 0 & 0 \\ {M_{{\rm f}} - 2\mu_{{\rm f}} } & {M_{{\rm f}} } & {M_{{\rm f}} - 2\mu_{{\rm f}} } & 0 & 0 & 0 \\ {M_{{\rm f}} - 2\mu_{{\rm f}} } & {M_{{\rm f}} - 2\mu_{{\rm f}} } & {M_{{\rm f}} } & 0 & 0 & 0 \\ 0 & 0 & 0 & {\mu_{{\rm f}} } & 0 & 0 \\ 0 & 0 & 0 & 0 & {\mu_{{\rm f}} } & 0 \\ 0 & 0 & 0 & 0 & 0 & {\mu_{{\rm f}} } \\ \end{array} } \right).$$

(21)

Assuming that the dry and fluid-saturated rocks share the same shear modulus, the shear modulus of fluid is set to $\mu_{{\rm f}} \sim 0$. Thus, the P-wave modulus of fluid M_f is the same as the bulk modulus of fluid K_f. The pore fluid is considered here as a mixed fluid. Thus, the fluid modulus can be calculated using Wood’s equation (Mavko et al. 2020):

$$K_{{\rm f}} = \left( {\frac{{S_{{\rm w}} }}{{K_{{\rm w}} }} + \frac{{1 - S_{{\rm w}} }}{{K_{{\rm g}} }}} \right)^{ - 1} ,$$

(22)

where $K_{{\rm w}}$ and $K_{{\rm g}}$ are the bulk moduli of water and gas, respectively. $S_{{\rm w}}$ is the water saturation.

In this model, the crack content $\alpha$, fluid pressure ratio $\gamma$, and water saturation $S_{{\rm w}}$ are three pending parameters. Next, the connection among them is investigated.

According to Papageorgiou and Chapman (2017), the average pressures of cracks and stiff pores are given by

$$\overline{\sigma }_{{\rm c}} = \frac{{\frac{{1 + i\omega \gamma \tau_{{\rm m}} }}{{3\left( {1 + D_{{\rm c}} } \right)}} + \gamma^{^{\prime}} }}{{1 + \left( {1 + i\omega \tau_{{\rm m}} } \right)\gamma }}\sigma^{*} ,$$

(23)

$$\overline{\sigma }_{{\rm p}} = \frac{{\frac{1}{{3\left( {1 + D_{{\rm c}} } \right)}} + \left( {1 + i\omega \tau_{{\rm m}} } \right)\gamma^{^{\prime}} }}{{1 + \left( {1 + i\omega \tau_{{\rm m}} } \right)\gamma }}\sigma^{*} ,$$

(24)

where $\gamma = \frac{{3\phi_{{\rm p}} \beta_{{\rm c}} \left( {1 + D_{{\rm p}} } \right)}}{{4\phi_{{\rm c}} \mu_{s} \left( {1 + D_{{\rm c}} } \right)}}$; $\gamma^{^{\prime}} = \gamma \frac{{1 - \nu_{s} }}{{\left( {1 + \nu_{s} } \right)\left( {1 + D_{{\rm p}} } \right)}}$; $D_{{\rm p}} = \frac{{4\mu_{s} }}{{3K_{{\rm f}} }}$; $D_{{\rm c}} = \frac{{\beta_{{\rm c}} }}{{K_{{\rm f}} }}$; $\beta_{{\rm c}} = \frac{{\pi \mu_{s} r}}{{2\left( {1 - \nu_{s} } \right)}}$, $\tau_{{\rm m}} = \left( {\frac{{\beta_{{\rm c}} l}}{{\alpha \phi^{*} \left( {1 + D_{{\rm c}} } \right)}}M_{{\rm f}} } \right)^{ - 1}$, and $M_{{\rm f}} = \left( {\frac{{\kappa_{{\rm w}} }}{{\eta_{{\rm w}} }} + \frac{{\kappa_{{\rm g}} }}{{\eta_{{\rm g}} }}} \right)\kappa$. r is the aspect ratio of cracks. σ* is the applied stress. $\nu_{s}$ and $\mu_{s}$ are Poisson’s ratio and shear modulus of solid grains. $l$ is the length scale characteristic of the pore network. $\kappa$ is the intrinsic permeability and $\omega$ is the angle frequency. $\kappa_{{\rm w}}$ and $\kappa_{{\rm g}}$ are the relative permeability of water and gas, respectively. $\eta_{{\rm w}}$ and $\eta_{{\rm g}}$ are the fluid viscosity of water and gas, respectively.

Combining Eqs. 23 and 24, the fluid pressure ratio can be expressed as follows:

$$\frac{{\overline{\sigma }_{{\rm p}} }}{{\overline{\sigma }_{{\rm c}} }} = \frac{{1 + 3\left( {1 + D_{{\rm c}} } \right)\left( {1 + i\omega \tau_{{\rm m}} } \right)\gamma^{^{\prime}} }}{{1 + i\omega \gamma \tau_{{\rm m}} + 3\left( {1 + D_{{\rm c}} } \right)\gamma^{^{\prime}} }}.$$

(25)

Thus, the parameter $\gamma$ is the function of $S_{{\rm w}}$ and $\alpha$:

$$\gamma = f\left( {S_{{\rm w}} ,\alpha } \right) = \hbox{real}\left( {\frac{{\overline{\sigma }_{{\rm p}} }}{{\overline{\sigma }_{{\rm c}} }}} \right).$$

(26)

If the pore space is filled with only gas and water, the gas saturation $S_{{\rm g}}$ satisfies

$$S_{{\rm g}} = 1 - S_{{\rm w}} .$$

(27)

From the above derivation, the three pending parameters become two, which helps to enhance the stability when solving the inverse problem. In the proposed model, the MT scheme is adopted and the fluid pressure ratio is considered, which helps to characterize the difference between cracks and stiff pores and finally improve the fluid identification. Next, both well-logging and seismic data are used to test the proposed method.

3 Applications of Well-Logging and Seismic Data

In this section, we first test the new model with well-logging data. The well-logging data contain a variety of information, including petrophysical properties, such as porosity and permeability, and elastic parameters, such as P- and S-wave velocities. The diversity of the well-logging data helps to test the validity of the proposed method. The field data are from a tight gas sandstone reservoir in north China. Figure 12 shows the measured curves of a typical well. In the study area, the porosity is very small, mostly less than 0.1. The permeability is also small, mostly concentrated within 1 mD (Fig. 13). The study area has the characteristics of low porosity and low permeability (Wang et al. 2020a, 2022). The lithology consists mainly of sand and shale. The moduli of them are estimated from sample points without pores.

Fig. 12

Well logs of a typical well. These logs are P-wave velocity (Vp), S-wave velocity (Vs), density (Den), porosity (Por), gas saturation (Sg) and lithology (Lith) from left to right

Fig. 13

Crossplot of porosity (Por) and permeability (Perm)

Figure 14 shows the comparison between the P-wave and shear moduli predicted using the FPRM and measured data. First, the crack content is adjusted by minimizing the misfit between the predicted P-wave modulus using the FPRM and the measured P-wave modulus. The used parameters are given in Table 2. Other parameters, including shale content, porosity, and water saturation, are obtained from the well-logging data. The prediction method is the McMC MH algorithm. The predicted crack content is shown in the first column of Fig. 14. Then, the corrected crack content is used to obtain the shear modulus based on the FPRM. It can be seen that both the predicted P-wave and shear moduli agree well with the measured data, which shows the viability of the proposed FPRM model. Due to the complexity of pore connectivity and the diversity of fluid in tight sandstones, Gassmann model shows limitations. Thus, the proposed model helps to understand the wave propagation and provides an approach to investigate the pore fluid.

Fig. 14

Comparison between predicted moduli and well-logging data. M and μ indicate the P-wave and shear moduli, respectively; α indicates the predicted crack content in the pore space

Table 2 Parameters for reservoir prediction by this work

Next, the FPRM is used for seismic prediction. Figure 15 shows the near-, mid-, and far-angle stack profiles with an interval of 5 traces, which are the input data for the AVO inversion. Wang et al. (2020a) pointed out that weak reflection events appear with a high probability for tight sandstone reservoirs. For tight sandstones, the reflectivity method is preferred as the forward-modeling engine to estimate the elastic parameters, such as P- and S-wave velocities and density (Liu et al. 2016). The forward theory of reflectivity method is given in Appendix 3. The AVO inversion is implemented in the framework of Bayesian inversion (Zong et al. 2012; Liu et al. 2016; Wang et al. 2020a). Figure 16 shows that the predicted elastic parameters match well with the well-logging data. They are the data base for the fluid discrimination.

Fig. 15

Near-, mid-, and far-angle stack profiles with an interval of 5 traces. a near-angle stack profile; b mid-angle stack profile; c far-angle stack profile. Well X is located at trace 91

Fig. 16

Predicted elastic parameters by using the reflectivity method. a P-wave velocity (Vp); b S-wave velocity (Vs); c density (Den). The displayed well-logging curves are smoothed

Then, 12 wells are used as the training data to obtain the shale content and porosity. In this step, we first derive the posterior probability of shale content by the Bayes discriminant method and then calculate the posterior probability of porosity. The detailed algorithm is shown in Appendix 4. By maximizing the posterior probability, the class that a set of elastic parameters belongs to can be evaluated and the shale content and porosity are finally determined. The predicted shale content and its posterior probability are shown in Fig. 17, and the predicted porosity and its posterior probability are shown in Fig. 18. It can be seen that in addition to obtaining inversion results, the Bayes discriminant method can be used for uncertainty evaluation. The well-logging curves are inserted in the profiles. The predicted results and well-logging data are in good agreement.

Fig. 17

Predicted shale content (SH) and its posterior probability (PP). a Maximum posterior probability result; b posterior probability. The displayed well-logging curve is smoothed

Fig. 18

Predicted porosity (Por) and its posterior probability (PP). a Maximum posterior probability result; b posterior probability. The displayed well-logging curve is smoothed

Finally, the gas saturation is estimated in Bayes posterior prediction framework by using FPRM, where the input data include the predicted elastic parameters, shale content and porosity. The Bayes posterior prediction framework is based on Bayes’ law. The posterior probability $P\left( {{\mathbf{m}}|{\mathbf{d}}} \right)$ is an approximate response of the prior probability $P\left( {\mathbf{m}} \right)$ and the conditional probability $P\left( {{\mathbf{d}}|{\mathbf{m}}} \right)$:

$$P\left( {{\mathbf{m}}|{\mathbf{d}}} \right) \propto P\left( {{\mathbf{d}}|{\mathbf{m}}} \right)P\left( {\mathbf{m}} \right),$$

(28)

where the conditional probability is obtained with the constraint of P-wave modulus (M) and shear modulus ($\mu$):

$$P\left( {{\mathbf{d}}|{\mathbf{m}}} \right) = P\left( {{d}_{M} |{\mathbf{m}}} \right)P\left( {{d}_{\mu } |{\mathbf{m}}} \right)$$

(29)

where d_M and d_μ are seismic P-wave and shear moduli, which can be calculated with the predicted P- and S-wave velocities and density shown in Fig. 16.

$P\left( {{d}_{M} |{\mathbf{m}}} \right)$ and $P\left( {{d}_{\mu } |{\mathbf{m}}} \right)$ share the same form:

$$P\left( {\left. {{\text{d}}_{i} } \right|{\mathbf{m}}} \right) = \frac{1}{{\sqrt {2\pi } \sigma_{i} }}\exp \left\{ { - \frac{{\left[ {{\text{d}}_{i} - {\mathbf{G}}_{i} \left( {\mathbf{m}} \right)} \right]^{2} }}{{2\sigma_{i}^{2} }}} \right\},$$

(30)

where $\sigma_{i}$ is the standard deviation of the noise. ${\mathbf{G}}_{i} \left( {\mathbf{m}} \right)$ is the forward modeling engine, here referring to the FPRM. ${\mathbf{m}}$ indicates the crack content ($\alpha$) and gas saturation ($S_{{\rm g}}$).

Then, Eq. 28 changes to

$$P\left( {\alpha , S_{{\rm g}} |{d}_{M} , {d}_{\mu } } \right) \propto P\left( {{d}_{M} |\alpha , S_{{\rm g}} } \right)P\left( {{d}_{\mu } |\alpha , S_{{\rm g}} } \right)P\left( \alpha \right)P\left( {S_{{\rm g}} } \right).$$

(31)

The McMC MH algorithm is used to obtain the maximum posterior probability solution. The Bayes posterior prediction framework is different from the Bayes discriminant method. The former is a model-driven algorithm, while the latter is a statistical rock physics method belonging to data-driven algorithm.

The predicted gas saturation and its posterior probability are shown in Fig. 19. The logging curve of gas saturation is projected onto the seismic profile at borehole-side trace (Fig. 19a). The predicted profile of gas saturation matches well with the logging data, especially for the gas-bearing zones indicated by the red arrows.

Fig. 19

Predicted gas saturation (Sg) and its posterior probability (PP). a Maximum posterior probability result; b posterior probability

In summary, in the seismic reservoir prediction, three steps are needed for this work: First, the elastic parameters are estimated from angle gathers, where the reflectivity method and the framework of Bayesian inversion are adopted; then, the shale content and porosity profiles are predicted from elastic parameters by using the Bayes discriminant method; finally, the gas saturation is obtained by using the FPRM, where the Bayes posterior prediction framework is applied. This prediction process is helpful for reservoir fluid identification based on seismic data. The seismic data application also demonstrates the validity of the proposed model. The direct reservoir prediction for tight sandstone reservoirs based on Bayesian theory will be our further investigation. Building the direct mapping relationship from angle gathers to petrophysical properties and investigating its high-accuracy approximation are crucial.

4 Conclusions

Fluid discrimination is critical to hydrocarbon exploration, which is the final goal of reservoir prediction. We first review the commonly used inclusion-based rock physics models, seismic elastic parameter inversion methods, and reservoir prediction strategies. Considering that rock physics modeling is the key to link petrophysical properties and elastic parameters, the theoretical background of the most popular modeling workflow for conventional sand-shale reservoirs is also reviewed. Due to the special petrophysical properties of tight sandstones, such as low porosity, low permeability and complex pore structure, conventional rock physics models show limitations. Based on the inclusion-based method, two models for tight sandstones are presented. Then, a new model is proposed by considering the fluid pressure ratio between cracks and stiff pores. Using this model, high-precision predictions of gas saturation can be obtained. Both well-logging data and seismic data demonstrate the feasibility of the proposed method. Fluid discrimination is a complex prediction technique, which needs to combine the connection between seismic data and elastic parameters, as well as the nonlinear relationship between elastic parameters and petrophysical properties. In this paper, a detailed prediction process for tight sandstones is given, where the reflectivity method is used for elastic parameters, and the Bayes posterior prediction framework is adopted for petrophysical properties. Other prediction strategies, especially machine learning algorithms and joint data and model-driven algorithms, are receiving increasing attention. They will be critical tools for the research on unconventional reservoirs such as tight sandstones.

Appendices

Appendix 1: The Coefficients $P^{{\rm mi}}$ and $Q^{{\rm mi}}$ for Randomly Oriented Ellipsoids

The coefficients $P^{{\rm mi}}$ and $Q^{{\rm mi}}$ are two important parameters, which associate the rock modulus with the pore aspect ratio $r$. Following Berryman (1980), these geometric factors are given by

$$P^{{\rm mi}} = \frac{{F_{1} }}{{F_{2} }},$$

(32)

$$Q^{{\rm mi}} = \frac{1}{5}\left( \frac{2}{{F_{3}}} + \frac{1}{{F_{4}}}+\frac{{F_{4} F_{5} + F_{6} F_{7} + F_{8} F_{9} }}{{F_{2} F_{4} }} \right),$$

(33)

where

$$F_{1} = 1 + A\left[ {\frac{3}{2}\left( {f + \vartheta } \right) - R\left( {\frac{3}{2}f + \frac{5}{2}\vartheta - \frac{4}{3}} \right)} \right],$$

(34)

$$\begin{aligned} F_{2} & = 1 + A\left[ {1{ + }\frac{3}{2}\left( {f + \vartheta } \right) - \frac{1}{2}R\left( {3f + 5\vartheta } \right)} \right] + B\left( {3 - 4R} \right) \\ {\kern 1pt} & \quad + \frac{1}{2}A\left( {A + 3B} \right)\left( {3 - 4R} \right)\left[ {f + \vartheta - R\left( {f - \vartheta + 2\vartheta^{2} } \right)} \right], \\ \end{aligned}$$

(35)

$$F_{3} = 1 + A\left[ {1 - \left( {f + \frac{3}{2}\vartheta } \right) + R\left( {f + \vartheta } \right)} \right],$$

(36)

$$F_{4} = 1 + \frac{1}{4}A\left[ {f + 3\vartheta - R\left( {f + \vartheta } \right)} \right],$$

(37)

$$F_{5} = A\left[ { - f + R\left( {f + \vartheta - \frac{4}{3}} \right)} \right] + B\vartheta \left( {3 - 4R} \right),$$

(38)

$$F_{6} = 1 + A\left[ {1 + f - R\left( {f + \vartheta } \right)} \right] + B\left( {1 - \vartheta } \right)\left( {3 - 4R} \right),$$

(39)

$$F_{7} = 2 + \frac{1}{4}A\left[ {3f + 9\vartheta - R\left( {3f + 5\vartheta } \right)} \right] + B\vartheta \left( {3 - 4R} \right),$$

(40)

$$F_{8} = A\left[ {1 - 2R + \frac{1}{2}f\left( {R - 1} \right) + \frac{1}{2}\vartheta \left( {5R - 3} \right)} \right] + B\left( {1 - \vartheta } \right)\left( {3 - 4R} \right),$$

(41)

$$F_{9} = A\left[ {\left( {R - 1} \right)f - R\vartheta } \right] + B\vartheta \left( {3 - 4R} \right),$$

(42)

with $A$, $B$, and $R$ given by

$$A = \frac{{\mu_{i} }}{{\mu_{{\rm m}} }} - 1,$$

(43)

$$B = \frac{1}{3}\left( {\frac{{K_{i} }}{{K_{{\rm m}} }} - \frac{{\mu_{i} }}{{\mu_{{\rm m}} }}} \right),$$

(44)

$$R = \frac{{1 - 2\nu_{{\rm m}} }}{{2\left( {1 - \nu_{{\rm m}} } \right)}},$$

(45)

and

$$\vartheta = \frac{r}{{\left( {1{ - }r^{2} } \right)^{\frac{3}{2}} }}\left[ {\cos^{ - 1} r - r\left( {1 - r^{2} } \right)^{\frac{1}{2}} } \right],$$

(46)

$$f = \frac{{r^{2} }}{{1 - r^{2} }}\left( {3\vartheta - 2} \right).$$

(47)

Appendix 2: The Theory of Inclusion-Based Method

According to Song et al. (2016b), the relation between the average stress over all inclusions ${{\varvec{\upsigma}}}_{ai}$ and the average stress of the entire fluid-saturated porous rock ${\overline{\mathbf{\sigma }}}$ is given by (Hill 1965)

$$\left( {{\mathbf{I}} - {\mathbf{L}}_{0} {\mathbf{L}}^{ - 1} } \right){\overline{\mathbf{\sigma }}} = \phi^{*} \left( {{\mathbf{I}} - {\mathbf{L}}_{0} {\mathbf{L}}_{ai}^{ - 1} } \right){{\varvec{\upsigma}}}_{ai} ,$$

(48)

where ${\overline{\mathbf{\sigma }}} = \sum\nolimits_{i = 0}^{N} {\phi_{i} {{\varvec{\upsigma}}}_{i} }$ and ${{\varvec{\upsigma}}}_{i}$ is the ith component of the rock. $i = 0$ indicates the solid grain; $1 \le i \le N$ refers to the ith inclusion and $N$ is the total number of inclusions. ${{\varvec{\upsigma}}}_{ai} = \frac{{\sum\nolimits_{i = 1}^{N} {\phi_{i} {{\varvec{\upsigma}}}_{i} } }}{{\phi^{*} }}$, ${{\varvec{\upsigma}}}_{ai} = {\mathbf{L}}_{ai} {\mathbf{e}}_{ai}$, ${\mathbf{e}}_{ai} = \frac{{\sum\nolimits_{i = 1}^{N} {\phi_{i} {\mathbf{e}}_{i} } }}{{\phi^{*} }}$, ${\overline{\mathbf{\sigma }}} = {\mathbf{L\overline{e}}}$, and ${\overline{\mathbf{e}}} = \sum\nolimits_{i = 0}^{N} {\phi_{i} {\mathbf{e}}_{i} }$.

For the ith inclusion, the inclusion deformation ${\mathbf{e}}_{i}$ is related to the strain of the reference medium ${\mathbf{e}}_{0}$:

$${\mathbf{e}}_{i} = \left( {{\mathbf{I}} - {\mathbf{S}}_{i} } \right)^{ - 1} {\mathbf{e}}_{0} - \left( {{\mathbf{I}} - {\mathbf{S}}_{i} } \right)^{ - 1} {\mathbf{S}}_{i} {\mathbf{L}}_{0}^{ - 1} {{\varvec{\upsigma}}}_{i} ,$$

(49)

where ${\mathbf{S}}_{i}$ is the Eshelby tensor for the ith inclusion.

Besides, the connection between microstructure and effective elastic properties is given by

$$\left( {{\mathbf{I}} - {\mathbf{L}}_{0} {\mathbf{L}}^{ - 1} } \right){\mathbf{\overline{\sigma } = }}\sum\limits_{i = 1}^{N} {\phi_{i} \left[ {{{\varvec{\upsigma}}}_{i} + {\mathbf{L}}_{0} \left( {{\mathbf{I}} - {\mathbf{S}}_{i} } \right)^{ - 1} {\mathbf{S}}_{i} {\mathbf{L}}_{0}^{ - 1} {{\varvec{\upsigma}}}_{i} - {\mathbf{L}}_{0} \left( {{\mathbf{I}} - {\mathbf{S}}_{i} } \right)^{ - 1} {\mathbf{e}}_{0} } \right]} .$$

(50)

Assuming that the average strain follows

$${\overline{\mathbf{e}}} = {\mathbf{T}}_{0} {\mathbf{e}}_{0} ,$$

(51)

where ${\mathbf{T}}_{0} = \left[ {{\mathbf{I}} + {\mathbf{S}}_{0} {\mathbf{L}}_{0}^{ - 1} \left( {{\mathbf{L}} - {\mathbf{L}}_{0} } \right)} \right]^{ - 1}$ and combing Eqs. 48, 49 and 50 leads to

$$\begin{aligned} \sum\limits_{i = 1}^{N} {\phi_{i} \left( {{\mathbf{I}} - {\mathbf{S}}_{i} } \right)^{ - 1} {\mathbf{L}}_{0}^{ - 1} {{\varvec{\upsigma}}}_{i} } & = \left[ {\sum\limits_{i = 1}^{N} {\phi_{i} \left( {{\mathbf{I}} - {\mathbf{S}}_{i} } \right)^{ - 1} } } \right]\left[ {\left( {{\mathbf{L}} - {\mathbf{L}}_{0} } \right)^{ - 1} + \left( {{\mathbf{L}}_{0} - {\mathbf{L}}_{{\rm d}} } \right)^{ - 1} } \right] \\ & \quad \times \left( {{\mathbf{I}} - {\mathbf{L}}_{0} {\mathbf{L}}_{ai}^{ - 1} } \right)\sum\limits_{i = 1}^{N} {\phi_{i} {{\varvec{\upsigma}}}_{i} } , \\ \end{aligned}$$

(52)

where ${\mathbf{L}}_{{\rm d}}$ is the elastic modulus of the dry rock frame:

$${\mathbf{L}}_{{\rm d}} = {\mathbf{L}}_{0} - {\mathbf{L}}_{0} \left\{ {{\mathbf{S}}_{0} + \left[ {\sum\limits_{i = 1}^{N} {\phi_{i} \left( {{\mathbf{I}} - {\mathbf{S}}_{i} } \right)^{ - 1} } } \right]^{ - 1} } \right\}^{ - 1} .$$

(53)

Appendix 3: Forward Theory of the Reflectivity Method

The reflectivity method is an effective way to investigate the wave propagation by considering the transmission loss and internal multiples. The reflectivity coefficient $R_{{{\text{PP}}}} (p,\omega )$ is the total response in the frequency-slowness domain for a given stack of $N$ layers (Wang et al. 2020a):

$$R_{{{\text{PP}}}} (p,\omega ) = \frac{{{\varvec{\nu}}_{0,4} }}{{{\varvec{\nu}}_{0,1} }},$$

(54)

where ${\varvec{\nu}}_{0,1}$ and ${\varvec{\nu}}_{0,4}$ are the first and fourth elements of the top-layer response ${\varvec{\nu}}_{0}$, respectively. The vector ${\varvec{\nu}}_{0}$ is given by

$${\varvec{\nu}}_{0} = {\mathbf{Q}}_{0} {\mathbf{Q}}_{1} \cdots {\mathbf{Q}}_{i} \cdots {\mathbf{Q}}_{N - 1} {\varvec{\nu}}_{N} ,$$

(55)

where ${\varvec{\nu}}_{N} = \left[ {1\begin{array}{*{20}l} {} \hfill \\ \end{array} 0\begin{array}{*{20}l} {} \hfill \\ \end{array} 0\begin{array}{*{20}l} {} \hfill \\ \end{array} 0\begin{array}{*{20}l} {} \hfill \\ \end{array} 0\begin{array}{*{20}l} {} \hfill \\ \end{array} 0} \right]^{{\rm T}}$ is the starting vector at the bottom and ${\mathbf{Q}}_{i}$ is the propagator matrix of each layer

$${\mathbf{Q}}_{i} = {\mathbf{T}}_{i}^{ + } {\mathbf{E}}_{i} {\mathbf{T}}_{i}^{ - } ,$$

(56)

in which ${\mathbf{E}}_{i}$ is the phase shift matrix of reflection

$${\mathbf{E}}_{i} = diag\left[ {\begin{array}{*{20}c} {e^{{ - i\omega \Delta h\left( {q_{{\rm p}} + q_{s} } \right)}} } & 1 & {e^{{ - i\omega \Delta h\left( {q_{{\rm p}} - q_{s} } \right)}} } & {e^{{i\omega \Delta h\left( {q_{{\rm p}} - q_{s} } \right)}} } & 1 & {e^{{i\omega \Delta h\left( {q_{{\rm p}} + q_{s} } \right)}} } \\ \end{array} } \right],$$

(57)

and ${\mathbf{T}}_{i}^{ + }$ and ${\mathbf{T}}_{i}^{ - }$ are the downward and upward interface crossing matrices. The matrices ${\mathbf{T}}_{i}^{ + }$ and ${\mathbf{T}}_{i}^{ - }$ are expressed as follows:

$${\mathbf{T}}_{i}^{ + } = \left[ {\begin{array}{*{20}c} { - \left( {p^{2} + q_{{\rm p}} q_{s} } \right)/\mu } & { - 2pq_{{\rm p}} /\mu } & { - \left( {p^{2} - q_{{\rm p}} q_{s} } \right)/\mu } \\ {iq_{s} /\beta^{2} } & 0 & { - iq_{s} /\beta^{2} } \\ { - ip\left( {\Gamma + 2q_{{\rm p}} q_{s} } \right)} & { - 4ip^{2} q_{{\rm p}} } & { - ip\left( {\Gamma - 2q_{{\rm p}} q_{s} } \right)} \\ { - ip\left( {\Gamma + 2q_{{\rm p}} q_{s} } \right)} & { - 2i\Gamma q_{{\rm p}} } & { - ip\left( {\Gamma - 2q_{{\rm p}} q_{s} } \right)} \\ { - iq_{{\rm p}} /\beta^{2} } & 0 & { - iq_{{\rm p}} /\beta^{2} } \\ { - \mu \left( {\Gamma^{2} + 4p^{2} q_{{\rm p}} q_{s} } \right)} & { - 4\mu \Gamma pq_{{\rm p}} } & { - \mu \left( {\Gamma^{2} - 4p^{2} q_{{\rm p}} q_{s} } \right)} \\ \end{array} } \right. \cdots \quad \cdots \left. {\begin{array}{*{20}c} {\left( {p^{2} - q_{{\rm p}} q_{s} } \right)/\mu } & { - 2pq_{s} /\mu } & { - \left( {p^{2} + q_{{\rm p}} q_{s} } \right)/\mu } \\ { - iq_{s} /\beta^{2} } & 0 & { - iq_{s} /\beta^{2} } \\ {ip\left( {\Gamma - 2q_{{\rm p}} q_{s} } \right)} & { - 2i\Gamma q_{s} } & { - ip\left( {\Gamma + 2q_{{\rm p}} q_{s} } \right)} \\ {ip\left( {\Gamma - 2q_{{\rm p}} q_{s} } \right)} & { - 4ip^{2} q_{s} } & { - ip\left( {\Gamma + 2q_{{\rm p}} q_{s} } \right)} \\ { - iq_{{\rm p}} /\beta^{2} } & 0 & {iq_{{\rm p}} /\beta^{2} } \\ {\mu \left( {\Gamma^{2} - 4p^{2} q_{{\rm p}} q_{s} } \right)} & { - 4\mu \Gamma pq_{s} } & { - \mu \left( {\Gamma^{2} + 4p^{2} q_{{\rm p}} q_{s} } \right)} \\ \end{array} } \right],$$

(58)

and

$${\mathbf{T}}_{i}^{ - } = \left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {{\text{T}}_{i}^{ + } (6,1)} & {{\text{T}}_{i}^{ + } (5,1)} & {{\text{T}}_{i}^{ + } (3,1)} \\ { - {\text{T}}_{i}^{ + } (6,5)} & 0 & { - {\text{T}}_{i}^{ + } (4,5)} \\ { - {\text{T}}_{i}^{ + } (6,3)} & { - {\text{T}}_{i}^{ + } (5,1)} & { - {\text{T}}_{i}^{ + } (3,3)} \\ {{\text{T}}_{i}^{ + } (6,3)} & { - {\text{T}}_{i}^{ + } (5,1)} & {{\text{T}}_{i}^{ + } (3,3)} \\ { - {\text{T}}_{i}^{ + } (6,2)} & 0 & { - {\text{T}}_{i}^{ + } (4,2)} \\ {{\text{T}}_{i}^{ + } (6,1)} & { - {\text{T}}_{i}^{ + } (5,1)} & {{\text{T}}_{i}^{ + } (3,1)} \\ \end{array} } & {\begin{array}{*{20}c} {{\text{T}}_{i}^{ + } (3,1)} & {{\text{T}}_{i}^{ + } (2,1)} & {{\text{T}}_{i}^{ + } (1,1)} \\ { - {\text{T}}_{i}^{ + } (3,5)} & 0 & { - {\text{T}}_{i}^{ + } (1,5)} \\ { - {\text{T}}_{i}^{ + } (3,3)} & {{\text{T}}_{i}^{ + } (1,2)} & { - {\text{T}}_{i}^{ + } (1,3)} \\ {{\text{T}}_{i}^{ + } (3,3)} & {{\text{T}}_{i}^{ + } (2,1)} & {{\text{T}}_{i}^{ + } (1,3)} \\ { - {\text{T}}_{i}^{ + } (3,2)} & 0 & { - {\text{T}}_{i}^{ + } (1,2)} \\ {{\text{T}}_{i}^{ + } (3,1)} & { - {\text{T}}_{i}^{ + } (2,1)} & {{\text{T}}_{i}^{ + } (1,1)} \\ \end{array} } \\ \end{array} } \right],$$

(59)

where $q_{{\rm p}} = \sqrt {1/\psi^{2} - p^{2} }$, $q_{s} = \sqrt {1/\beta^{2} - p^{2} }$, $\Gamma = 2p^{2} - 1/\beta^{2}$, and $\mu = \rho \beta^{2}$. $\Delta h$, $\psi$, $\beta$, $p$, and $\rho$ are the thickness, P- and S-wave velocities, ray parameter, and density of the ith layer, respectively.

Then, $R_{{{\text{PP}}}} (p,\omega )$ is changed to $R_{{{\text{PP}}}} (\theta ,\omega )$ by using Snell’s law $p = {{\sin (\theta )/ \psi}}$ and $R_{{{\text{PP}}}} (\theta ,\omega )$ is finally changed to $R_{{{\text{PP}}}} (\theta ,t)$ by inverse Fourier transform. $R_{{{\text{PP}}}} (\theta ,t)$ belongs to the time-angle domain and can be used as the forward modeling engine in the seismic AVO inversion.

Appendix 4: Bayes Discriminant Method

According to Bayes’ law, the posterior probability that for a given ${\mathbf{X = x}}$, ${\mathbf{Y}}$ belongs to the class $c_{k}$ is given by

$$P\left( {{\mathbf{Y}} = c_{k} \left| {{\mathbf{X}} = {\mathbf{x}}} \right.} \right)\frac{{P\left( {{\mathbf{X}} = {\mathbf{x}}\left| {{\mathbf{Y}} = c_{k} } \right.} \right)P\left( {{\mathbf{Y}} = c_{k} } \right)}}{{\sum\nolimits_{k} {P\left( {{\mathbf{X}} = {\mathbf{x}}\left| {{\mathbf{Y}} = c_{k} } \right.} \right)P\left( {{\mathbf{Y}} = c_{k} } \right)} }},$$

(60)

where $P\left( {{\mathbf{Y}} = c_{k} } \right)$ and $P\left( {{\mathbf{X}} = {\mathbf{x}}\left| {{\mathbf{Y}} = c_{k} } \right.} \right)$ are prior and conditional probabilities, respectively, which are obtained by training well-logging data. By maximizing Eq. 60, the class that ${\mathbf{X = x}}$ belongs to can be estimated.

Considering the independent identically distributed assumption between the different parameters of ${\mathbf{X}}$, the conditional probability reduces to

$$P({\mathbf{X}} = {\mathbf{x}}|{\mathbf{Y}} = c_{k} ) = \mathop \prod \limits_{{i = 1}}^{n} P({\mathbf{X}} = x^{i} |{\mathbf{Y}} = c_{k} ),$$

(61)

where $x^{i}$ is the ith feature of ${\mathbf{X}}$.

In this work, the posterior probability of shale content is expressed as follows:

$$P({\mathbf{Y}} = sh|{\mathbf{X}} = V_{{\rm P}} ,V_{{\rm S}} ,\rho ),$$

(62)

where ${\rm sh}$, $V_{{\rm P}}$, $V_{{\rm S}}$, and $\rho$ are shale content, P- and S-wave velocities, and density, respectively.

The posterior probability of porosity (por) is given in an integral form:

$$\begin{aligned} P({\mathbf{Y}} = {\text{por}}|{\mathbf{X}} = V_{{\rm P}} ,V_{{\rm S}} ,\rho ) & = \int {P({\mathbf{Y}} = {{\rm sh}}|{\mathbf{X}} = V_{{\rm P}} ,V_{{\rm S}} ,\rho )} \\ \quad \times P({\mathbf{Y}} = {\text{por}}|{\mathbf{X}} = V_{{\rm P}} ,V_{{\rm S}} ,\rho ,{\rm sh}){d}\left( {{\rm sh}} \right). \\ \end{aligned}$$

(63)