Pressure-dependent dynamic elastic properties and wave velocities in saturated porous rocks with aligned fractures

SUMMARY

1 INTRODUCTION

Upper crustal rocks are usually permeated with fractures or heterogeneities of various scales (e.g. Guéguen & Kachanov 2011; Liu et al. 2021), which significantly influence the elastic properties of rocks (e.g. Sevostianov & Kachanov 1999; Mear et al. 2007; Xu et al. 2020; Li et al. 2023). Knowledge of the elastic properties of fractured rock is of considerable interest in many areas of research, including rock mechanics, geophysics, underground engineering, petroleum exploration, CO₂ sequestration, mining and the safe storage of nuclear waste (e.g. Hu & Huang 1993; Viete & Ranjith 2006; Zhou et al. 2020; Xu et al. 2021; Wu et al. 2023; Zhang et al. 2023a).

An important issue in geophysical exploration is the elastic characterization of fractured rocks (e.g. Schubnel & Guéguen 2003; Xu et al. 2018). When an elastic wave passes through a fractured rock, a fluid pressure gradient will be created between compliance cracks (microcracks and fractures) and rigid pores (equant pores), causing wave-induced fluid flow (WIFF; Yang et al. 2024). In the case of natural single-phase saturated porous rocks with aligned fractures, there are two main WIFFs including fracture-background WIFF (FB-WIFF; Chapman 2003) and microcrack squirt flow (SF) (Mavko & Nur 1975). To study the FB-WIFF, many models have been proposed based on Biot's poroelastic theory (e.g. Galvin & Gurevich 2009; Gurevich et al. 2010; Fu et al. 2018). Guo et al. (2018a,b) investigated the dynamic (frequency-dependent) elastic properties of FB-WIFF based on the branch function and the full stiffness matrix assumption. In addition to the WIFF mechanism, another important attenuation mechanism is the elastic scattering from macroscopic fractures (e.g. Krenk & Schmidt 1982; Gurevich et al. 1997; Sato et al. 2012; Mavko et al. 2020). By solving the mixed-valued boundary problems of the fracture scattering wavefield in the poroelastic medium, many models have been proposed unified models that consider both the FB-WIFF and elastic scattering mechanisms (e.g. Guo & Gurevich 2020a, b; Song et al. 2020, 2021; Wang et al. 2024b). In addition, Guo et al. (2022a, b) investigated the wave-induced response of intersecting fractures by modifying the boundary conditions and considering the WIFF between intersecting fracture elements (Fracture-Fracture WIFF, FF-WIFF). However, ignoring the irregularity of the background pores is a shortcoming of these models. In this regard, based on previous FB-WIFF and elastic scattering models, Li et al. (2023a) considered the SF mechanism with microcracks in the background medium by modifying Biot's constitutive equations with the Tang et al. (2012) model. Li & Yan (2023) proposed a multiscale model considering FB-WIFF and SF for seismic wavefields based on the full stiffness matrix model (Galvin & Gurevich 2015). Wang et al. (2024a) obtained the effective dynamic stiffness matrix of fractured rocks by solving the scattered wavefields of P and SV waves incident perpendicular to the fracture surface in the pore-microcracked background. However, microcracks in their models have only one shape (with the same aspect ratio). This differs significantly from scanning electron microscopic (SEM) observations of natural porous rocks (e.g. Burns 1985; Nasseri et al. 2007). Therefore, the background needs to be treated as a porous medium containing microcracks of multiple shapes. This model is the basis for pressure-dependent dynamic elastic modelling, and the extension of the background medium to contain multishaped microcracks is the key to this work.

As a factor influencing rock elasticity and velocity, pressure-dependent rock elasticity or wave motion characteristic has always been an important area (e.g. Guéguen & Kachanov 2011; Pimienta et al. 2017). The above studies of arbitrarily complex fractures are unstressed. Especially in engineering, hydrocarbon production can cause changes in reservoir pressure. Moreover, due to wellbore conditions and different reservoir pressure conditions, the influence of formation pressure should be considered when calculating porosity using in-situ acoustic logging data (Khaksar & Griffiths 1999). Under stress loading, the poroelastic rocks will undergo nonlinear elasticity due to the gradual closure of soft pores, hyperelasticity due to residual stress and inelastic deformation due to fracture growth (e.g. Brace et al. 1966; Scholz 1968; Wawersik & Fairhurst 1970; Tapponnier & Brace 1976). All compliant pores will be closed even under small effective pressures, involving nonlinear elasticity. In this aspect, theoretical models have been developed in three aspects: (i) the microstructure-dependent models, such as the dual-porosity model (e.g. Shapiro 2003, 2017; Chen et al. 2023) and the David-Zimmerman model (David & Zimmerman 2012); (ii) The theory of acoustoelasticity, which describes the stress-induced anisotropy of the background medium (e.g. Pao et al. 1984; Meegan et al. 1993; Johnson & Rasolofosaon 1996; Winkler & McGowan 2004; Tromp et al. 2019; Chen & Zong 2022; Chen et al. 2022); (iii) the other theories of nonlinear elasticity (e.g. Cheng 1978; Cheng & Toksöz 1979; Hokstad 2004; TenCate et al. 2016; Delorey et al. 2021).

Based on these theoretical developments, Yan et al. (2014), Wang & Tang (2021) and Tang et al. (2021) proposed the pore structure inversion methods based on laboratory ultrasonic velocities. In addition, Zhang et al. (2022) incorporated ultrasonic velocities and electrical resistivity into the inversion. Ba et al. (2023) inverted the pore aspect ratio and crack radius spectra and obtained the permeability from the pore geometry. Ba et al. (2024) studied the dispersion of rock velocity under stress and frequency by combining acoustoelasticity and squirt flow. For rock containing aligned fractures, Hoenig (1978) obtained the fracture thickness deformation when the fracture axis coincides with the transverse isotropic (TI) background axis. Meanwhile, the background medium is required to be independent of effective pressure (e.g. pure mineral). Yan et al. (2022) modified the model to account for the pressure-dependent background, but the variation of the background parameters (e.g. elastic modulus, porosity and permeability) depends on the measurement at effective pressure loading. However, their model gave the effective static modulus, which did not reflect the influence of WIFF.

In this study, we develop a pressure-dependent dynamic stiffness matrix model of fractured rocks within the hyperelasticity stage. The rocks are composed of equant stiff pores, multishaped microcracks and aligned fractures at different effective pressures. The attenuation mechanisms of FB-WIFF, fracture surface elastic scattering, MMSF and their coupling effects are considered. The pressure-dependent dynamic fracture normal and shear compliances, effective vertical transverse isotropy (VTI) stiffness coefficients and P- and S-wave anisotropy were analysed. Furthermore, the pressure-dependent dynamic wave velocities and attenuation were also investigated. By comparing to the previous laboratory data under effective pressure loading, the model is validated. The present model provides a useful tool for characterizing pressure-dependent elastic properties, wave velocities and attenuation at different frequencies, which will be useful for geological fluid detection, fracture identification, seismic depth conversion and reservoir development. It provides theoretical support for interpreting the field seismic and borehole acoustic logging observations of fractured reservoirs.

