Inversion Correction Method for Pore and Fracture Rock Physics Model and Azimuthal Seismic Reflection Characteristic Analysis

1. Introduction

Rock physics models act as an essential bridge connecting seismic elastic parameters with reservoir physical parameters. Generally, conventional rock physics modeling can be implemented via three standard procedures. Firstly, different minerals are homogenized to determine the elastic moduli of equivalent solid matrices. Secondly, pores and fractures are introduced into the solid matrices to calculate the elastic moduli of dry rock. Thirdly, various pore fluids are mixed and injected into dry rock to acquire the elastic moduli of fluid-saturated rocks. Among the whole modeling process, equivalent medium theory lays a theoretical foundation and plays a decisive role.

Equivalent medium theory is the fundamental theoretical support for rock physics modeling. According to whether the pore geometry effect is taken into consideration, existing equivalent medium models can be classified into two categories: models ignoring pore shape effects (i.e., bound models) and models considering pore shape effects (i.e., self-consistent models). Bound models were first proposed in the 1920s. Among them, the upper bound of elastic modulus for mixed mineral compositions pro-posed by Voigt [1] and the lower bound proposed by Reuss [2] are two most classical theoretical frameworks. Apart from verifying the rationality of newly established rock physics models, the average value of Voigt and Reuss bounds can be adopted to calculate the average elastic modulus of minerals and further estimate the elastic properties of equivalent solid matrices. This widely used averaging method is defined as the Voigt-Reuss-Hill (VRH) average theory [3].

When calculating the elastic moduli of dry rock, it is generally assumed that pores within the solid matrix dominate the isotropic elastic properties of media. Three classic and representative isotropic equivalent medium models that account for pores and their morphological effects are widely adopted in rock physics research. Specifically, the Kuster-Toksöz (K-T) model is derived based on the first-order wave scattering theory [4]. The self-consistent approximation (SCA) model requires no predefined background medium and treats all mineral components equally without priority [5]. Different from the SCA model, the differential effective medium (DEM) model distinguishes primary and secondary mineral components according to the sequence of mineral addition during modeling [6].

In contrast, oriented fractures will induce anisotropic elastic characteristics of rock media during the calculation of elastic moduli for dry rock frameworks. For anisotropic media with ellipsoidal inclusions, Hudson established a classic rock physics model applicable to thin penny-shaped fractures and first proposed the definition of fracture density [7]. Nevertheless, the original Hudson model suffers from obvious limitations, particularly its poor performance for fractures with small aspect ratios. Subsequent theoretical derivations have removed unreasonable assumptions of the original model and greatly extended its applicable scope [8]. On the other hand, when fractures are simplified as vertically oriented infinite planar cracks, the fractured rock medium can be equivalent to HTI (Horizontal Transverse Isotropy) medium, and the Schoenberg theory suitable for slender fractures has been developed correspondingly [9]. With continuous improvement, the Schoenberg theory can further characterize the variation of elastic parameters caused by changing fracture aspect ratios and has achieved favorable application performance in field reservoir evaluation. By integrating linear slip theory with matrix pore-based equivalent medium theory, a novel pore-fracture medium model with a single set of vertically aligned fractures is constructed, which can quantitatively analyze the respective contributions of pores and fractures to rock elastic properties. Currently, rock physics theories coupling pores and fractures have been well developed, and such integrated models have been widely applied to cope with complex actual geological conditions [10,11,12,13,14,15].

On the basis of previous research, we develop a novel rock physics model considering both pore shape and the azimuth of oriented fractures. To improve the accuracy of the established model, model correction is performed via inversion for elastic moduli of minerals and oriented fracture parameters, which realizes a fine matching between the theoretical model and logging data. Thereafter, an inversion objective function is constructed based on the corrected rock physics model to quantitatively predict the azimuth of oriented fractures. Ultimately, combined with the Rüger approximate equation for seismic re-flection coefficients [16], the proposed rock physics model is adopted to theoretically clarify the response mechanism of reservoir physical parameters to seismic reflection characteristics.

2. Materials and Methods

2.1. Rock Physics Model with Pores and Oriented Fractures

Rock physics modeling is to establish a functional relationship between seismic elastic parameters and reservoir physical parameters. Rock physics modeling mainly consists of three steps: Mix minerals to calculate the equivalent elastic modulus of the matrix; Introduce pores and fractures to compute the elastic modulus of the dry rock; Incorporate fluids to obtain the elastic parameters of fluid-saturated rock. The detailed modeling workflow is as follows:

