Abstract
Elastic interactions between fractures will greatly affect the effective elasticity, which, in turn, will reshape the effective fracture parameters. The disturbance will be more complex in the fault zone due to the complicated fracture distributions. This problem is addressed by the comparison of two types of solutions: one containing the stress interaction while the other one excluding the stress interaction. The gap between the two solutions allows the quantitative estimation of stress interactions on elasticity. Furthermore, based on the orthorhombic assumption for fracture clustering in the damage zone, the effect of stress interaction on the equivalent fracture parameter is estimated. We first characterize the fracture parameters in the fault damage zone considering more realistic distributions of fractures. Then, a series of numerical simulations are conducted to study the effective parameters of the fractured model. Finally, assuming the orthorhombic system of the fracture clustering, we invert the crack density and validate the accuracy of the inversion through the incidence angle seismic velocities. Our numerical results suggest that the size of fractures will determine the dominant type of stress interactions, and thus significantly reshape the effective properties of the models regardless of the spatial distribution of the fracture. Furthermore, the stress interactions tend to underestimate the fracture density for models containing long fractures but generate a relatively satisfactory inverted fracture density for short fractures.
Introduction
In general, local stress distribution generated by a single crack hardly influences their neighbors for a sparse concentration for cracks, which, however, would be significant as the crack density exceeds the dilute limit (Zhao et al., ; Cao et al., ). That is, for a small fracture density, the resultant effective compliance tensor depends linearly on the crack density but non-linearly for a high concentration for cracks, implying a non-negligible effect of stress interactions.
Actually, there are two types of stress interactions, namely shielding and amplification, with opposite signs (Cao et al., ). Stress shielding considerably stiffens rocks while stress amplification appreciably reduces the effective elasticity (Zhao et al., ). Furthermore, the shielding and amplification effect dominate for stacked cracks and co-planar cracks, respectively (Grechka and Kachanov, ). The effect of stress interactions has been described by both numerical simulations and analytical expressions of the effective elasticity for the fractured media in a series of publications (Kachanov, ; Hopkins, ; Lapin et al., ; Cao et al., ).
Many effective analytical medium theories have been developed for the fractured media, including the self-consistent theory (SC), differential effective medium (DEM) theory, and T-matrix solution (Jakobsen, ), and a comprehensive review of these theories can be found in Mavko et al. (). These methods could be used to characterize the mechanical properties of an effective sample considering the stress interactions. Besides the analytical solutions, a series of numerical solutions have also been developed for characterizing the fractured media (Masson and Pride, ). Using the finite-element (FE) solution, Wenzlau et al. () obtain five independent elastic moduli in a heterogeneous layered medium. Similarly, Quintal et al. () propose a computationally efficient method to solve Biot's quasi-static equations of consolidation. Furthermore, based on the solution, Quintal et al. () also obtain S-wave attenuation for a medium containing periodically distributed circular inclusions. Due to the tectonic movement, the underground medium often behaves like a triclinic medium with the highest degree of anisotropy, which, however, is hard to describe using the aforementioned numerical methods. Therefore, a novelty numerical method, developed by Rubino et al. (), is used to solve the problem by generating complex-valued, frequency-dependent equivalent stiffness tensors, the advantage of which lies in the quantification of six elastic tensors for the two-dimensional (2D) case. However, neither the analytical nor the numerical solutions could explicitly identify the dominant type of stress interactions. Cao et al. () proposed a workflow to address the problem through the gaps between the analytical results without stress interactions and the numerical one with stress interactions, making the quantification of stress interaction possible.
The non-interaction approximation (NIA) or the linear slip (LS) theory treats the compliance tensors of the fractured rock as a sum of the compliance tensors of the host matrix plus individual fractures, which works well for a low concentration for cracks, and thus provides a benchmark for a quantitative estimation of stress interactions. However, when applied to fractured solids, NIA has a larger-than-expected range of applicability in the cases of dilute limit as well as the high fracture densities (Grechka and Kachanov, ). The reason lies in that, the opposite effects of shielding and amplification largely cancel one another in the medium with the fractures distributed randomly (Lapin et al., ), thus making the result similar to that of the NIA solution.
Strictly speaking, only the equivalence between the normal fracture compliance and the shear fracture compliance could ensure the orthorhombic symmetry of the fractured medium (Schoenberg and Sayers, ). However, Kachanov () believes that, in the NIA approximation, effective elasticity for the cracked model could also be orthorhombic. Based on this, Grechka () extends the case to the transversely isotropic with a host with a vertical symmetry axis (VTI), which also satisfies the aforementioned orthorhombic assumption. Therefore, it is possible to retrieve some key parameters related to the fractured information based on seismic recordings such as fracture density (Barbosa et al., ).
In the damage zone, fractures often exhibit spatial variation in the petro-physical properties, such as characteristics of the fracture length distribution (Lei and Gao, ) and the spatial distribution for the fractures inside (Harris et al., ). In recent years, there have been some attempts at describing the fracture length distribution based on the power law. Hunziker et al. () adopted the power law to describe the fracture length distribution for the systematic exploration of the attenuation sensitivity to the properties characterizing the fracture network. Lei and Gao () also characterize the stress variability in geological media based on the assumed synthetic fracture networks following power law length scaling. On the other hand, the fracture also exhibits in some forms of spatial distribution (Odling et al., ). The aspect ratio (Ξ), proposed by Jakobsen et al. (), is used to characterize the fracture spatial distribution, which, however, is more complex in nature.
For example, the fracture clustering in the damage zone, often distributed in a complex pattern and, hence, characterizing realistic spatial distributions of fractures quantitatively, is necessary for an accurate interpretation of seismic measurements. Savage and Brodsky () investigate the development of fracture spatial distributions as a function of displacement to determine whether damage around small and large faults is governed by the same process. Harris et al. () also capture the characteristics of the damage zone using distance from the fault surface. Although the fracture clustering in the damage zone is a rather common scenario in nature, description of the effects of stress interaction on the elasticity and effective fracture parameters of the damage zone remains largely unexplored, partly due to the limitations of the analytical models as well as the high computational cost because of the complex fracture distribution.
In this study, the authors address this problem by the following procedures. First, we develop a numerical model for the fracture clustering in the 2D fault damage zone. This model includes the statistical parameters of the fault damage zones observed in outcrop, followed by the investigation of effective elasticity as well as anisotropic parameters of the models. Finally, assuming the damage zone possesses the orthorhombic symmetry in the NIA, the fracture parameters for two sets of orthogonal fractures are inverted, based on which we could obtain the incidence angle dependency of the seismic wave.
Generation of Fracture Clustering in the Fault Damage Zone
For modeling the fracture clustering in the fault damage zones, the characteristics of the fracture system must be quantified in advance, including fracture length distribution, spatial distribution, and orientation distributions. Therefore, in the following section, detailed information about these characteristics, and how they are systematically incorporated into numerical models, is introduced.
Rules for the Fracture Length
The power law for fracture length distribution (Bonnet et al., ) is often used to characterize the natural fracture system. The stochastic distribution m for the fracture length l (Hunziker et al., ) is expressed as:
where L is the scale of the modeling domain while Imin and Imax are the minimum and maximum fracture lengths, respectively. Furthermore, the maximum fracture length (Imax) is only one-third of the square length (L/3 = 0.33 m), because having fractures greater than half the sample size would mechanically weaken the sample. In this study, the size L of the sample was fixed at 1.0 m. Throughout the study, the aperture of the fracture is set to be 3 mm. The exponent a is the power law length exponent, which affects the relative probability of long and short fractures. A smaller a prefers the long fractures at the expense of the short fractures. Here, we would like to set two values for the length exponent, a = 1.5 and 3, which covers two scenarios, one dominated by long fractures and the other by short fractures, respectively.
The parameter dc is the fracture density term, defined as the number of fracture centers per unit area. However, throughout the study, the fracture density, e, (Equation 1) is used as defined by Guo et al. (), which is quite different from that defined by de Dreuzy et al. (). As for dc, we set it as a function of the fracture number as follows:
where n is the fracture number inside the 2D sample, with four different values (50, 100, 150, and 200) in our study. The constant 0.03 is an optimized parameter. For keeping the units unchanged, the factor L2 is imposed into equation 1 as the denominator.
Rules for Fracture Spatial Distribution
In addition to the aforementioned length distribution, particular focus is also given to characteristics of fracture spatial distribution, one of the most challenging tasks that remain to be unexplored.
For fracture clustering, the fracture density often decays away from the main fault. Moreover, for any fracture, except the fracture with the smallest length, there is a chance that a smaller fracture will be clustered about it. That is, each fracture is placed depending on the location of the long fracture.
To construct such a kind of fractured model, an ellipsoidal volume with the aspect ratio Ξ is generated. This ellipsoid is supposed to include all centers of the clustering fractures (Figure 1A). Then, the volume is divided into a series of concentric ellipsoidal shells (Figure 1B), which cuts the vertical principal axis into equal parts. The fractures are sorted by size scale and further divided equally into concentric shells. The centers of the largest fracture are located inside the smallest ellipsoidal shell. The successive clusters of fractures with decreasing scales are then located in successive shells, which are of increasing volume (Harris et al., ).
Figure 1
Based on the extended ellipse around the major fault, a random proportion of the area (black area in Figure 2) swept by a radius of this ellipsoidal volume is used to define an azimuthal angle θ and the radius OB (Figure 2). Compared with the adoption of the random angle, the use of the swept area (black zone in Figure 1) would produce a bias toward the major axis of the ellipse, which is more representative of the natural fault zone.
Figure 2
For the location of the fracture center C, the authors would like to introduce the normalized function, t(r), as the function of normalized distance, r, along with the radius OB, with r = 1 being the shell boundary and r = 0 being the shell center. The normalized function has the form given in equation 3:
where the value p = 0.15, which is used throughout this study, resembles a pyramid-like displacement profile, making the fractures distributed more evenly over the major fault.
For a better description of the generation of the numerical model, a workflow chart to describe the generation procedure is presented in Figure 3.
Figure 3
Rules for Fracture Orientation
According to Harris et al. (), fractures inside both the Moab and Ninety Fathom fault damage zones are generally oriented sub-parallel to the major fault. Therefore, in the following models presented here, the strike of the fractures is set according to the Gaussian distribution law, each with a standard deviation of 10°.
A smaller Ξ results in a locally dense fracture clustering, thus leading to more fracture intersections, which, however, would bias the realistic fracture densities. Therefore, in this study, none of the fractures inside the models have intersections with each other.
Numerical Methodology
Methodology of the Numerical Simulation
Due to tectonic movement, the rock mass often shows a certain degree of deformation, resulting in the triclinic crystal models with the highest degree of anisotropy, which could characterize a homogeneous distribution of cracks.
In order to study the elastic properties of this complex medium, three numerical simulations with various boundary conditions are conducted here (Rubino et al., ) to obtain the elasticity matrix of the model. Detailed information about the process is given in Table 1 and as follows: (1) a fixed displacement is imposed on the upper boundary while keeping the vertical displacement to the other boundaries zero (Figure 4A); (2) the displacement is applied on the lateral side while keeping the displacement vertical to the other boundaries zero (Figure 4B); (3) a simple shear test is conducted (Figure 4C). The displacements Δu are all static rather than oscillatory.
Table 1
The boundary conditions for the three models (Figure 4).
Where L, R, B, and T are the left, right, bottom, and top boundaries of the domain, respectively.
Figure 4
In this work, for each test (Figure 4) with specified boundary conditions, the stress–strain relation (Equation D-2) is solved based on the finite-element method (FEM). Detailed information on the process can be found in Appendix D.
A cost function based on the three sets of stress–strain relations could be used to determine the equivalent stiffness tensors (six unknowns) by using the least square method. Details describing the process can be found in the publication reported by Rubino et al. ().
Two-Dimensional NIA Solution
According to the NIA theory, the effective elastic compliance tensor S of a cracked medium is determined by two additive components: the compliance tensor of the background medium (Sb) and that of the cracks (ΔS).
where the crack compliance tensor ΔS represents the accumulative contributions of the fractures to the effective compliance tensors, which is a function of the Eshelby () tensors:
where φ represents the crack porosity, J is the fourth-order symmetric identity tensor, Si are the compliance tensors of the individual fractures, and Sbcorresponds to the tensors of host matrix surrounding the fractures, respectively; the components SEshelby represents the Eshelby tensor (Eshelby, ; Guo et al., ).
For the 2D case, the ellipsoidal fracture inside the volume could be treated as an infinite cylinder (a3 → ∞); therefore, we have the Eshelby tensors in Appendix B (Masson and Pride, ). Therefore, the input parameters include the tensors of the host matrix Sb, semiaxis of the fractures, fracture orientations, Poisson ratio, and so on.
In order to get the compliant tensors, a series of parameters are input in advance, including the host matrix Sb; all fracture apertures are assumed to be 3 mm; and the fracture length and the fracture number for two orthogonal fractures are selected by using ant colony algorithms.
The NIA, considering no stress interactions, works well-under the dilute assumption of crack densities. Moreover, according to equation (6), there is no parameter describing the fracture spatial distribution or fracture size distribution; therefore, it could not provide detailed information about fracture clustering.
Evaluation of the Stress Interactions
According to Hudson's theory (Hudson et al., ), the effective stiffness tensors of the medium are composed of three parts:
where the summation of the first and second terms corresponds to the inverse matrix of SNIA (equation 6), while the third term ΔΔC is the result of the stress interaction between different fractures. C equals the CNUM (Equation D-2), which could be obtained through the numerical simulation directly. Therefore, a comparison between the and CNUM, allows for the quantification of the effect of stress interactions.
Numerical Results
In this part, a series of fractured models are introduced, with various spatial and length distributions. The corresponding elasticity and anisotropic parameters, affected by the stress interactions, are then studied using the FEM, respectively. Finally, assuming the orthotropy of the fracture clustering, we get the two sets of inverted fracture parameters3.
Fracture Parameter Setup
For a better illustration of the considered fracture networks, we vary the fracture density (e), fracture size distribution (a in equation 1), and aspect ratio (Ξ in Figure 1A) of the fracture clustering. For each parameter combination, 20 stochastic fracture networks are generated. Three typical examples are presented in Figure 5.
Figure 5
For each of the fracture densities (Table 2), a series of numerical models with a given fracture density is presented, with the fracture centers located in different areas for each model (Figure 5).
Table 2
| Case 1–2 | Case 3–4 | Case 5–6 | Case 7–8 | |
|---|---|---|---|---|
| a = 1.5 | 0.1056 | 0.1806 | 0.2938 | 0.3815 |
| a = 3.0 | 0.0083 | 0.0144 | 0.0258 | 0.0341 |
Fracture densities for eight cases with various values of a.
We choose the hypothetical dry sample with isotropic background medium, which has the following parameters: λ = 7.83 GPa; μ = 19.74 GPa (Guo et al., ). Moreover, the filling material is the gas whose Young's modulus is 0.066 GPa and the Poisson ratio υ = 0. Additionally, the definition of the fracture density in the 2D sample is introduced as follows:
where n defines the total fracture number inside the 2D sample, r is the major radius of the elliptical 2D fracture, and A is the area of the sample (Guo et al., ). Detailed information about the fracture densities for different cases is listed in Table 2.
Effective Elasticity Affected by Stress Interactions
To gain insight into the dominant stress interactions in the fracture clustering, the elastic NIA modulus (Equations B-1, B-2), which neglects stress interactions, is used to compare with the numerical results considering the stress interaction, the gap between them allows for a quantitative description of the dominant stress interactions.
For C33, opposite signs of the gaps between the NIA and the numerical results can be obtained for both long fractures (a = 1.5 in Figure 6A) and short fractures (a = 3.0 in Figure 6B), implying different dominant stress interactions determined by the fracture size. For instance, a smaller C33 produced by NIA compared with the numerical one implies a dominant amplification effect, as given in Figure 6A, while a greater C33 obtained from NIA suggests a weakly dominant shielding effect for case b (Figure 6B; Cao et al., ). Similar rules could also be applied to the shear modulus C44 in Figure 7.
Figure 6
Figure 7
Moreover, modulus discrepancies between different models with various aspect ratios (Ξ) of the bounding ellipsoid are quite different for C33. In Figure 6A, the discrepancies for C33 with different Ξ are negligible for long fractures but more significant for short fractures (a = 3.0 in Figure 6B). Conversely, C44 with various Ξ is greater for long fractures but almost negligible for short fractures.
This can be explained by the effect of the fracture interactions. For long fractures, it is believed that a smaller aspect ratio (Ξ) leads to a smaller distance between the fracture surfaces, which thus leads to a greater shielding effect, especially for the models containing long fractures. However, for the fracture clustering in this study, most of the long fractures concentrate at the core part, and therefore, Ξ variation contributes little to the distance between the close fracture, and thus, the additional shielding effect due to the increment in Ξ is comparably negligible. For short fracture, due to the cancellation of stress interactions (Kachanov, ), the net effect of stress interaction is small, and therefore, the additional stress interaction caused by a small variation in Ξ will lead to a more significant change in C33. On the other side, for the shear test (Figure 4C), the shear stress distribution is quite different from the compressive stress in the compressive test (Figures 4A,B; Cao et al., ); it is more than a function of the distance between the ajacent fracture surface; further exploration is expected in the future.
Anisotropic Properties
Besides the effective stiffness tensors (C44 and C33), the anisotropic properties of the fractured samples (Figure 5) could be described by ε and δ. For orthorhombic media, the velocity anisotropic parameters in the 2D can be computed according to equations C-1, C-2.
For the models with long fractures (a = 1.5 in Figure 8A), the shielding effect corresponds to a larger C33. Meanwhile, since the fractures are almost parallel with the X-axis, variation in C22 due to the stress interaction is negligible. Therefore, the stress interaction minimizes the contrast between C22 and C33, corresponding to a decreasing ε, according to equation C-1.
Figure 8
However, for ε in the models with short fractures (a = 3 in Figure 8B), similar discrepancies could be observed, but in a reverse pattern, that is, the numerical result is greater than the NIA result. This is understandable since the dominant amplification effect leads to a smaller C33 (Figure 6B).
According to the definition of δ by Woodruff et al. (), δ means the difference between the short-offset velocity and the vertical velocity. For the model containing long fractures, the vertical velocity, mainly determined by the shielding effect, leads to a hardening effect that compensates for the shrinking elasticity due to the individual fracture, corresponding to a smaller gap between the short-offset velocity and vertical velocity, which, in turn, leads to a smaller δ in Figure 9A. Furthermore, a greater Ξ means a less dense fracture clustering, and thus a weaker shielding effect, corresponding to a greater δ in Figure 9B.
Figure 9
In contrast to the model with long fractures, δ for the model containing short fractures is greater compared with the NIA solution (Figure 9B), due to the dominant amplification effect.
Inversion for Fracture Parameters
Actually, a model containing multiple sets of vertical dry fractures in an isotropic background matrix behaves more like orthotropic media, suggesting that the accumulative contribution of variously oriented fractures to the effective elastic properties is equivalent to that of only two principal fracture sets (Lapin et al., ). However, previous publications focused on the application of the approximation in the randomly located models where the net stress interaction is negligible (Grechka and Kachanov, , Grechka and Kachanov, ; Lapin et al., ). As for the damage zone with significant stress interaction, the effects of two types of stress interactions on the effective fracture densities for the two orthogonal sets have not been fully explored yet. We first evaluate the incidence angle dependence of the seismic velocity, so as to check the accuracy of the inverted result. Then, the fracture densities, ex and ey, associated with these two principle fracture sets are inverted from the effective elasticity.
Inversion Methodology
The problem of characterizing multiple fractures has been studied extensively by Grechka and Tsvankin () and Grechka (). Using 2D NIA, the authors of the previous publications observed that the accumulative effect of variously oriented fractures on the effective elasticity approaches that of just two orthogonal sets of fractures. Thus, we prefer to employ the NIA solution (Appendix B) to invert the equivalent orthogonal fracture parameters for two orthogonal fractures.
where the unknown vector for fracture characterization contains four parameters:
where n is the fracture number while r is the radius of the fracture; the subscripts x and y refer to the principal axes of a 2D Cartesian coordinate system. In accordance with equation 8, ex and ey for the two orthogonal fracture sets can be predicted based on the assumed orthorhombic effective elasticity.
Based on the stiffness tensors through the numerical solution (Figures 6, 7), we adopt the 2D NIA to compute for a trial vector . Using the ant colony algorithm, we could get the expected fracture parameters, which could minimize the discrepancies (Equation 11) between the trial result and the numerical result. Based on the unknowns in the vector ., the fracture density could be deduced according to equation 6. Detailed information about the inversion procedure is displayed in Figure 10.
Figure 10
Besides, in NIA, the fracture densities for two sets of orthorhombic fractures could also be obtained directly like the eigenvalues of the second-order fracture–density tensor estimated in equation 9.
The difference between the two types of fracture densities lies in that, both densities fit well with each other under the dilute fracture density assumption. However, at high fracture densities, the two densities will be quite different due to the stress interactions, which also allow for a quantitative description of the stress interaction on the inverted fracture densities.
Incidence Angle Dependence of the Inverted Seismic Velocity
So far, the inversion methodology based on the orthorhombic assumption for fracture clustering has been displayed. In the following, by employing the inverted fracture parameters (Equation 8), the comparability of the inverted result through the incidence angle dependence of the seismic wave is to be checked.
The inverted P- and SV-wave velocities, as functions of the inverted parameters for the two equivalent sets of fractures (blue solid/dashed lines in Figure 11), fit well with the numerical seismic velocities based on the real fracture parameters (Equation A1–A4), suggesting that the inverted parameters obtained are reasonable and that, as expected, the accuracy of the orthorhombic approximation is acceptable.
Figure 11
However, some minor gaps in the magnitude of the SV-wave anisotropy exist, especially for the case with long fractures (Figure 11A). These are expected for the following reasons: for long fractures (a = 1.5, Figures 11A–C), the medium is the one with the triclinic system, and thus, the nonzero C24 and C34 (Equation 5), due to the stress interactions, lead to greater SV-wave fluctuations.
Effect of Stress Interactions on the Inverted Fracture Densities
Considering the good agreement between the numerical velocities and the inverted data (Figure 11), the authors would like to further explore the effect of stress interactions on the inverted fracture parameters.
For long fractures (a = 1.5 in Figure 12A), the inverted ey is gradually biased toward the smaller values, especially at the high fracture densities, where the relative error reaches nearly 50%. A similar bias that was observed by Grechka and Kachanov () corresponds to the hardening effect of the shielding effect. On the other side, for short fractures (a = 3 in Figure 12B), the greater slope of the inverted fracture density ey (solid lines) suggests a slightly dominant amplification effect of fracture interactions, whose relative error is about 10%, which is much smaller compared with that of the long fracture.
Figure 12
Moreover, for both long (Figures 12A–C) and short fractures (Figures 12D–F), similar gaps between the blue solid and blue dashed lines for the models with various Ξ suggest that aspect ratio (Ξ) contributes little to inverted fracture densities.
On the other hand, for the wave propagates along with the fracture surface, the effect of stress interaction is negligible; therefore, the inverted fracture density ex has an overall agreement with the NIA ones generally.
Actually, the inverted fracture densities based on NIA share the same size, which, however, bias the real models with various length scales and spatial distributions. Furthermore, due to the limitation of the NIA solution, the inverted result could not provide information about the spatial distribution of fractures.
Discussions
According to many previous publications about the T-matrix method (Jakobsen, , Jakobsen and Chapman, ), the aspect ratio (Ξ) characterizing the spatial distribution of fractures (Jakobsen et al., , Zhao et al., ) represents the conditional probability of finding another inclusion given the position of a known inclusion. However, the definition of such a parameter limits its application in describing fracture clustering. In the present work, the aspect ratio is used to characterize the boundary of the fracture clustering. Such a definition change for the aspect ratio (Ξ) allows for the description of the properties for the fracture clustering.
According to Lapin et al. (), for the NIA, the effective elastic properties of an isotropic matrix with fractures distributed in certain spatial patterns always possess the orthorhombic symmetry, regardless of the orientations of the fractures without the same symmetry. However, due to the effect of stress interactions, the model behaves like a triclinic medium rather than an orthorhombic one. Thus, the nonzero moduli C24 and C34, due to the stress interaction, contribute more to the SV-wave anisotropy, leading to more obvious SV-wave velocity fluctuation (such as the SV-wave velocity in Figure 11A).
Both the fracture size and the aspect ratio of the fracture cluster contribute to the effective elasticity. However, an agreement between C33 for the models with various aspect ratios suggests that the contribution of aspect ratio to the stress interaction is negligible, partly due to the fact that variation in the aspect ratio influences the distance between the short fractures distributed at the outer part of the damage zone more when compared to the long fractures at the core part.
It should be noted that all the fractures inside the cluster have no intersection with each other, which, however, is counterintuitive to reality. Nonetheless, according to the conclusion of Grechka and Kachanov (), the fracture intersections impose negligible influence on the effective tensors. Therefore, the conclusion presented in this study still holds acceptable accuracy for the fractured clustering in the field.
Garboczi and Berryman () believe that the error of the numerical simulation may come from statistical variation due to the inhomogeneous distribution of fractures. For the inhomogeneous distribution with a given concentration of cracks, there are multiple spatial arrangements for the cracks, which may have somewhat different elastic moduli. The error size due to inhomogeneous distribution could be assessed by computing the elastic parameters for various realizations of the same system. Therefore, all numerical results as well as the inverted parameters are determined based on the 20 different realizations by changing the crack locations with the same crack density, which could guarantee the accuracy of the conclusions in this study.
Throughout the article, the fracture aperture is set to be 3 mm. However, the aperture varies a lot in the field, which would also contribute to the effect of stress interactions. According to current knowledge, the aperture will also affect both types of stress interactions; further development in this aspect is expected in the future.
Conclusion
Aiming at the damage zone in the fault, numerical models based on the characteristics of the fracture clustering were built, based on which the authors analyzed the effective elastic tensors of the fracturing clusters and estimated the fracture parameters based on the orthorhombic assumption.
Long parallel fractures tend to generate a strong shielding effect, which, therefore, contribute more to the C33 and C44, and a smaller ε. Conversely, the weakly dominant amplification effect induced by the short fractures leads to smaller C33 and C44, and a greater ε. Moreover, the effect of the boundary aspect ratio (Ξ) on the stress interaction for fracture clustering is almost negligible for long fractures (a = 1.5) but relatively more significant for the short fractures (a = 3), where a greater Ξ corresponds to a greater amplification effect.
The inverted fracture density is also strongly influenced by the stress interaction. Amplification tends to weakly overestimate the fracture density while the shielding effect underestimates the fracture density. Furthermore, the incidence angle dependency for the P-wave velocity in fracture clustering is similar to that for the medium containing two orthogonal or principal fracture sets; however, the similarity is somewhat less satisfactory for the SV-wave velocity.
Statements
Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Author contributions
L-YF and CC conceived the research. CC wrote the manuscript and prepares the figures. L-YF reviewed and supervised the manuscript. B-YF and QG are involved in the modeling and inversion programme, respectively. All authors finally approve the manuscript and thus agree to be accountable for this work.
Funding
The Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDA14010303).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2021.643372/full#supplementary-material
References
1
BarbosaN. D.RubinoJ. G.CaspariE.HolligerK. (2018). Impact of fracture clustering on the seismic signatures of porous rocks containing aligned fracturesFracture clustering effects on WIFF. Geophysics83:A65–A68. 10.1190/geo2017-0799.1
2
BonnetE.BourO.OdlingN. E.DavyP.MainI.CowieP.et al. (2001). Scaling of fracture systems in geological data. Rev. Geophys. 39, 347–383. 10.1029/1999RG000074
3
CaoC.ChenF.FuL.-Y.BaJ.HanT. (2020a). Effect of stress interactions on anisotropic P-SV-wave dispersion and attenuation for closely spaced cracks in saturated porous media. Geophys. Prospect. 68, 2536–2556. 10.1111/1365-2478.13007
4
CaoC.FuL.-Y.BaJ.ZhangY. (2019). Frequency- and incident-angle-dependent P-wave properties influenced by dynamic stress interactions in fractured porous media. Geophysics. 84, 1–12. 10.1111/1365-2478.12446
5
CaoC.FuL.-Y.FuB.-Y. (2020b). An elastic numerical method based on the 3D complex medium. Chin. J. Geophys.63, 2836–2845. 10.6038/cjg2020N0035
6
de DreuzyJ.-R.DavyP.BourO. (2001). Hydraulic properties of two-dimensional random fracture networks following a power law length distribution: 1. Effective connectivity. Water Res. Res.37, 2065–2078. 10.1029/2001WR900011
7
EshelbyJ. D. (1957). The Determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc. R. Soc. Lond.241, 376–396. 10.1098/rspa.1957.0133
8
GarbocziE. J.BerrymanJ. G. (2001). Elastic moduli of a material containing composite inclusions: effective medium theory and finite element computations. Mech. Mater.33, 455–470. 10.1016/S0167-6636(01)00067-9
9
GrechkaV. (2007). Multiple cracks in VTI rocks: effective properties and fracture characterization. Geophysics72, D81–D91. 10.1190/1.2751500
10
GrechkaV.KachanovM. (2006a). Effective elasticity of fractured rocks: a snapshot of the work in progress. Geophysics71, W45–W58. 10.1190/1.2360212
11
GrechkaV.KachanovM. (2006b). Effective elasticity of rocks with closely spaced and intersecting cracks. Geophysics71, D85–D91. 10.1190/1.2197489
12
GrechkaV.TsvankinI. (2003). Feasibility of seismic characterization of multiple fracture sets. Geophysics68, 1399–1407. 10.1190/1.1598133
13
GuoJ.HanT.FuL.XuD.FangX. (2019). Effective elastic properties of rocks with transversely isotropic background permeated by aligned penny-shaped cracks. J. Geophys. Res. Solid Earth124, 400–424. 10.1029/2018JB016412
14
GuoJ.RubinoJ. G.BarbosaN. D.GlubokovskikhS.GurevichB. (2018). Seismic dispersion and attenuation in saturated porous rocks with aligned fractures of finite thickness: theory and numerical simulations — part 1: P-wave perpendicular to the fracture plane. Geophysics83, WA49–WA62. 10.1190/geo2017-0065.1
15
HarrisS.McAllisterE.KnipeR.OdlingN. E. (2003). Predicting the three-dimensional population characteristics of fault zones: a study using stochastic models. J. Struct. Geol.25, 1281–1299. 10.1016/S0191-8141(02)00158-X
16
HopkinsD. L. (2000). The implications of joint deformation in analyzing the properties and behavior of fractured rock masses, underground excavations, and faults. Int. J. Rock Mech. Min. Sci.37, 175–202. 10.1016/S1365-1609(99)00100-8
17
HudsonJ. A.LiuE.CrampinS. (1996). The mechanical properties of materials with interconnected cracks and pores. Geophys. J. Int.124, 105–112. 10.1111/j.1365-246X.1996.tb06355.x
18
HunzikerJ.FavinoM.CaspariE.QuintalB.RubinoJ. G.KrauseR.et al. (2018). Seismic attenuation and stiffness modulus dispersion in porous rocks containing stochastic fracture networks. J. Geophys. Res. Solid Earth123, 125–143. 10.1002/2017JB014566
19
JakobsenM. (2012). T-matrix approach to seismic forward modelling in the acoustic approximation. Stud. Geophys. Et Geodaet.56, 1–20. 10.1007/s11200-010-9081-2
20
JakobsenM.ChapmanM. (2009). Unified theory of global flow and squirt flow in cracked porous media. Geophysics74, WA65–WA76. 10.1190/1.3078404
21
JakobsenM.HudsonJ. A.JohansenT. A. (2003). T-matrix approach to shale acoustics Geophys. J. Int.154, 533–558. 10.1046/j.1365-246X.2003.01977.x
22
KachanovM. (1993). Elastic solids with many cracks and related problems. Adv. Appl. Mech.30, 259–445. 10.1016/S0065-2156(08)70176-5
23
LapinR.KuzkinV.KachanovM. (2018). On the anisotropy of cracked solids. Int. J. Eng. Sci.124, 16–23. 10.1016/j.ijengsci.2017.11.023
24
LeiQ.GaoK. (2018). Correlation between fracture network properties and stress variability in geological media. Geophys. Res. Lett. 45, 3994–4006. 10.1002/2018GL077548
25
MassonY. J.PrideS. R. (2007). Poroelastic finite difference modeling of seismic attenuation and dispersion due to mesoscopic-scale heterogeneity. J. Geophys. Res. Atmosph.112, 1642–1642. 10.1029/2006JB004592
26
MassonY. J.PrideS. R. (2014). On the correlation between material structure and seismic attenuation anisotropy in porous media. J. Geophys. Res. Solid Earth119, 2848–2870. 10.1002/2013JB010798
27
MavkoG.MukerjiT.DvorkinJ. (2009). The Rock Physics Handbook, 2nd Edn. Cambridge: Cambridge University Press.
28
OdlingN.HarrisS.VasziA.KnipeR. (2005). Properties of fault damage zones in siliclastic rocks: a modelling approach. Geolog. Soc. Lond.249, 43–59. 10.1144/GSL.SP.2005.249.01.04
29
QuintalB.SteebH.FrehnerM.SchmalholzS. (2011). Quasi-static finite element modeling of seismic attenuation and dispersion due to wave-induced fluid flow in poroelastic media. J. Geophys. Res. 116:B01201. 10.1029/2010JB007475
30
QuintalB.SteebH.FrehnerM.SchmalholzS. M.SaengerE. H. (2012). Pore fluid effects on S-wave attenuation caused by wave-induced fluid flow. Geophysics77, L13–L23. 10.1190/geo2011-0233.1
31
RubinoJ. G.CaspariE.MüllerT. M.MilaniM.BarbosaN. D.HolligerK. (2016). Numerical upscaling in 2-D heterogeneous poroelastic rocks: Anisotropic attenuation and dispersion of seismic waves. J. Geophys. Res. Solid Earth121, 6698–6721. 10.1002/2016JB013165
32
SavageH.BrodskyE. (2011). Collateral damage: Evolution with displacement of fracture distribution and secondary fault strands in fault damage zones. J. Geophys. Res. 116:B03405. 10.1029/2010JB007665
33
SchoenbergM.SayersC. (1995). Seismic anisotropy of fracture rock. Geophysics60, 204–211. 10.1190/1.1443748
34
WenzlauF.AltmannJ. B.MüllerT. M. (2010). Anisotropic dispersion and attenuation due to wave-induced fluid flow: Quasi-static finite element modeling in poroelastic solids. J. Geophys. Res. 115:B07204. 10.1029/2009JB006644
35
WoodruffW. F.RevilA.PrasadM.Torres-VerdÃnC. (2015). Measurements of elastic and electrical properties of an unconventional organic shale under differential loading. Geophysics80, D363–D383. 10.1190/geo2014-0535.1
36
ZhaoL.YaoQ.HanD. H.YanF.NasserM. (2015). Characterizing the effect of elastic interactions on the effective elastic properties of porous, cracked rocks. Geophys. Prospect.64, 157–169. 10.1111/1365-2478.12243
Summary
Keywords
stress interaction, fracture clustering, damage zone, inversion, modeling
Citation
Cao C, Fu L-Y, Fu B-Y and Guo Q (2021) Effect of Stress Interactions on Effective Elasticity and Fracture Parameters in the Damage Zones. Front. Earth Sci. 9:643372. doi: 10.3389/feart.2021.643372
Received
18 December 2020
Accepted
21 April 2021
Published
13 May 2021
Volume
9 - 2021
Edited by
Pascal Audet, University of Ottawa, Canada
Reviewed by
Ke Gao, Southern University of Science and Technology, China; Junxin Guo, Southern University of Science and Technology, China
Updates
Copyright
© 2021 Cao, Fu, Fu and Guo.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Li-Yun Fu lfu@upc.edu.cn
This article was submitted to Solid Earth Geophysics, a section of the journal Frontiers in Earth Science
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.