2 THEORETICAL FORMULATIONS

2.1 Modelling description

The work aims to build the relationship between effective pressure and multiple attenuation mechanisms in porous fractured rocks, such as FB-WIFF, elastic scattering, MMSF and their coupling effects. For this purpose, a pressure-dependent dynamic stiffness matrix model in fractured rocks is proposed, and the following issues need to be discussed:

• The properties of the background rock under effective pressure loading, in particular, the background modulus: To solve this issue, we followed the process of Wang & Tang (2021) and combined the multishaped microcracks wave motion model [extended from Tang et al. (2012)] with Toksöz et al.’s (1976) crack closure model to obtain the pressure-dependent dynamic background modulus;

• The deformation of inserted aligned fractures under effective pressure loading: This step used the modified fracture deformation model (Hoenig 1978) by Yan et al. (2022). This work further extended it to dynamic due to the dynamic background modulus [see in issue (i)];

• The dynamic stiffness matrix model that considered both microscale randomly oriented microcracks and macroscale aligned fractures (shown in Fig. 1): This step used the dynamic stiffness matrix model with aligned fractures proposed by Wang et al. (2024a), based on Tang et al.’s (2012) modified Biot's model, which includes the FB-WIFF, elastic scattering and SF. However, because compliance microcracks (irregular soft pores) gradually close under effective pressure loading, the background medium of the stiffness model should consider multishaped microcracks.

Figure 1.

A saturated, porous and microcracked background medium permeated by aligned penny-shaped fractures. The microcracks (pore-crack structures) have M group shapes (with different densities and aspect ratios) randomly oriented in the background medium. Therefore, the background medium is still isotropic.

Open in new tabDownload slide

We first extend the stiffness matrix model proposed by Wang et al. (2024a) to a background medium with multishaped coexisting microcracks. Thus, the scattering solutions of the normal incident P and SV waves in the fracture need to be investigated. To do this, a cylindrical coordinate system (rOz) was created for each fracture (d and h denote the diameter and thickness, respectively), as shown in Fig. 2. The fracture centre plane (z = 0) coincides with the centre of the cylindrical system. For simplicity, we consider a P- or SV- wave normal incidence in the x-z plane, where the wave displacements (u_p and u_sv) are:

where |${u_0}$| denotes the displacement amplitude, |${k_1}$| and |${k_3}$| denote the background complex-valued P and S wavenumbers and the time factor |${e^{ - i\omega t}}$| is omitted, with ω and t being the angular frequency and time, respectively.

Figure 2.

Cylindrical coordinate system (rOz) for each fracture with normal incident P or SV wave (u_p or u_sv) in the x-z plane (with φ = 0), where d and h are the fracture diameter and thickness, respectively.

Open in new tabDownload slide

2.2 Dynamic stiffness matrix of microcracked porous rocks with aligned fractures

Under the derivation procedure of Song (2017), through the extended Tang–Biot theory (details in Appendix  A), the fracture general scattering solutions can also be expressed as a superposition of infinitely m-order components. As discussed by Wang et al. (2024a), the scattered wavefields are mainly determined by the zero-order components, so we take the zero-order components (Appendix  B) into the calculation. The difference is that we have incorporated the coupling of the MMSF mechanism into their model.

2.2.1 Dynamic fracture compliances

As investigated by Song et al. (2021) and Wang et al. (2024a), the dynamic saturated fracture normal and shear compliances [|${Z_{n,\rm sat}}( \omega )$| and |${Z_{t,\rm sat}}( \omega )$|] can be obtained from fracture scattered wavefields by Foldy–Lax approximation (Foldy 1945; Lax 1951):

where |$H = {K_d} + {4 {/ {\vphantom {4 3}} } 3}\mu + {\alpha ^2}\bar{M}$| and μ denote the saturated background P-wave low-frequency modulus and shear modulus, respectively, with |${K_d}$| being the modulus of dry background, |$\alpha = 1 - {{{K_d}} {/ {\vphantom {{{K_d}} {{K_s}}}} } {{K_s}}}$| and |${K_s}$| being the Biot's coefficient and the bulk modulus of solid grain, respectively; |${f_P}( \omega )$| and |${f_{SV}}( \omega )$| denote the scattering amplitudes of the normal incident P and SV waves, respectively; and |${n_0} = {{8\varepsilon } {/ {\vphantom {{8\varepsilon } {{d^3}}}} } {{d^3}}}$| denotes the fracture number density, with ε being the fracture density. Particularly, in the low-frequency limit, dry forms of eqs (3) and (4) are consistent with the solutions given by Schoenberg & Douma (1988).

Differently from Li et al. (2024) and Wang et al. (2024a), Biot's modulus |$\bar{M}$| for the case of multishaped microcracks is given by (Tang et al. 2021; Wang & Tang 2021):

where ϕ denotes the rock porosity, |${K_f}$| denotes the fluid modulus and |$\sum\limits_{m = 1}^M {{S_m}( \omega )} $| denotes the summary of the squirting function of M-shaped microcracks.

For each shape (with a different aspect ratio), the microcrack squirting function |${S_m}( \omega )$| can be expressed as:

where η denotes the fluid viscosity, |${\varepsilon _{{\rm mic},m}}$| and |${\gamma _{{\rm mic},m}}$| denote the m-th shape microcrack density and aspect ratio, respectively; |${K_0}$|⁠, |${\nu _0}$| and |${\mu _0}$| denotes the background bulk modulus, Poisson's ratio and shear modulus at |${S_m}( \omega ) = 0$|⁠, respectively, obtained by the Biot-consistent theory (Thomsen 1985) and parameters λ_m and N_m are given by:

Then, the effective background-saturated bulk modulus and shear modulus are given by (Wang & Tang 2021):

where modulus |$K = {K_d} + {\alpha ^2}\bar{M}$| and K(0) denotes the bulk modulus K at zero frequency.

Combining eqs (5)–(9), the background complex-valued wavenumbers (⁠|${k_1}$|⁠, |${k_2}$| and |${k_3}$|⁠) of fast P, slow P and S waves can be obtained through Biot's formulae (specify see Tang et al. 2012, 2021).

2.2.2 Far-field scattering amplitudes

In the Foldy–Lax approximation [eqs (3) and (4)], the far-field scattering amplitudes for a single fracture controlled the fracture normal and shear compliances. For normal incident cases, the far-field scattering amplitudes [|${f_P}( \omega )$| and |${f_{SV}}( \omega )$|] for P and SV waves, written as (Wang et al. 2024a):