• The Voigt-Reuss-Hill average method [3] is adopted to mix minerals and calculate the equivalent elastic moduli of the rock matrix. It is assumed that the minerals of sandstone consist of sandy minerals and clay minerals. The expressions for equivalent bulk modulus, shear modulus and density of rock matrix are given as follows:

  $$K_{mat}=\left[\left(V_{sand}{K}_{sand}+V_{clay}{K}_{clay}\right)+{\left(V_{sand}/K_{sand}+V_{clay}/K_{clay}\right)}^{-1}\right]/2,$$

  (1)

  $${\mu }_{mat}=\left[\left(V_{sand}{\mu }_{sand}+V_{clay}{\mu }_{clay}\right)+{\left(V_{sand}/{\mu }_{sand}+V_{clay}/{\mu }_{clay}\right)}^{-1}\right]/2,$$

  (2)

  $${\rho }_{mat}=V_{sand}{\rho }_{sand}+V_{clay}{\rho }_{clay},$$

  (3)
  In Equations (1–3): $V_{clay}$denotes clay content, $V_{sand}$denotes sand fraction. We assume $V_{sand}=1-V_{clay}$, meaning the rock matrix only contains two mineral components: sandy minerals and clay minerals. $K_{sand}$and ${\mu }_{sand}$ are the bulk modulus and shear modulus of sandy minerals, respectively, and ${\rho }_{sand}$ is the density of sandy minerals; $K_{clay}$ and ${\mu }_{clay}$ are the bulk modulus and shear modulus of clay minerals, respectively, and ${\rho }_{clay}$ is the density of clay minerals. In addition, $K_V$ in Equation. (1) and ${\mu }_V$ in Equation (2) represent the Voigt upper bounds of the bulk modulus and shear modulus of the rock matrix, while $K_R$ and ${\mu }_R$ stand for the Reuss lower bounds of the rock matrix bulk modulus and shear modulus, respectively;
• The equivalent bulk modulus and density of fluids are calculated using the Wood equation [17]:

  $$K_{fl}=S_g/K_g+S_w/K_w,$$

  (4)

  $${\rho }_{fl}=S_g{\rho }_g+S_w{\rho }_w,$$

  (5)

  where $S_w$ denotes water saturation, $S_g$ denotes gas saturation. We assume $S_g=1-S_w$, meaning the fluid is a mixture of water and natural gas; $K_w$ and $K_g$are the bulk moduli of water and gas, respectively; ${\rho }_w$ and ${\rho }_g$ are the densities of water and gas, respectively; $K_{fl}$ is the equivalent bulk modulus of mixed pore fluids, and ${\rho }_{fl}$ is the equivalent density of mixed pore fluids;
• The SCA model [6] is adopted to incorporate pores with variable pore aspect ratios into the equivalent rock matrix:

  $$K_{SC}=K_{mat}+{\phi }_p\left(K_{fl}-K_{mat}\right)P^{*}\left({\alpha }_p\right),$$

  (6)

  $${\mu }_{SC}={\mu }_{mat}-{\phi }_p\cdot {\mu }_{mat}\cdot Q^{*}\left({\alpha }_p\right),$$

  (7)

  where $K_{SC}$ and ${\mu }_{SC}$ denote the bulk modulus and shear modulus of rock after embedding fluid-saturated pores, respectively; ${\phi }_p$ represents matrix porosity, and ${\alpha }_p$ denotes the aspect ratio of pores; $P^{*}\left({\alpha }_p\right)$ and $Q^{*}\left({\alpha }_p\right)$ are geometric coefficients accounting for variations in pore aspect ratios;
• The Schoenberg model [9] is applied to introduce fractures with variable strike azimuths into the rock containing pores. First, the Lamé constant ${\lambda }_{SC}$ is calculated using the bulk modulus $K_{SC}$ and shear modulus ${\mu }_{SC}$ of the pore-bearing rock:

  $${\lambda }_{SC}=K_{SC}-{\displaystyle \frac{2}{3}}{\mu }_{SC},$$

  (8)

Assuming the oriented fractures are vertical with variable strike azimuths, the stiffness matrix of the rock physical model accounting for pore geometries and strike angles of aligned fractures can be expressed as:

$$C=C_0+C_S,$$

(9)

where $C_0$ and $C_S$ represent the stiffness matrix of the isotropic background rock medium containing only pores and the stiffness perturbation induced by oriented fractures, respectively:

$$C_0=\left[\begin{array}{cccccc}{\lambda }_{SC}+2{\mu }_{SC}& {\lambda }_{SC}& {\lambda }_{SC}& 0& 0& 0\\ {\lambda }_{SC}& {\lambda }_{SC}+2{\mu }_{SC}& {\lambda }_{SC}& 0& 0& 0\\ {\lambda }_{SC}& {\lambda }_{SC}& {\lambda }_{SC}+2{\mu }_{SC}& 0& 0& 0\\ 0& 0& 0& {\mu }_{SC}& 0& 0\\ 0& 0& 0& 0& {\mu }_{SC}& 0\\ 0& 0& 0& 0& 0& {\mu }_{SC}\end{array}\right],$$

(10)

$$C_S=\left[\begin{array}{cccccc}\left({\lambda }_{SC}+2{\mu }_{SC}\right)\left(-\mathsf{\Delta }N\right)& {\lambda }_{SC}\left(-\mathsf{\Delta }N\right)& {\lambda }_{SC}\left(-\mathsf{\Delta }N\right)& 0& 0& 0\\ {\lambda }_{SC}\left(-\mathsf{\Delta }N\right)& \left({\lambda }_{SC}+2{\mu }_{SC}\right)\left(-{\left({\displaystyle \frac{{\lambda }_{SC}}{{\lambda }_{SC}+2{\mu }_{SC}}}\right)}^2\mathsf{\Delta }N\right)& {\lambda }_{SC}\left(-{\displaystyle \frac{{\lambda }_{SC}}{{\lambda }_{SC}+2{\mu }_{SC}}}\mathsf{\Delta }N\right)& 0& 0& 0\\ {\lambda }_{SC}\left(-\mathsf{\Delta }N\right)& {\lambda }_{SC}\left(-{\displaystyle \frac{{\lambda }_{SC}}{{\lambda }_{SC}+2{\mu }_{SC}}}\mathsf{\Delta }N\right)& \left({\lambda }_{SC}+2{\mu }_{SC}\right)\left(-{\left({\displaystyle \frac{{\lambda }_{SC}}{{\lambda }_{SC}+2{\mu }_{SC}}}\right)}^2\mathsf{\Delta }N\right)& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& {\mu }_{SC}\left(-\mathsf{\Delta }T\right)& 0\\ 0& 0& 0& 0& 0& {\mu }_{SC}\left(-\mathsf{\Delta }T\right)\end{array}\right],$$

(11)

The expressions for the oriented fracture parameters $\mathsf{\Delta }N$ and $\mathsf{\Delta }T$ in Equation (11) are given as follows:

$$\mathsf{\Delta }N=4e/\{3g_{SC}(1-g_{SC})[1+{\displaystyle \frac{1}{\pi {\alpha }_f(1-g_{SC})}}{\displaystyle \frac{K_{fl}}{{\mu }_{SC}}}]\},$$

(12)

$$\mathsf{\Delta }T=16e/\left[3(3-2g_{SC})\right],$$

(13)

where ${\alpha }_f$ denotes the aspect ratio of oriented fractures, whose value is generally much smaller than that of pores. $g_{SC}$ is the squared ratio of S-wave velocity to P-wave velocity of the rock with pores, $e$ represents the oriented fracture density, and the expressions for $g_{SC}$ and $e$ are written as:

$$g_{SC}={\displaystyle \frac{{\mu }_{SC}}{{\lambda }_{SC}+2{\mu }_{SC}}},$$

(14)

$$e={\displaystyle \frac{3{\phi }_f}{4\pi {\alpha }_f}},$$

(15)

where ${\phi }_f$ denotes the porosity of oriented fractures, whose magnitude is generally far smaller than pore porosity. Assuming the strike azimuth of oriented fractures is ${\varphi }_0$, the stiffness matrix $C_{sat}$ of fluid-saturated rock is written as:

$$C_{sat}=M_{{\varphi }_0}C{M}_{{\varphi }_0}^T,$$

(16)

In Equation (16):

$$M_{{\varphi }_0}=\left[\begin{array}{cccccc}{\cos }^2\left({\varphi }_0\right)& {\sin }^2\left({\varphi }_0\right)& 0& 0& 0& -\sin \left(2{\varphi }_0\right)\\ {\sin }^2\left({\varphi }_0\right)& {\cos }^2\left({\varphi }_0\right)& 0& 0& 0& \sin \left(2{\varphi }_0\right)\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& \cos \left({\varphi }_0\right)& \sin \left({\varphi }_0\right)& 0\\ 0& 0& 0& -\sin \left({\varphi }_0\right)& \cos \left({\varphi }_0\right)& 0\\ {\displaystyle \frac{1}{2}}\sin \left(2{\varphi }_0\right)& -{\displaystyle \frac{1}{2}}\sin \left(2{\varphi }_0\right)& 0& 0& 0& \cos \left(2{\varphi }_0\right)\end{array}\right],$$

(17)

Given the stiffness matrix $C_{sat}$ of fluid-saturated rock, the P-wave velocity $V_P$ and S-wave velocity $V_S$, as well as three HTI anisotropy parameters [16] ${\epsilon }^{\left(V\right)}$, ${\gamma }^{\left(V\right)}$ and ${\delta }^{\left(V\right)}$, can be calculated therefrom. In addition, the bulk density $\rho$ is obtained via weighted averaging, are written as:

$$V_P=\sqrt{{\displaystyle \frac{C_{sat}\left(3,3\right)}{\rho }}},$$

(18)

$$V_S=\sqrt{{\displaystyle \frac{C_{sat}\left(4,4\right)}{\rho }}},$$

(19)

$$\rho =\left(1-{\phi }_p-{\phi }_f\right){\rho }_{mat}+\left({\phi }_p+{\phi }_f\right){\rho }_{fl},$$

(20)

$${\epsilon }^{\left(V\right)}={\displaystyle \frac{C_{sat}\left(1,1\right)-C_{sat}\left(3,3\right)}{2C_{sat}\left(3,3\right)}},$$

(21)

$${\gamma }^{\left(V\right)}={\displaystyle \frac{C_{sat}\left(6,6\right)-C_{sat}\left(4,4\right)}{2C_{sat}\left(4,4\right)}},$$

(22)

$${\delta }^{\left(V\right)}={\displaystyle \frac{{\left[C_{sat}\left(1,3\right)+C_{sat}\left(6,6\right)\right]}^2-{\left[C_{sat}\left(3,3\right)-C_{sat}\left(6,6\right)\right]}^2}{2C_{sat}\left(3,3\right)\left[C_{sat}\left(3,3\right)-C_{sat}\left(6,6\right)\right]}},$$

(23)

Thus far, the construction of the rock physical model accounting for pores and oriented fractures has been completed. The mapping relationship between the anisotropic elastic parameters of rock and physical parameters $f_{RPM}^{HTI}$ is written as:

$$\left[V_P,V_S,\rho ,{\epsilon }^{\left(V\right)},{\gamma }^{\left(V\right)},{\delta }^{\left(V\right)}\right]=f_{RPM}^{HTI}\left({\phi }_p,{\alpha }_p,{\phi }_f,{\alpha }_f,{\varphi }_0,V_{clay},S_w\right),$$

(24)

2.2. Stochastic Inversion Algorithm with Variable Sampling Space

Given logging data containing anisotropic elastic parameters (P-wave and S-wave velocities, density, three anisotropy parameters) and three physical parameters (porosity, clay content and water saturation), we aim to invert unknown parameters including pore aspect ratio and fracture strike azimuth. To address this complicated nonlinear multi-parameter inversion issue, a variable sampling space stochastic inversion method [18] is introduced. This algorithm is applied in two stages to correct the rock physics model for well log datasets.

Figure 1 shows the flowchart of the stochastic inversion algorithm with variable sampling space, and its detailed procedures are described as follows:

• Firstly, assign the prior probability density function, then generate N random sample points following this distribution;
• During every iteration, we retain the top $aqa\cdot N$ samples with the smallest objective function values. Assuming these samples obey Gaussian statistics, we compute their new mean $\mu$ and variance ${\sigma }^2$;
• Resample N sample points based on the updated Gaussian distribution$X\sim N\left(\mu ,apa\cdot {\sigma }^2\right)$, and then launch the next iteration. This procedure is repeated until the algorithm converges.

Figure 1 assumes lower objective function values represent smaller calculation errors. If higher objective function values correspond to smaller errors, samples are sorted descendingly. Input parameters of the inversion algorithm: sample count N, maximum iterations, shrinkage factor $aqa$ ($aqa\cdot N$ must be integer), amplification factor $apa$ . Factor $aqa\in \left(0,1\right)$ extracts well-fitted $aqa\cdot N$ samples for statistics. Factor $apa$ enlarges variance to expand sample space and enrich candidate cases in resampling.

3. Results

The logging dataset contains P-wave velocity, S-wave velocity, bulk density, three anisotropy parameters (${\epsilon }^{\left(V\right)}$, ${\delta }^{\left(V\right)}$and ${\gamma }^{\left(V\right)}$), and three physical parameters (porosity, clay content, and water saturation). Taking Well A as an example, Figure 2 illustrates its full logging profiles, encompassing nine elastic and physical parameters in total.

3.1. Inversion-Based Correction of Rock Physics Model