where the subscripts ‘P'’ and ‘SV’ indicate the P and SV waves, respectively; |${A_p}( {\omega ,0} )$| and |${C_{sv}}(\omega ,y)$| denote the corresponding amplitude coefficients that are determined by the fracture mixed boundary conditions (the details are given in the next section), respectively. The second variable in the coefficients |${A_p}( {\omega ,0} )$| and |${C_{sv}}(\omega ,y)$| denote the wave incident angle.

2.2.3 Boundary conditions

The P and SV waves incident normal to the fracture surface cause normal symmetric and shear antisymmetric oscillations (e.g. Krenk & Schmidt 1982; Guo & Gurevich 2020a; Wang et al. 2024a), respectively. Therefore, the boundary conditions are given as follows:

For normal incident P wave:

For normal incident SV wave:

where |$\sigma _{zz}^{\rm in}$|⁠, |$\sigma _{zr}^{\rm in}$|⁠, |$\sigma _{z\varphi }^{\rm in}$| and |$p_f^{\rm in}$| denote the scattered normal stress, radial and torsion shear stresses and fluid pressure induced by incidence waves, respectively; |${\sigma _{zz}}$|⁠,|${\sigma _{zr}}$| and |${\sigma _{z\varphi }}$| denote the normal stress and shear stresses, respectively; |${u_z}$|⁠, |${u_r}$| and |${u_\varphi }$| denote the scattered solid displacements in z-, r- and φ- axis directions; |${w_z}$| denotes the scattered fluid displacement relative to the solid in the z-axis direction; and |${p_f}$| denotes the scattered fluid pressure.

2.2.4 Solutions for scattered wavefields

Substitution of the general solutions (details in Appendix  B) into the boundary conditions [eqs (12)–(18)] gives the undetermined coefficients [|${A_p}( {\omega ,0} )$| and |${C_{sv}}(\omega ,y)$| in eqs (10) and (11)] for normal incident P and SV waves, respectively. The procedures for solving the coefficients are consistent with Wang et al. (2024a). For simplicity, we provided the derivation in Appendices  C and  D for the P and SV waves, respectively. For the details, please refer to the work of Wang et al. (2024a).

2.3 Pressure effects in microcracked porous rocks with aligned fractures

Once effective pressure on the fractured rock, it will simultaneously compress the rock background and the aligned fractures (e.g. Han et al. 2022; Yan et al. 2022). Especially for impure mineral rock backgrounds (e.g. sandstone and limestone), the increase in effective pressure will lead to an increase in the elasticity of the rock background (e.g. Coyner 1984; Yin et al. 2017; Zhang et al. 2022; Ba et al. 2023). Therefore, it is necessary to take into account both the pressure-dependent background properties and the fracture deformation.

2.3.1 Pressure-dependent dynamic background properties

Following the preview studies (e.g. David & Zimmerman 2012; Wang & Tang 2021; Zhang et al. 2022), the internal cause for the increase in elasticity of the sandstone is the decrease in microcracks (soft pores) under effective pressure loading. To complete the modelling, we make the following assumptions:

• Under the initial effective pressure condition (⁠|${P_0}$| = 0), there are M shapes/types of microcracks in the rock background, which are randomly oriented and distributed as shown in Fig. 1;

• Assume there are N calculation points in the calculated effective pressure loading interval, and the effective pressure at each calculation point is expressed as |${P_n}$|⁠. By doing this, the background properties in each small effective pressure loading interval can be considered to be consistent;

• For each effective pressure point P_n, the m-th shaped unclosed microcrack density and aspect ratio are expressed as |${\varepsilon _{{\rm mic},nm}}$| and |${\gamma _{{\rm mic},nm}}$|⁠, respectively.

According to Toksöz et al. (1976) and Wang & Tang (2021), the microcrack porosity change rate for m-th shape microcracks under P_n relative to P₀ is given by:

Assuming that the long axis of the microcrack is independent of P, therefore, the m-th microcrack aspect ratio change rate is equal to that for microcrack porosity. The microcrack aspect ratio under P_n can be expressed as:

The microcrack closure condition is satisfied:

In addition, the initial microcrack distribution in unfractured can be obtained rock by mercury injection, nuclear magnetic resonance (NMR), SEM and laboratory ultrasonic P and S wave velocity inversion.

2.3.2 Pressure-dependent dynamic fracture deformation

As discussed by Hoenig (1978) and Yan et al. (2022), the fracture thickness deformation U (also assuming the fracture diameter is independent of P) for a single fracture can be expressed as follows:

where the fracture deformation intensity |$L( P )$| in an isotropic background (due to the randomly oriented microcracks) is given by:

and

with modulus |${L_a} = H - 2\mu $|⁠.

Similarly, if the fracture thickness deformation (U) is greater than or equal to the fracture's initial thickness (h), the fracture will close:

Finally, based on the non-interaction approximation (NIA) theory (e.g. Grechka & Kachanov 2006; Li et al. 2023), the fracture deformation of the aligned fracture set is consistent with the single fracture. Therefore, the fracture porosity at P_n is given by:

It should be noted that the original Hoenig (1978) model was presented for the case of pressure-independent background. The model describes the compressive deformation of the microcracks, which was characterized by a constant background stiffness parameter. This work considers a porous background containing a set of aligned penny-shaped fractures. Therefore, the extension of Hoenig's model by Yan et al. (2022) for aligned fractures is used, where the fracture deformation is determined by the integral of the fracture deformation during the pressurization process. The background modulus is different at each point, which is determined by the closure of the microcracks and is a nonlinear changing process.

2.4 Dynamic elastic properties of fractured rock

2.4.1 Dynamic fracture compliance tensor

Eqs (3) and (4) obtained the dynamic saturated fracture normal and shear compliances. For the case where the symmetry axis of the aligned fractures corresponds to that of the background medium, the saturated fracture excess compliance tensor Z has the following form (Guo et al. 2024):

2.4.2 Dynamic stiffness tensor and attenuation factor

Through linear slip theory (Schoenberg 1980), the effective saturated stiffness tensor C can be easily obtained, as follows:

where |${{{\bf C}}_{{{\bf bk}}}}$| denotes saturated background stiffness tensor.

In addition, the stiffness coefficients |${C_{ij}}( {i,j = 1,2,3,...,6} )$| obtained by eq. (28) are also complex valued. Therefore, the stiffness attenuation factor is given by (Collect & Gurevich 2016):

where ‘Im(∗)’ and ‘Re(∗)’ denote the imaginary and real parts, respectively.

2.5 Phase wave velocities and attenuation

According to the Christoffel equation (Auld 1990), the complex-valued wave velocities (details in Appendix  E) can be obtained by complex-valued stiffness coefficients. Finally, the phase wave velocities and attenuation can be calculated through the complex-valued velocities (Carcione 2007):

where |${\bar{V}_i}( {i = P,SV,SH} )$| denotes the velocities for P, SV and SH waves, respectively.

3 NUMERICAL RESULTS

3.1 Sample parameters

To approximate the natural rock, we chose Weber sandstone as the background medium. According to Burns (1985), the Weber sandstone is a type of fine-grained sandstone containing clay and calcite, accompanied by a large number of microcracks. Following the laboratory data from Coyner (1984), the Weber sandstone is saturated with benzene. The porosity and density of the rock frame are 9.5 per cent and 2392 kg m⁻³, respectively. The density, bulk modulus and viscosity of benzene are 880 kg m⁻³, 1.21 GPa and 0.647 × 10⁻³ Pa s, respectively. And the static permeability is 63 × 10⁻¹⁵ m². Tang et al. (2021) inverted the pore structure of the Weber sandstone and found that there are microcracks with densities of 0.0032, 0.0084, 0.01276, 0.006, 0.01453, 0.01681, 0.01028, 0.00863, 0.01489, 0.01805, 0.00936, 0.01123, 0.01433, 0.01432, 0.0246, 0.02456, 0.0227, 0.03456 and 0.01817, corresponding to aspect ratios of 5.83 × 10⁻³, 4.87 × 10⁻³, 3.91 × 10⁻³, 2.74 × 10⁻³, 2.25 × 10⁻³, 1.88 × 10⁻³, 1.48 × 10⁻³, 1.38 × 10⁻³, 1.27 × 10⁻³, 1.16 × 10⁻³, 1.08 × 10⁻³, 9.6 × 10⁻⁴, 8.6 × 10⁻⁴, 7.9 × 10⁻⁴, 6.9 × 10⁻⁴, 6.4 × 10⁻⁴, 4.9 × 10⁻⁴, 2.5 × 10⁻⁴ and 2.0 × 10⁻⁴. In addition, the solid bulk and shear modulus are 28 and 26.9 GPa, respectively.

Treating the Weber Sandstone as a background medium, we hypothesize that a group of aligned penny-shaped fractures permeated it. The fractures are hydraulically in communication with the porous background, allowing fluid flow between the fractures and the background. Meanwhile, the fracture diameter and thickness (under the initial effective pressure) are set to 1 and 0.01 m, respectively, and the fracture density to 0.05.

3.2 Pressure influence on dynamic fracture compliances

Under effective pressure increasing, both the background medium and the fracture are compressed. Fig. 3 shows the real and imaginary parts of the normal (solid lines) and shear (dashed lines) fracture compliances, with the real part representing the elastic properties and the imaginary part for the attenuation. It can be seen that both the real parts of the fracture normal and shear compliances decrease (stiffness increase) with increasing effective pressure at frequencies lower than 10⁴ Hz. This is mainly due to the hardening of the background medium as microcracks close under effective pressure. Otherwise, for the real part of the fracture normal compliance, the relaxation amplitude of both FB-WIFF and elastic scattering mechanisms decreases with effective pressure increasing. The imaginary part also shows a significant decrease. Especially for the elastic scattering mechanism (around 10³–10⁴ Hz), the increase in effective pressure causes a decrease in the amplitude of the elastic scattering mechanism. This is mainly because of the decrease in fracture thickness under effective pressure. Simultaneously within this frequency band, effective pressure has a greater effect on fracture shear compliance than normal compliance.

Figure 3.

Real (a) and imaginary (b) parts of fracture normal (solid lines) and shear (dashed lines) compliances [eqs (3) and (4)] with frequency under effective pressures of 5, 10, 15, 20 and 90 MPa. ‘FB-WIFF’ and ‘Elastic scattering’ refer to the fracture-background WIFF and elastic scattering, respectively.

Open in new tabDownload slide

In addition, at the high-frequency limit (above 10⁴ Hz), both the fracture normal and shear compliances are close to zero. This is because the fractures are hydraulically separated from the background and the frequencies cross the scattering domain (wavelength comparable to the fracture size) in these frequencies. Moreover, due to the low fracture density, the elasticity of the rock is now mainly controlled by the background medium. In other words, the effect of effective pressure on fractures at the high-frequency limit can be ignored and the model is degenerated into modelling the compression of porous rocks (e.g. Wang & Tang 2021).

3.3 Pressure influence on dynamic stiffness coefficients

Figs 4–8 show the dynamic stiffness coefficients [C₁₁, C₁₃, C₃₃, C₄₄ and C₆₆] at effective pressures of 5, 10, 15, 20 and 90 MPa respectively, both real part stiffness coefficients and attenuation factors [eq. (29)] are calculated. For stiffness coefficients C₁₁ and C₆₆ (Figs 4 and 8), these two stiffness coefficients are almost unaffected by the fracture-related attenuation mechanisms (FB-WIFF and elastic scattering), especially for stiffness coefficient C₆₆. Therefore, the relaxation of stiffness coefficients C₁₁ and C₆₆ are mainly controlled by the background MMSF mechanism and the MMSF coupling effect. Due to the gradual closure of multishaped microcracks under effective pressure loading, the MMSF mechanism in the background medium gradually weakens. Similar phenomena also occur with other coefficients. When the effective pressure reaches 90 MPa, the microcracks in the background medium are almost completely closed. Therefore, the real part of stiffness coefficients C₁₁ and C₆₆ are almost independent of frequency. This indicates that the velocities of horizontally propagating (angle between incident direction and the z-axis is 90°) P wave (e.g. Guo & Gurevich 2020b) and S wave (e.g. SH-wave dispersion curves, see Fig. 16) are almost non-dispersive with frequency.

Figure 4.

Real part (a) and attenuation factor (b) of stiffness coefficient C₁₁ with frequency under effective pressures of 5, 10, 15, 20 and 90 MPa. ‘MMSF’ and ‘Coupling of MMSF’ refer to ‘multishaped microcracks' squirt flows and the coupling effects, respectively.

Open in new tabDownload slide

Figure 5.

Real part (a) and attenuation factor (b) of stiffness coefficient C₁₃ with frequency under effective pressures of 5, 10, 15, 20 and 90 MPa. ‘FB-WIFF’, ‘Elastic scattering’ and ‘MMSF’ refer to fracture-background WIFF, elastic scattering and multishaped microcracks' squirt flows, respectively.

Open in new tabDownload slide

Figure 6.

Real part (a) and attenuation factor (b) of stiffness coefficient C₃₃ with frequency under effective pressures of 5, 10, 15, 20 and 90 MPa. ‘FB-WIFF’, ‘Elastic scattering’ and ‘MMSF’ refer to fracture-background WIFF, elastic scattering and multishaped microcracks' squirt flows, respectively.

Open in new tabDownload slide

Figure 7.

Real part (a) and attenuation factor (b) of stiffness coefficient C₄₄ with frequency under effective pressures of 5, 10, 15, 20 and 90 MPa. ‘Elastic scattering’, ‘MMSF’ and ‘Coupling of MMSF’ refer to elastic scattering, multi-shaped microcracks' squirt flows and the coupling effects, respectively.

Open in new tabDownload slide

Figure 8.

Real part (a) and attenuation factor (b) of stiffness coefficient C₆₆ with frequency under effective pressures of 5, 10, 15, 20 and 90 MPa. ‘MMSF’ and ‘Coupling of MMSF’ refer to multishaped microcracks' squirt flows and the coupling effects, respectively.