From the principle of rock physics modeling, the rock physics model behaves isotropically when the effect of oriented fractures is not incorporated. In reservoir conditions, the porosity of oriented fractures is generally much smaller than the total porosity of the rock. Under the assumption that fracture porosity is far lower than total porosity, P-wave velocity, S-wave velocity, density, porosity, clay content and water saturation are combined with the isotropic component of the rock-physical model (i.e., the isotropic rock physics model that only accounts for pores) to invert the elastic moduli of minerals for model correction.

In this step, the three anisotropy parameters are temporarily excluded, and the other six logging parameters are treated as known inputs. A stochastic inversion algorithm with variable sampling space is utilized to invert mineral elastic moduli and correct the isotropic component of the rock physics model. Taking logging-derived porosity, clay content and water saturation as inputs, Figure 3 compares the elastic parameters computed from the corrected isotropic rock physics model with the elastic parameters acquired from well logging measurements.

In Figure 3, the red curves denote logging data, and the blue curves stand for P-wave velocity, S-wave velocity and density calculated from the corrected isotropic component of the rock physics model. The correction is realized by inverting mineral elastic moduli using a stochastic inversion algorithm with variable sampling space, where logging porosity, clay content and water saturation serve as input variables. Density is calculated via Equation (20). The equivalent densities of sand and clay minerals are predefined ${\rho }_{sand}=2.56g/c{m}^3$ and ${\rho }_{clay}=2.58g/c{m}^3$. Pore fluids are treated as a binary mixture of natural gas and water, where the shear modulus of fluid is zero. The bulk modulus and density of natural gas $K_g=1.0\mathrm{GPa}$ and ${\rho }_g=0.8g/c{m}^3$, as well as the bulk modulus and density of water $K_w=2.38\mathrm{GPa}$ and ${\rho }_w=1.0g/c{m}^3$, are specified. Using logging data to invert the elastic moduli of minerals, the objective function for correcting the isotropic component of the rock physics model is established as follows:

$$E={\left[V_{P,well}-V_{P,iso}\left(K_{sand},K_{clay},{\mu }_{sand},{\mu }_{clay}\right)\right]}^2+{\left[V_{S,well}-V_{S,iso}\left(K_{sand},K_{clay},{\mu }_{sand},{\mu }_{clay}\right)\right]}^2,$$

(25)

where $V_{P,well}$ and $V_{S,well}$ denote the measured P- and S- wave velocities from well logging, respectively; $V_{P,iso}$ and $V_{S,iso}$ are the P- and S-wave velocities calculated by the isotropic component of the rock physics model, which are written as:

$$V_{P,iso}=\sqrt{{\displaystyle \frac{K_{SC}+\frac{4}{3}{\mu }_{SC}}{\rho }}},$$

(26)

$$V_{S,iso}=\sqrt{{\displaystyle \frac{{\mu }_{SC}}{\rho }}},$$

(27)

In the objective function expressed by Equation (25), the variables are the elastic moduli of minerals, namely the parameters to be inverted. It is assumed that the prior probability density distribution of mineral elastic moduli follows a Gaussian distribution. The mean values of the initial model for mineral elastic moduli are referenced to the rock physics handbook [17]: the initial mean of bulk modulus for sand minerals $K_{sand}$ is 37.00 GPa with a variance of 4.00; the initial mean of shear modulus for sand minerals ${\mu }_{sand}$ is 44.00 Gpa with a variance of 4.00; the initial mean of bulk modulus for clay minerals $K_{clay}$ is 21.00 Gpa with a variance of 4.00; the initial mean of shear modulus for clay minerals ${\mu }_{clay}$ is 7.00 Gpa with a variance of 4.00. The stochastic inversion algorithm with variable sampling space is utilized to invert the elastic moduli of minerals. The preset algorithm parameters are set as: the number of sampling points N=500, maximum iteration number of 10, shrinkage coefficient aqa=0.3 and expansion coefficient apa=1.2. The bulk modulus of sand minerals $K_{sand}=39.24\mathrm{GPa}$, shear modulus of sand minerals ${\mu }_{sand}=43.36\mathrm{GPa}$, bulk modulus of clay minerals $K_{clay}=23.85\mathrm{GPa}$ and shear modulus of clay minerals ${\mu }_{clay}=8.53\mathrm{GPa}$ are acquired after inversion. The inverted elastic moduli of minerals are substituted into Equations (1) and (2) to establish the inversion-corrected rock physics model. The elastic parameters of the isotropic component of the corrected rock physics model in Figure 3 are calculated subsequently. Then, error analysis is conducted between the simulated elastic parameters obtained from the corrected model and the logging elastic parameters, and the corresponding analysis results are presented in Figure 4.