Open in new tabDownload slide

For the stiffness coefficient C₁₃ (Fig. 5), the effects of the three attenuation mechanisms can be observed. Under low effective pressure conditions (5, 10, 15 and 20 MPa), the effect of effective pressure on the MMSF mechanism is greater than that on the FB-WIFF and elastic scattering. This indicates that effective pressure has a greater effect on the microcracks in the background medium than on the fractures for stiffness coefficient C₁₃. However, at higher effective pressures (90 MPa), the background microcracks are almost completely closed. In this case, the deformation of the fractures becomes the main factor. In addition, the FB-WIFF mechanism is more sensitive to effective pressure than the elastic scattering mechanism.

For stiffness coefficient C₃₃ (Fig. 6), the real part increases significantly with the effective pressure compared to C₁₃. The other variation characteristics at effective pressure are similar to C₁₃. In addition, for normal S-wave stiffness coefficient C₄₄ (Fig. 7), only elastic scattering, MMSF and MMSF coupling effects influence the dynamic modulus. The increase in effective pressure only increases the amplitude of the elastic scattering mechanism effect on the real part of C₄₄ but does not affect its attenuation. This is because the elastic scattering mechanism in the normally propagating S wave (e.g. Song et al. 2021) is largely unaffected by the thickness reduction (effective pressure increase). The simulation results of dynamic stiffness coefficients C₄₄ and C₆₆ indicate that the S-wave anisotropy parameter (Thomsen 1986) is controlled only by elastic scattering and is almost independent of the effective pressure. This phenomenon is discussed in detail in the next section.

3.4 Pressure influence on dynamic anisotropy

A major advantage of obtaining the effective stiffness of fracture-induced effective vertical transverse isotropic (VTI) rocks is the ability to simulate anisotropy parameters directly. Following the anisotropic parameters defined by Thomsen (1986), the two P-wave anisotropy parameters and one S-wave anisotropy parameter are given by:

where |${\varepsilon _{\rm TH}}$| and |${\gamma _{\rm TH}}$| describe the differences between horizontal incident P and S waves and vertical incident, respectively; and |${\delta _{\rm TH}}$| describes the variation in directions of P waves near the fracture normal.

Due to the complex-valued stiffness coefficients, the anisotropy parameters are also complex, with the real part representing the wave velocity anisotropy and the imaginary part representing the attenuation anisotropy. Fig. 9(a) shows the P-wave velocity anisotropic parameters, the solid and dashed lines denote the parameters |${\varepsilon _{\rm TH}}$| and |${\delta _{\rm TH}}$|⁠, respectively. In the low-frequency range limit of FB-WIFF, the effective pressure does not affect the real part of the anisotropic parameters. However, at the higher frequencies of FB-WIFF (10⁰–10⁴ Hz), both |${\varepsilon _{\rm TH}}$| and |${\delta _{\rm TH}}$| decrease with effective pressure increase. Ultimately, this is due to the decrease in crack thickness, which is consistent with the analysis of Guo & Gurevich (2020b). In the elastic scattering regime, the P-wave velocity anisotropy decreases with effective pressure increasing. Furthermore, effective pressure increasing has a relatively small effect on S-wave velocity anisotropy (Fig. 9b), only the elastic scattering mechanism is included but no FB-WIFF mechanism occurs. This is because horizontally and vertically propagating S waves do not undergo the FB-WIFF mechanism (e.g. Song et al. 2021; Li et al. 2024; Wang et al. 2024b). For the attenuation anisotropy (Fig. 10), effective pressure increasing in the FB-WIFF regime reduces the attenuation anisotropy, while that in the elastic scattering regime increases the attenuation anisotropy.

Figure 9.

Real part of P (a) [|${\varepsilon _{\rm TH}}$| (solid lines, labelled as Re(⁠|${\varepsilon _{\rm TH}}$|⁠)) and |${\delta _{\rm TH}}$| (dashed lines, labelled as Re(⁠|${\delta _{\rm TH}}$|⁠))] and S (b) [|${\gamma _{\rm TH}}$| (solid lines, labelled as Re(⁠|${\gamma _{\rm TH}}$|⁠))] wave anisotropy parameters with frequency under effective pressures of 20, 40, 60, 80 and 100 MPa. Note that the real part indicates the velocity anisotropy.

Open in new tabDownload slide

Figure 10.

Imaginary part of P (a) [|${\varepsilon _{\rm TH}}$| (solid lines, labelled as Im(⁠|${\varepsilon _{\rm TH}}$|⁠)) and |${\delta _{\rm TH}}$| (dashed lines, labelled as Im(⁠|${\delta _{\rm TH}}$|⁠))] and S (b) [|${\gamma _{\rm TH}}$| (solid lines, labelled as Im(⁠|${\gamma _{\rm TH}}$|⁠))] anisotropy parameters with frequency under effective pressures of 20, 40, 60, 80 and 100 MPa. Note that the imaginary part indicates the attenuation anisotropy.

Open in new tabDownload slide

For P wave, attenuation anisotropy is more sensitive to effective pressure than velocity anisotropy (compare Figs 9a and 10a). However, for S-wave anisotropy, velocity anisotropy seems to be more sensitive than attenuation anisotropy (compare Figs 9b and 10b). In addition, although effective pressure increases directly reduce fracture thickness, variations in background elasticity also influence fracture-related mechanisms. In short, although isotropic backgrounds do not produce anisotropy, they can influence the geometry of the fracture and thus the anisotropy of the fractured rock.

3.5 Pressure influence on wave velocity and attenuation

The dynamic stiffness coefficients can be used to obtain the wave velocities and attenuation [via eqs (E1)–(E3) and (30)–(31)] in the y-z plane. Fig. 11 shows the normal incidence of P-wave velocity and attenuation with frequency at effective pressures of 5, 10, 15, 20 and 90 MPa (solid lines). For comparison, the background P-wave velocity and attenuation calculated by Tang et al. (2021) model are also included (dashed lines). Three types of attenuation mechanisms can be seen, including FB-WIFF, fracture elastic scattering, MMSF, coupling of MMSF and their coupling effects. The P-wave velocity increases as the effective pressure increases, mainly due to the compression of the background medium. The reduction deformation of the fracture thickness mainly influences the characteristic frequencies of these mechanisms. Due to the lower fracture density, the reduction in fracture porosity caused by effective pressure increase has a negligible effect on the velocity. As the microcracks close under effective pressure, the MMSF mechanism gradually disappears. This can also be seen in terms of velocity and attenuation, especially between 10⁵ and 10⁶ Hz. In addition, we also show the P-wave velocity measured by Coyner (1984) (red dots) for the Weber sandstone at 10⁶ Hz under effective pressure increasing. The good prediction compared to the measured data further validates the model. When the incidence angle is increased to 45° (Fig. 12), it can be seen that the variations of fracture-related mechanisms affect their coupling effects with MMSF. The incidence angle is also a significant factor affecting wave propagation under effective pressure increase.