It can be seen from Figure 4 that the maximum frequency of errors of both P- and S-wave velocity appears around zero, which indicates that the inversion process for mineral elastic moduli based on the inversion algorithm with variable sampling space is feasible. At this point, the correction for the isotropic component of the rock physics model is completed.

Taking P-wave velocity, S-wave velocity, density and three anisotropic parameters as known parameters, five physical parameters including pore porosity, pore aspect ratio, oriented fracture porosity, azimuth of oriented fractures and clay content are inverted by using the rock physics model with corrected isotropic component. The inversion results of physical parameters are shown in Figure 5.

Figure 5e compares the inverted clay content with logging-measured clay content, and Figure 5f compares the sum of inverted pore porosity and oriented fracture porosity with logging total porosity. The red curves represent logging data, while the blue curves denote inversion results. It can be seen that both inverted clay content and total porosity match well with logging data. The correlation coefficient of clay content is 0.66, and that of porosity is 0.78, both indicating strong correlations.

Figure 6 shows the comparison between logging data and P-wave velocity, S-wave velocity as well as three anisotropic parameters calculated by the rock physics model using physical parameters inverted from Figure 5.

It can be observed from Figure 6 that the errors of P-wave and S-wave velocities are lower than those of anisotropic parameters. This demonstrates that the inversion accuracy of physical parameters for the isotropic component of the rock physics model (including pore porosity, pore aspect ratio and clay content) is higher than that of physical parameters for the anisotropic component (including aligned fracture porosity and aligned fracture strike azimuth). Meanwhile, it also proves that the inversion of anisotropic physical parameters is more difficult than that of isotropic physical parameters.

3.2. Azimuthal Seismic Reflection Characteristics Analysis Based on Rock Physics

The rock physics model is further combined with the AVAZ reflection coefficient equation to analyze the influences of variations in physical parameters on reflection coefficients and synthetic seismic records. In horizontally transverse isotropic (HTI) media, the reflection coefficient varies with incident angle and azimuth angle [16].

$$\begin{array}{l}{R}_{PP}(\theta ,\phi )={\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathsf{\Delta }IP}{IP}}\\ +{\displaystyle \frac{1}{2}}\left\{{\displaystyle \frac{\mathsf{\Delta }{V}_P}{V_P}}-{\left({\displaystyle \frac{2V_S}{V_P}}\right)}^2{\displaystyle \frac{\mathsf{\Delta }\mu }{\mu }}+\left[\mathsf{\Delta }{\delta }^{(V)}+2{\left({\displaystyle \frac{2V_S}{V_P}}\right)}^2\mathsf{\Delta }{\gamma }^{\left(V\right)}\right]{\cos }^2\left(\phi -{\phi }_{sym}\right)\right\}{\sin }^2\theta \\ +\left\{{\displaystyle \frac{\mathsf{\Delta }{V}_P}{2V_P}}+{\displaystyle \frac{1}{2}}\left[\mathsf{\Delta }{\epsilon }^{(V)}{\sin }^2\left(\phi -{\phi }_{sym}\right)+\mathsf{\Delta }{\delta }^{(V)}{\cos }^2\left(\phi -{\phi }_{sym}\right)\right]{\sin }^2\left(\phi -{\phi }_{sym}\right),\right\}{\sin }^2\theta {\tan }^2\theta \end{array}$$

(28)

In Equation (28) IP stands for P-wave impedance, $\mu$denotes shear modulus, ${\phi }_{sym}$refers to the included angle between the symmetry axis of the medium and the principal direction (generally ${\phi }_{sym}=0^{\circ }$), $\theta$ represents the seismic wave incident angle, and $\phi$ is the observation azimuth angle. The reflection coefficients at different azimuth angles under fixed incident angles $\theta ={30}^{\circ }$ are shown in Figure 7a, and the reflection coefficients at different incident angles under fixed azimuth angles $\phi =0^{\circ }$are shown in Figure 7b:

Figure 7 presents the reflection coefficients at different azimuth angles and incident angles calculated from the logging data in Figure 2 based on Equation (28). It can be observed that the variation of incident angles exerts a more significant influence on reflection coefficients compared with azimuth angles. This phenomenon is caused by the small absolute values of anisotropic parameters derived from logging data. As shown in Figure 2 and Figure 6, the absolute values of anisotropic parameters are relatively small, approximately equal to 0.1.