Figure 11.

Variations of normal incident P-wave velocity and attenuation with frequency at effective pressures of 5, 10, 15, 20 and 90 MPa. ‘Fractured’ (solid lines) and ‘Unfractured’ (dashed lines) denote our model results and the P-wave results of the porous multishape microcracked background medium calculated by the Tang et al. (2021) model, respectively. In addition, the red dots are the Weber sandstone P-wave velocity data measured by Coyner (1984). ‘FB-WIFF’, ‘Elastic scattering’, ‘MMSF’ and ‘Coupling of MMSF’ refer to fracture-background WIFF, elastic scattering, multishaped microcracks' squirt flows and the coupling effects, respectively.

Open in new tabDownload slide

Figure 12.

Variations of 45° incident P-wave velocity and attenuation with frequency at effective pressures of 5, 10, 15, 20 and 90 MPa. ‘Fractured’ (solid lines) and ‘Unfractured’ (dashed lines) denote our model results and the P-wave results of the porous multishape microcracked background medium calculated by the Tang et al. (2021) model, respectively. ‘FB-WIFF’, ‘Elastic scattering’ and ‘MMSF’ refer to fracture-background WIFF, elastic scattering and multishaped microcracks' squirt flows, respectively.

Open in new tabDownload slide

Fig. 13 shows the results for the normal incident SV wave (solid lines), the background Weber sandstone S-wave velocity measured by Coyner (1984) (red dots) is also shown. The background S-wave velocity and attenuation calculated by Tang et al. (2021) model are dashed lines. Because the vertical propagation SV wave velocity is controlled by shear modulus C₄₄ (Fig. 7). Therefore, the velocity and attenuation characteristics are the same for the C₄₄. In addition, the good agreement between the S-wave velocity of the Weber sandstone measured data and the model results indicates the accuracy of our model in predicting shear waves. Compared to the P-wave velocity prediction, the difference between the model prediction and the measured S-wave data can be attributed to the microcracks shear stiffening effect and the clay shear weakening effect. Similarly, we also calculated the incidence of SV waves at 45° (Fig. 14), where the maximum FB-WIFF mechanism controlled the wave dispersion in the seismic frequencies.

Figure 13.

Variations of normal incident SV-wave velocity and attenuation with frequency at effective pressures of 5, 10, 15, 20 and 90 MPa. ‘Fractured’ (solid lines) and ‘Unfractured’ (dashed lines) denote our model results and the P-wave results of the porous multishape microcracked background medium calculated by the Tang et al. (2021) model, respectively. In addition, the red dots are the Weber sandstone S-wave velocity data measured by Coyner (1984). ‘Elastic scattering’, ‘MMSF’ and ‘Coupling of MMSF’ refer to elastic scattering, multishaped microcracks' squirt flows and the coupling effects, respectively.

Open in new tabDownload slide

Figure 14.

Variations of 45° incident SV-wave velocity and attenuation with frequency at effective pressures of 5, 10, 15, 20 and 90 MPa. ‘Fractured’ (solid lines) and ‘Unfractured’ (dashed lines) denote our model results and the P-wave results of the porous multishape microcracked background medium calculated by the Tang et al. (2021) model, respectively. ‘FB-WIFF’, ‘Elastic scattering’ and ‘MMSF’ refer to fracture-background WIFF, elastic scattering and multishaped microcracks' squirt flows, respectively.

Open in new tabDownload slide

In addition, Figs 15 and 16 show the numerical results for the SH wave at an incidence angle of 45° and 90°, respectively. We also compared the model-predicted S-wave velocity at different effective pressures with measured S-wave data (red dots), which showed good agreement. In contrast to the 45° incident P and SV waves, the SH wave is not affected by the FB-WIFF mechanism. This is because the SH-wave velocity and attenuation are controlled by modulus C₄₄ and C₆₆ (Figs 8 and 7). Therefore, for SH wave, the effective pressure is independent of the FB-WIFF attenuation mechanism (e.g. Guo et al. 2018a; Wang et al. 2024a). In particular, when the incidence angle is 90° (Fig. 16), the effect of the effective pressure on the fractured rock SH wave can be equivalent to the effect on the background medium.

Figure 15.

Variations of 45° incident SH-wave velocity and attenuation with frequency at effective pressures of 5, 10, 15, 20 and 90 MPa. ‘Fractured’ (solid lines) and ‘Unfractured’ (dashed lines) denote our model results and the P-wave results of the porous multishape microcracked background medium calculated by the Tang et al. (2021) model, respectively. In addition, the red dots are the Weber sandstone S-wave velocity data measured by Coyner (1984). ‘Elastic scattering’, ‘MMSF’ and ‘Coupling of MMSF’ refer to elastic scattering, multishaped microcracks' squirt flows and the coupling effects, respectively.

Open in new tabDownload slide

Figure 16.

Variations of horizontally incident SH-wave velocity and attenuation with frequency at effective pressures of 5, 10, 15, 20 and 90 MPa. ‘Fractured’ (solid lines) and ‘Unfractured’ (dashed lines) denote our model results and the P-wave results of the porous multishape microcracked background medium calculated by the Tang et al. (2021) model, respectively. ‘MMSF’ and ‘Coupling of MMSF’ refer to multishaped microcracks' squirt flows and the coupling effects, respectively.

Open in new tabDownload slide

4 COMPARISON OF MODEL PREDICTION AND LABORATORY DATA

To verify the model, the model predicted P-, SV- and SH-wave velocities under effective pressure loading (5–50 MPa) are compared with laboratory ultrasonic (500 kHz) data measured by Han et al. (2022).

4.1 Inversion of pore structure of artificial brine-saturated rock without fracture

According to the modelling, to predict the velocity of fractured rock under effective pressure loading, it is necessary to know the pore structure of unfractured background rock. Therefore, the unfractured artificial brine-saturated sandstone P and S wave velocities measured by effective pressure loading are first used to inverse the background pore structure. The inversion follows Wang & Tang's (2021) process, and since only saturation data are reported, the following objective function is used:

where is the stiff pore porosity; |$V_{p,n}^{\rm {Meas}}$| and |$V_{s,n}^{\rm {Meas}}$| denote the measured P and S wave velocities under effective pressure P_n; |$V_{p,n}^{\rm {Model}}$| and |$V_{s,n}^{\rm {Model}}$| denote the modelled P and S wave velocities; the subscript ‘sat’ indicates saturation.

Fig. 17 shows the inversion results of the unfractured sample. In Fig. 17(a), the scatters and lines are measurement velocity data and model fitting, respectively. The good agreement between the two results indicates that the inversion results are accurate. The inverted pore structure spectrum is plotted on a double logarithmic scale where the blue and red bars represent stiff pores and microcracks (soft pores), respectively, as shown in Fig. 17(b). The inverted microcrack aspect ratio is distributed between 7 × 10⁻⁴ and 10⁻². Fig. 17(c) shows the model porosities (red star line) using the inversion results with effective pressure loading that is close to the measured results (scattered points). The total microcrack density decreases with effective pressure loading (Fig. 17d), and the closure of microcracks is found under almost every effective pressure, which further verifies that the inversion of the pore structure is accurate. In addition, the other two inverted parameters bulk and shear moduli of the grain are K_s = 26.98 GPa and μ_s = 24.48 GPa, respectively. Furthermore, the resistivity measurements can be combined with the corresponding model to improve the inversion results (e.g. Zhang et al. 2022).

Figure 17.

Inversion results for unfractured sample: (a) comparison of laboratory-measured (Han et al. 2022) velocity data (scattered points) and model fitting (lines) by inversion results with effective pressure loading; (b) inverted pore structure spectrum including microcracks (soft pores) and stiff pores; (c) comparison of laboratory-measured porosity (scattered points) and model fitting (red star line) by inversion results with effective pressure loading; (d) variations of total microcrack density with effective pressure loading.

Open in new tabDownload slide

4.2 Comparison with laboratory velocity data in the fractured sample

Based on the inversion of the pore structure of the unfractured sample, it is further used to predict the P-, SV- and SH-wave velocities of the fractured sample under three different angles of incidence (0°, 45° and 90°, respectively) and effective pressure loading. The fracture density is 0.062 with fracture diameter and thickness under initial zero effective pressure are 3 and 0.063 mm, respectively. Fig. 18 compares model predictions (lines) with laboratory measurements (scattered points). For the P wave, the model has a poor prediction for 0° and 45° incidences, but a good prediction for 90° incidences. For both SV- and SH-shear waves, the prediction results of the model are close to the laboratory results, especially for SV waves.

Figure 18.

Prediction results for fractured sample: predicted P- and SV-wave velocities with effective pressure loading for the angles of incidence are 0° (a), 45° (b) and 90° (c); predicted SH-wave velocities with effective pressure loading for the same three angles. The scattered points and lines denote the measured (Han et al. 2022) and modelled results, respectively.

Open in new tabDownload slide

There are two reasons for the difference between the predicted results of the model and the measured results. One is that the synthetic samples use stratification technology, and the other is that there is a certain difference between the unfractured sample and the background of fractured samples. In all, the model's accuracy is verified by the successful prediction of laboratory data.

5 DISCUSSIONS

5.1 Advantages and application value

In this work, we developed a theoretical modelling approach for the pressure-dependent dynamic elastic properties and wave velocities of porous rocks with aligned penny-shaped fractures. By considering the pressure response of background properties and fracture deformation in combination with a scattering-based fractured rock stiffness model. The elastic properties of fractured porous rocks under varying effective pressures can be determined. The advantage of our model is that it establishes the relationship between fractured rocks' effective pressure and dynamic attenuation mechanisms, including fracture-background WIFF, fracture elastic scattering and microcrack squirt flow. Hence, our model connects the data measured at different frequencies [seismic exploration (∼Hz), acoustic logging (∼kHz) and laboratory ultrasonic measurement (∼MHz)] under different in-suit pressures. This helps to integrate different measurements for in-situ exploration of fractured reservoirs. Simultaneously, in the natural case, the microstructures of host rocks for fractures can be treated as a pore aspect ratio spectrum. The closure of microcracks (irregular soft pores) is the main expression of compression in the host rock. Therefore, the fractured rock models constructed against the background of pores and single-shape pore microcrack structures cannot be directly used for effective pressure modelling (e.g. Li et al. 2024; Wang et al. 2024a). To solve this problem, our model extends the background of fractured rock to consider multishaped microcracks. This makes our model superior in simulating the elastic and wave properties of fractured rock, as the pore structure of the real background rock can be considered. The extended dynamic full stiffness model completely degenerates to Wang et al.’s (2024a) model when considering the single-shaped microcracks. Although this work is based on Weber sandstone, due to the universality of background wave motion theory (Tang et al. 2012, 2021; Wang & Tang 2021), this model can be applied to arbitrarily porous rocks (from large porosity sandstone to tight limestone). This is because the multishaped microcracks can be used to model more general and complex pore structures. In addition, when the complex pore structure acts together on the pressure, the coupling between different mechanisms will cause the complex attenuation behaviours of the wave propagation. In this condition, our model can accurately describe the wave velocity and attenuation and thus is more suitable for application in such reservoirs. Furthermore, by combining our model with the electrical model for fractured rocks with porous backgrounds (e.g. Hatta & Taya 1986; Shafiro & Kachanov 2000), it is possible to invert the rock properties by combining seismic and electrical data for the fractured reservoirs. This would benefit several geoscience fields, such as CCUS (Kumar et al. 2023) and the exploration and development of unconventional reservoirs (e.g. Burton & Wood 2013; Gale et al. 2014).

5.2 Assumptions and limitations

Although our model constructs the effective pressure response relationship of the elastic properties of fractured rock, there are still limitations. This model assumes that the aligned fractures are sparsely distributed in the porous background. The background consists of pores and randomly oriented, multishaped microcracks, and the fractures are considered to be connected to the porous background. Pores, microcracks and fractures are saturated with single-phase fluid. Isotropic/uniform confining and pore pressures are applied to such rocks. In fracture deformation modelling, the use of the NIA assumption requires that fractures are sparsely distributed which caused our model to be valid for the case of low fracture density—neglecting the effects of stress amplification and shielding between the fractures (e.g. Grechka & Kachanov 2006; Zhang et al. 2023b). This is consistent with the assumption of the dynamic full-stiffness model for fractured rocks (Wang et al. 2024a). In addition, for layered fractured rocks, the pressure response modelling of the background medium must consider intrinsic anisotropy. The interlayer wave-induced fluid flow needs to be considered (e.g. Krzikalla & Müller 2011; Liao et al. 2023; He & Guo 2024). For fracture deformation, the Hoenig (1978) model can be used in the TI background (Yan et al. 2022). The pressure-dependent TI rock elastic properties can be established by introducing microcracks into existing models (e.g. Sayers & Kachanov 1995) and establishing a layered rock model with microcracks as reference variables. Due to the isotropic confining and pore pressures assumed to be applied to the rock, the closure of microcracks is not aligned. When uniaxial, multiaxial or off-axis stress is applied, the microcracks that are perpendicular or nearly perpendicular to the stress axis will close, resulting in aligned microcracks (producing anisotropy). In those cases, our modelling approach will not apply. Take the uniaxial stress as an example, the elasticity of the aligned closure of the background microcracks is the key to modelling. For fracture deformation, the uniaxial stress can be decomposed into its normal and tangential directions, and the stress component along the fracture normal is the major contributor to fracture deformation. Meanwhile, the effects of stress and the elasticity of the background rock must be considered to obtain the background (anisotropy modulus) and fracture parameters under uniaxial stress. Finally, the elastic properties or wave propagation characteristics are studied by combining them with the fractured rock model considering the anisotropy background (e.g. Xu et al. 2015, Xu et al. 2020; Li et al. 2023; 2024). In addition, the proposed model is limited within the hyper-elasticity stage where further compressive deformation of stiff pores is not considered. Stress-induced anisotropy is also not considered (e.g. Fu & Fu 2017). We intend to investigate these limits shortly.