Since the overall time-domain range of reflection coefficients calculated from time-depth converted logging data in Figure 7 is less than 30 ms, a Ricker wavelet with a dominant frequency of 160 Hz and a sampling interval of 1 ms is adopted to convolve with reflection coefficients to highlight seismic events. Figure 8a shows synthetic seismic records at incident angles of 0°, 15° and 30° when the azimuth angle is 0°; Figure 8b shows synthetic seismic records at incident angles of 0°, 15° and 30° when the azimuth angle is 45°; Figure 8c shows synthetic seismic records at incident angles of 0°, 15° and 30° when the azimuth angle is 90°. Nevertheless, the dominant frequency of actual seismic data can hardly reach 160 Hz. To fit the actual field conditions better, another Ricker wavelet with a dominant frequency of 40 Hz and a sampling interval of 1 ms is used for convolution. Figure 8d displays synthetic seismic records at incident angles of 0°, 15° and 30° under the azimuth angle of 0°; Figure 8e displays synthetic seismic records at incident angles of 0°, 15° and 30° under the azimuth angle of 45°; Figure 8f displays synthetic seismic records at incident angles of 0°, 15° and 30° under the azimuth angle of 90°.

It can be observed from Figure 8a that when the azimuth angle is fixed at 0°, the incident angle imposes an obvious effect on synthetic seismic records, and the amplitude of synthetic seismic records gradually decreases with the increase of incident angle. Figure 8c presents a similar variation law for an azimuth angle of 90°, consistent with that for azimuth angles of 0° and 45°, namely the seismic amplitude declines gradually as the incident angle rises. However, the influence of incident angle on synthetic seismic records is weaker at the azimuth angle of 90° than that at 0° and 45°. Similarly, Figure 8d indicates that incident angle changes have a prominent impact on synthetic seismic records at the azimuth angle of 0°, while Figure 8f shows a fairly weak impact of incident angle at the azimuth angle of 90°. Overall, within the azimuth angle range from 0° to 90°, the smaller the azimuth angle, the more significant the influence of incident angle variation on synthetic seismic records generated by wavelets with different dominant frequencies.

Since the oriented fracture porosity is generally much lower than the pore porosity, the total porosity can be approximately regarded as pore porosity when ignoring the influence of oriented fracture porosity. Under this assumption, we change the pore porosity values and analyzes the effects of pore porosity variations on synthetic seismic records at different incident angles by adopting the inversion-corrected rock physics model. It can be seen from Figure 2 that the total porosity presents a low magnitude, with most values less than 0.1. Keeping other physical parameters unchanged, we discuss the influence of increasing pore porosity on synthetic seismic records at incident angles of 0°, 15° and 30° under the azimuth angle of 0°, and the relevant results are displayed in Figure 9. Consistent with the processing workflow for Figure 8, a Ricker wavelet with a dominant frequency of 160 Hz and a sampling interval of 1 ms is used for convolution of reflection coefficients to highlight seismic events. Figure 9a shows synthetic seismic records at different incident angles calculated by taking the total logging porosity as matrix porosity; Figure 9b,c show synthetic seismic records at different incident angles with assumed matrix porosity of 0.2 and 0.3, respectively. Furthermore, another Ricker wavelet with a dominant frequency of 40 Hz and a sampling interval of 1 ms is adopted for convolution. The impacts of rising matrix porosity on synthetic seismic records at incident angles of 0°, 15° and 30° under the azimuth angle of 0° are illustrated in Figure 9d, Figure 9e and Figure 9f, respectively. Specifically, Figure 9d presents synthetic seismic records at different incident angles calculated using total logging porosity and the 40 Hz dominant-frequency wavelet; Figure 9e corresponds to synthetic seismic records with matrix porosity of 0.2 and the 40 Hz wavelet; Figure 9f corresponds to synthetic seismic records with matrix porosity of 0.3 and the 40 Hz wavelet.

It can be observed from Figure 9 that the amplitudes of synthetic seismic records at different incident angles all increase with the rising pore porosity. The variation law of seismic events obtained from the 160 Hz dominant-frequency wavelet is basically consistent with that from the 40 Hz dominant-frequency wavelet. Compared with synthetic seismic records generated by the 160 Hz wavelet, the records based on the 40 Hz wavelet show a more prominent phenomenon: when the pore porosity is at a low level, the change of pore porosity has a greater influence on seismic amplitude. Specifically, the difference of seismic events between Figure 9e,e is more significant than that between Figure 9e,f.