6 CONCLUSIONS

A model to simulate the pressure-dependent dynamic elastic properties of fractured porous rocks has been developed in this work. The effects of pressure-dependent background dynamic properties and fracture deformation have been investigated. In addition, we have also introduced the MMSF. Using this model, we investigated dynamic fracture compliances, dynamic stiffness, P- and S-wave anisotropy, P-, SV- and SH-wave velocity and attenuation under effective pressure increasing. The results indicated that pressure significantly reduces fracture compliances and increases rock stiffness. However, for fracture compliances, the influence of pressure in the ultrasonic frequency range can be ignored. Pressure affects the attenuation anisotropy parameters more than the velocity anisotropy parameters. In addition, the S-wave anisotropy parameter included only the fracture elastic scattering mechanism and was almost independent of effective pressure. By comparing the model predictions with the previous laboratory wave velocities measured at ultrasonic frequency, we verified our model. It should be noted that the current validation is based on a limited number of laboratory samples. Further validation will require additional dynamic experimental and field data, which we plan to investigate in the near future.

ACKNOWLEDGEMENTS

This study was supported by the Open Research Fund of State Key Laboratory of Deep Oil and Gas (No. SKLDOG2024-KFZD-03, SKLDOG2024-ZYTS-04), the National Natural Science Foundation of China (grant no. 42174145, 42274146, 41821002), Shandong Provincial Natural Science Foundation (Grant No. ZR2024YQ062), Laoshan National Laboratory Science and Technology Innovation Project (grant no. LSKJ202203407)).

DATA AVAILABILITY

This is a theoretical paper. The data in this paper are available on reasonable request from the corresponding author. The Weber sandstone laboratory data are available from Coyner (1984). The Weber sandstone pore microstructure spectrum inverted from velocity data is available from Tang et al. (2021). The measured velocity data and details for unfractured and fractured artificial samples can be found in Han et al. (2022).

REFERENCES

APPENDIX A: WAVE GOVERNING EQUATIONS

By summing the multishaped microcrack squirt flow terms in the constitutive equations of the Tang et al. (2012) model, Tang et al. (2021) proposed the wave motion theory of a porous microcracked medium with multiple morphologies. The governing equations can be expressed as follows:

where σ, |${{{\bf u}}_{{\bf s}}}$| and w denote the total stress tensor, the solid displacement and the fluid displacement relative to the solid, respectively; |${p_f}$| denotes the fluid pressure; |$\rho = ( {1 - \phi } ){\rho _s} + \phi {\rho _f}$|⁠, with |${\rho _s}$| and |${\rho _f}$| being the density of solid and fluid, respectively; and the dynamic permeability is given by (Johnson 1987):

where |${\kappa _0}$| and |${\alpha _\infty }$| denote the static permeability and tortuosity, respectively.

The constitutive relationships are written as follows:

where modulus |$C = \alpha \bar{M}$|⁠; and I denotes the identity tensor.

The above gives equations of wave motion considering multishaped microcracks. It should be noted that those equations are consistent with the classical Biot's theory in terms of forms. The multishaped microcrack squirt flow affects the modulus |$\bar{M}$| [eqs (5) and (6)].

APPENDIX B: GENERAL SOLUTIONS FOR FRACTURE SCATTERING

The general solutions for fracture scattering in a saturated porous background with multishaped microcracks can be obtained following Song's (2017) derivation. According to Wang et al. (2024a), the zero-order components of the general solutions are given by:

where |$A( {\omega ,k} )$|⁠, |$B( {\omega ,k} )$| and |$C( {\omega ,k} )$| denote the scattered wavefield coefficients that are determined by fracture boundary conditions; J(∗) denotes the 1st class Bessel function. In addition, the background body waves amplitude ratios |${\chi _i}\,\,( {i = 1,2,3} )$| between infiltration displacement and solid displacement are given by:

and parameters |${Q_i}\,\,( {i = 1,2,3} )$|⁠:

APPENDIX C: METHOD FOR SOLVING THE NORMAL INCIDENT P-WAVE SCATTERED WAVEFIELD COEFFICIENT |${{\rm{A}}_{\rm{P}}}( {{\rm{\omega , 0}}} )$|

By substituting the general solutions into the normal incident P-wave boundary, the integrals for the scattering coefficient |${A_p}( {\omega ,k} )$| can be expressed as follows (Wang et al. 2024a):

where

with

Then convert integral eq. (C1) into Noble (1963) form:

where

Finally, the integral eq. (C5) can be converted into the form of the 2nd class Fredholm integral:

where

Eq. (C9) can be solved linearly as follows:

where |${\delta _{ij}}$| and N_dis denote the Kroencker delta function and the linear discretized points, respectively, and |${N_d} = {d {/ {\vphantom {d {( {2 \cdot {N_{\rm dis}}} )}}} } {( {2 \cdot {N_{\rm dis}}} )}}$|⁠.

After obtaining |${\Phi _p}( {\omega ,j{N_d}} )$| in eq. (C13) and using it in eq. (C12), |${U_p}( {\omega ,k} )$| in eq. (C6) can be obtained. Therefore, the coefficient |${A_p}( {\omega ,0} )$| is given as:

APPENDIX D: METHOD FOR CALCULATING THE NORMAL INCIDENT SV WAVE FAR-FIELD SCATTERING AMPLITUDE |${{\rm{f}}_{{\rm{SV}}}}( {\rm{\omega }} )$|

Similarly, to the P wave cases, by substituting the general solutions into the normal incident SV wave boundary, the integrals for the scattering coefficient |${A_{SV}}( {\omega ,k} )$| can be expressed as follows (Wang et al. 2024a):

where

and

where y is set as the default precision for numerical calculations.

Following the same steps in eqs (C9)–(C13), the coefficient |${C_{SV}}( {\omega ,k} )$| can then be obtained through the relationship with |${A_{SV}}( {\omega ,k} )$|⁠:

Therefore, the scattering amplitude in eq. (11) can be calculated by:

where |${\Phi _{sv}}( {\omega ,\tilde{\tau }} )$| is the same as calculated in eq. (C9), but for the normal incident SV wave [eq. (2)].

APPENDIX E: COMPLEX-VALUED WAVE VELOCITIES FOR FRACTURE-INDUCED EFFECTIVE VTI ROCKS IN THE Y-Z PLANE

Complex-valued wave velocities [in eqs (30) and (31)] can be obtained through eq. (28). The formulae for wave velocities in the y-z plane are calculated as follows (Mavko et al. 2020):

where

and θ denotes the wave incidence angle to the z-axis.