This section further analyzes the influence of azimuth variation of oriented fractures on synthetic seismic records with different observation azimuth angles. Given that the oriented fracture porosity is assumed to be much smaller than pore porosity during rock physics modeling, the variation of fracture porosity has limited significance for seismic response analysis. In contrast, the change of fracture azimuth imposes a more prominent impact on synthetic seismic records with different observation azimuth angles. In addition, the fractures established in the rock physics model are set to have variable strike azimuths and vertical dips. Theoretically, the variation of fracture azimuth has no influence on seismic waves under vertical incidence (incident angle = 0°), and such influence becomes more significant as the incident angle increases. Therefore, the incident angle is fixed at 30° to explore the effects of varying fracture azimuths on synthetic seismic records at observation azimuth angles of 0°, 45° and 90°, and the corresponding results are shown in Figure 10. To highlight the response difference of seismic events caused by fracture azimuth variation, a Ricker wavelet with a dominant frequency of 40 Hz and a sampling interval of 1 ms is adopted for reflection coefficient convolution. Figure 10a shows synthetic seismic records at different observation azimuth angles calculated based on the predicted fracture azimuth from logging data in Figure 5; Figure 10b, Figure 10c and Figure 10d display synthetic seismic records at different observation azimuth angles with assumed fracture strike azimuths of 0°, 45° and 90°, respectively.

As illustrated in Figure 10, the azimuth of fractures exerts dual effects on synthetic seismic records: it not only modulates the amplitude characteristics of seismic responses, but also changes the spatial position of seismic events. Specifically, fractures with an azimuth of 0° impose a significant influence on synthetic seismic records corresponding to different observation azimuths. Among all tested fracture schemes, the synthetic seismic record established under the fracture azimuth of 45° achieves the optimal matching effect with well logging data, with only slight deviations existing in the position of seismic events. In comparison with other fracture azimuths involved in this study, the most remarkable difference of the 90° strike azimuth lies in its unique amplitude variation law: the amplitude of seismic events presents a continuous decreasing trend with the increase of observation azimuth.

4. Discussion

In the process of rock physics modeling, the influences of pores and oriented fractures on the elastic properties of geological media are considered separately. Variations in pore porosity affect the isotropic elastic parameters of the media, while changes in oriented fractures dominate the variations of anisotropic elastic parameters. The purpose of separately introducing pores and oriented fractures during rock physics modeling is to eliminate the influence of oriented fractures conveniently when correcting the rock physics model via inverted mineral elastic moduli, which can be realized by directly removing the fracture addition procedure from the original rock physics model.

When combining the corrected rock physics model with the approximate reflection coefficient equation of HTI media to investigate how variations in oriented fracture azimuth alter synthetic seismic records, a seismic incident angle of 30° is adopted in numerical simulations under the controlled-variable principle. In this framework, the azimuth of fractures serves as the sole variable, and the rock physics model assumes vertical dip for oriented fractures. At an incident angle of 0° (vertical incidence), changes in fracture azimuth impose no observable effect on seismic responses. In general, the influence of fracture azimuth becomes more pronounced as the incident angle increases. To prominently characterize the azimuthal effect of oriented fractures on synthetic seismic records, we perform simulations at a relatively large incident angle of 30°. This result further implies that high-quality azimuthal seismic data with large incident angles are essential prerequisites for the prediction of steeply dipping fractures.

5. Conclusions

We focus on the critical issue of applying anisotropic rock physics models to analyze logging data and azimuthal seismic data. A rock physics model suitable for logging data is established, and a stochastic inversion algorithm with variable sampling space is pro-posed to invert the parameters of the anisotropic rock physics model for model correction. To improve the adaptability of the rock physics model to actual well logging data, the stochastic inversion algorithm is adopted to invert key parameters required for rock phys-ics modeling, including the elastic moduli of minerals and fracture porosity, so as to cor-rect the rock physics model for better matching with logging measurements. The corrected rock physics model is employed to construct an objective function for petrophysical pa-rameter inversion. The inverted petrophysical parameters are compared with well logging data to verify their accuracy. By combining the corrected rock physics model with the ap-proximate azimuthal reflection coefficient equation and adjusting petrophysical parame-ter values, this study further analyzes the influences of variations in pore porosity and azimuth of oriented fractures on synthetic seismic records. We provide a theoretical foun-dation for monitoring pores and fractures in anisotropic reservoirs containing vertically dipping fractures with arbitrary azimuths using azimuthal seismic data. Simulating the effects of petrophysical parameter variations on synthetic seismic records via theoretical rock physics models enables theoretical characterization of how reservoir petrophysical property changes impact time-lapse seismic data at different azimuths during energy production, which lays a theoretical foundation for monitoring energy exploitation with time-lapse seismic technology.