Fracture spacing behavior in layered rocks subjected to different driving forces: a numerical study based on fracture infilling process

Abstract

Natural layered rocks subjected to layer-parallel extension typically develop an array of opening-mode fractures with a remarkably regular spacing. This spacing often scales with layer thickness, and it decreases as extension increases until fracture saturation is reached. To increase the understanding of how these opening-mode fractures form in layered rocks, a series of 2D numerical simulations are performed to investigate the infilling process of fractures subjected to different driving forces. Numerical results illustrate that any one of the following could be considered as a driving force behind the propagation of infilling fractures: thermal effect, internal fluid pressure, direct extension loading, or pure compressive loading. Fracture spacing initially decreases with loading process, and at a certain ratio of fracture spacing to layer thickness, no new fractures nucleate (saturated). Both an increase in the opening of the infilled fractures and interface delamination are observed as mechanisms that accommodate additional strain. Interface debonding stops the transition of stress from the neighboring layers to the embedded central layer, which may preclude further infilling of new fractures. Whatever the driving force is, a large overburden stress and a large elastic contrast between the stiff and soft layers (referred to as a central or fractured layer and the top and bottom layers) are key factors favoring the development of tensile stress around the infilled fractures in the models. Fracture spacing is expected to decrease with increasing overburden stress. Numerical results highlight the fracturing process developed in heterogeneous and layered sedimentary rocks which provides supplementary information on the stress distribution and failure-induced stress redistribution., It also shows, in detail, the propagation of the fracture zone and the interaction of the fractures, which are impossible to observe in field and are difficult to consider with static stress analysis approaches.

Keywords

opening-mode fractures / fracture spacing / driving force / numerical simulation / heterogeneity

Introduction

Opening-mode fractures (joints or veins) are prevalent structural features of the Earth’s crust occur in many different rock types and tectonic environments, and are often periodically distributed with spacing linearly related to the thickness of the fractured layer. Figure 1 and Figure 2 show the natural fractured behavior at both surface and depth in the Earth’s crust. As these fractures play a major role in the subsurface flow of fluids, such as hydrocarbons and groundwater contaminants in most geologic environments, it is necessary to investigate which, and then how, different possible driving forces lead to such fracture formation.

The process of fracture formation in layered rocks was described as “sequential infilling” ( Gross, 1993). Kelly and Tyson (1965) and Price (1966) first presented the important concept of “fracture saturation,” which is closely related to fracture spacing. Fracture saturation is a phenomenon where no further infilling of fractures takes place once this critical ratio (fracture spacing to fractured layer thickness ratio) is reached. The additional strain is accommodated by the opening of existing fractures or other types of deformation ( Aveston et al., 1971; Wu and Pollard, 1995). There are numerous possible sources within the Earth’s crust, including lithostatic forces (removal or addition of material due to deposition or erosion), fluid pressure, tectonic forces, and thermal energy, in addition to other geological processes, such as magmatic intrusion or folding, which can produce stress strong enough to form a new fracture ( Ramsay and Huber, 1987; Engelder and Peacock, 2001).

The fracture formation driven by a far-field fracture-normal tensile stress is the most widely studied (e.g., Lachenbruch, 1961; Pollard and Segall, 1987; Olson, 1993; Gross et al., 1997; Tang et al., 2008). As the tension strain increases, the fracture spacing decreases approximately as the inverse of the remote strain, by fractures nucleating and propagating between earlier formed fractures. Eventually the fractures reach such a close spacing that no more fractures can infill, even with increasing strain. Instead, the existing fractures continue to open to accommodate the applied strain. In addition, the stress transition from tensile to compressive (e.g., Bai and Pollard, 2000a, b) and the slip effect between the fractured layer and adjacent layers (e.g., Ji et al., 1998; Jain et al., 2007) were considered as the reasons for precluding further infilling of new fractures. Although absolute tensile stresses can occur under certain geological circumstances (e.g., above the neutral fiber in a buckle fold), the predominant state of stress in the brittle crust is either triaxial or polyaxial compression. Such opening-mode fracturing may be induced by other driving forces, such as drying or heating effect ( Helgeson and Aydin, 1991; Fischer et al., 1995; Gross and Engelder, 1995; Engelder and Fischer 1996; Li et al., 2009). Contraction or expansion, coupled with adhesion to neighboring layers, leads to a buildup of stress in materials, and when the stress exceeds the local tensile strength, fractures nucleate in the material. In such cases, thermal mismatch stress (TMS) plays an important role in the fracturing process, especially for heterogeneous materials ( Kingery, 1955; Zuo et al., 2010). Moreover, many researchers believe that fracture propagation is driven by fluids, whether by pore pressure changes ( Secor, 1965; Engelder and Lacazette, 1990; Fischer et al., 1995; Sibson, 1996; Larsen et al., 2010; Li et al., 2012) or by sub-critical crack growth following fluid-induced rock weakening ( Savalli and Engelder, 2005). To explain such fracture formation, Ladeira and Price (1981) and Price and Cosgrove (1990) conceptually proposed that the joints in thick beds are produced by hydraulic fracturing. Fischer et al. (1995) investigated the stress distribution around the joints in layered rocks and found that the spacing of fluid-driven joints should depend on lithology and fluid pressure, at least for homogeneous fluid pressure distributions in joints that are under compression. In general, the fracture formation induced by highly pressurized fluid is divided into two types: internal hydraulic fracture and hydraulic intrusion fracture ( Mandl, 2005). In addition to thermal effect and fluid pressure, pure compression can also induce opening-mode fracturing in layered rocks. For example, Bourne (2003) showed that tensile stress could develop in multilayers under uniform remote compression provided that the contrast in the elastic properties between the layers is sufficient. De Joussineau and Petit ( 2007) investigated the development of tensile stress and fracture spacing in elastic multilayers under biaxial compressive strain conditions.

In this study, to increase the current understanding of how opening-mode fractures initiate, are developed, and to what extent fracture saturation occurs in layered rocks under different possible driving forces, a series of 2D numerical simulations, based on a uniform model, are performed. First, the typical three-layer model without pre-existing fractures is applied to investigate the progressive evolution of fractures. The typical three-layer model holds all of the material parameters constant while varying boundary conditions simulate different possible sources of driving fracture infilling. Second, the stress states between two adjacent fractures for a typical three-layer model with pre-existing fractures are numerically simulated. The term ‘‘fracture’’ in this study describes a planar discontinuity that shows predominant opening-mode displacement and cuts through the fractured layer, i.e., extends from one of the interfaces of the fractured layer to the other. The term ‘‘flaw’’ is used for a primarily opening-mode discontinuity that does not fully extend the thickness of the fractured layer. Fractures form by propagation of flaws or a coalescence of multiple flaws.

A brief introduction to the numerical code RFPA2D

In this study, all the simulations are conducted by using RFPA2D ( Tang et al., 2002; Zhu and Tang, 2004; Li et al., 2013). RFPA2D represents a two-dimensional finite element code that can simulate the fracture and the failure processes of quasi-brittle rocks. By modelling with RFPA2D, the rock model is assumed to be composed of many small-scale elements. For heterogeneous rocks, the material properties, including Young’s modulus, Poisson’s ratio, and strength properties, are randomly distributed throughout the analysis domain. The material parameters of these elements are assumed to follow a Weibull distribution as defined by the probability density function ( Weibull, 1951):

(1) \varphi ={\displaystyle \frac{m}{{\mu }_{\mathrm{0}}}}{\left({\displaystyle \frac{\mu }{{\mu }_{\mathrm{0}}}}\right)}^{m-\mathrm{1}}\exp \left[-{\left({\displaystyle \frac{\mu }{{\mu }_{\mathrm{0}}}}\right)}^m\right],

where μ is the macroscopic magnitude of parameters (such as the Young’s modulus, Poisson’s ratio, strength properties), i.e., the real values obtained from laboratory tests, while the scale parameter μ₀ is related to the average of the element parameter corresponding to μ, i.e., the input values for the numerical simulation. The homogeneity index m is a parameter defined by the shape of the distribution function that, in turn, defines the degree of material heterogeneity; a larger m implies a more homogeneous material, and vice versa. A heterogeneous material can be produced numerically in a computer simulation for a material composed of many elements. Each element is assumed to be isotropic and homogeneous. Higher homogeneity indices lead to more homogeneous numerical samples. The analysis and micro-structural observations conducted by Wong et al. (2006) require the Weibull parameter m to be greater than 2.0, but it falls into the typical range of m = 2.0-6.0 reported for engineering materials ( McClintock and Argon, 1966). Systematic studies of the homogeneity index m have been published previously (e.g., Tang et al., 2000; Zhu and Tang, 2004; Wong et al., 2006).

In RFPA2D, the finite element method is employed to obtain the stress field. Elastic damage mechanics are used to describe the constitutive law of each element. Initially, the element is considered to be elastic. The stress-strain relationship of each element is considered to be linear elastic until the given damage threshold is attained. The maximum tensile stress (strain) criterion and Mohr-Coulomb criterion are chosen as damage thresholds. The maximum tensile stress criterion is used first to determine whether an element is damaged. If the element is not damaged in tensile mode, the Mohr-Coulomb criterion is then used to determine whether the element is damaged in shear. The sign convention used throughout this paper is that compressive stresses and strains are positive.

In elastic damage mechanics, the elastic modulus may degrade gradually as damage progresses, defined as:

(2) E=\left(\mathrm{1}-D\right)E,_{\mathrm{0}}

where D represents the damage variable, and E and E₀ are the elastic moduli of the damaged and undamaged elements, respectively. The element and its damage are assumed to be isotropic, and therefore, E, E₀, and D are all scalar. D ranges from zero (0.0) for the undamaged material to one (1.0) to represent full failure.

When the element is under uniaxial tension, the constitutive relationship illustrated in Fig. 3(a) is adopted. Fig. 3(a) presents the elastic-brittle damage constitutive relation with given specific residual strength. When the tensile stress in an element reaches the tensile strength f_t0, i.e., {\sigma }_{\mathrm{3}}\le -f_{t\mathrm{0}} , the damage variable D can be calculated as:

(3)

where λ is the residual tensile strength coefficient, which is given as f_tr=λf_t0 and f_tr is called residual tensile strength. ϵ_t0 is the strain at the elastic limit, which is the so-called threshold strain for tensile damage, while ϵ_tu is the ultimate tensile strain, at which the element would be completely damaged in tension as shown in Fig. 3(a). The ultimate tensile strain is defined as ϵ_tu=ηϵ_t0 , where η is the ultimate strain coefficient. Under multiaxial stress states, the element still becomes damaged in tensile mode when the equivalent major tensile strain \overline{ϵ} increases above the threshold strain, ϵ_t0. The equivalent principal strain, \overline{ϵ} , is defined as follows:

(4) \overline{ϵ}=-\sqrt{{\langle -{ϵ}_{\mathrm{1}}\rangle }^{\mathrm{2}}+{\langle -{ϵ}_{\mathrm{2}}\rangle }^{\mathrm{2}}+{\langle -{ϵ}_{\mathrm{3}}\rangle }^{\mathrm{2}}},

where ϵ₁, ϵ₂, and ϵ₃ are the three principal strains, and〈 〉is a function defined as follows:

(5)

The constitutive law for an element subjected to multiaxial stresses can be obtained by substituting the strain ϵ in Eq. (3) with the equivalent strain \overline{ϵ} . The damage variable is then expressed as:

(6)

The degradation of an element can then be calculated by Eqs. (2) and (6). However, it must be emphasized that when D = 1, it can be calculated from Eq. (2) with the damaged elastic modulus as zero, which would make the system of equations ill-posed. Therefore, in this model, a relatively small number, i.e., 1.0×10^-5, is specified to the elastic modulus under this condition.

When the element is under uniaxial compression, the constitutive relationship illustrated in Fig. 3(b) is adopted. To describe the element damage under compressive or shear stress conditions, the Mohr-Coulomb criterion is chosen as the second damage criterion

(7) {\sigma }_{\mathrm{1}}-{\displaystyle \frac{\mathrm{1}+\sin \phi }{\mathrm{1}-\sin \phi }}{\sigma }_{\mathrm{3}}\ge f_{c\mathrm{0}},

where σ₁ is the major principal stress, σ₃ is the minor principal stress, \phi is the friction angle, and f_c0 is the uniaxial compressive strength. The damage variable under uniaxial compression is described as follows:

(8)

where ϵ_c0 is the compressive strain at the elastic limit, λ is the residual strength coefficient, and f_cr/ f_c0= f_tr/ f_t0=λ is assumed to be true when the element is under uniaxial compression or tension.

When the element is in a multi-axial stress state and its strength satisfies the Mohr-Coulomb criterion, the maximum principal strain (maximum compressive principal strain) may be evaluated at the peak value of the maximum principal stress (maximum compressive principal stress), ϵ_c0.

(9) {ϵ}_{c\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{E_{\mathrm{0}}}}\left[f_{c\mathrm{0}}+{\displaystyle \frac{\mathrm{1}+\sin \phi }{\mathrm{1}-\sin \phi }}{\sigma }_{\mathrm{3}}-\nu \left({\sigma }_{\mathrm{1}}+{\sigma }_{\mathrm{2}}\right)\right],

where ν is Poisson’s ratio.

In this respect, we assume that the shear damage evolution is only related to the maximum compressive principal strain, ϵ₁. Correspondingly, we substitute the maximum compressive principal strain ϵ₁ of the damaged element for the uniaxial compressive strain ϵ₀ in Eq. (8). Thus, Eq. (8) may be straightforwardly extended to triaxial stress states for shear damage, as follows:

(10)

From the above derivation of damage variable D, which is generally called the damage evolution law in damage mechanics, as well as Eq. (2), the elastic modulus of the damaged element at different stress or strain levels can be calculated. The unloaded element keeps its original elastic modulus and strength. That is to say, the element will unload elastically and no residual deformation will be incorporated into the element.

A single damage event represents a micro-crack (flaw) forming event that assesses the damage evolution, and the damage energy release is related to the strain energy of the element before and after its damage. Therefore, the number of damage events is counted by the number of damaged elements. The damage energy release e_f={\displaystyle \frac{{\sigma }_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{2}E}}{v}_e and the cumulative damage energy \mathrm{\Sigma }{e}_f={\displaystyle \frac{v_e}{\mathrm{2}}}\mathrm{\Sigma }{\displaystyle \frac{{\sigma }_{\mathrm{0}}^{\mathrm{2}}}{E}} are calculated from the strain energy release of damaged elements, in which σ₀ is the uniaxial compressive strength of individual element, E is the elastic modulus of individual element, and ν_e is the volume of individual element. After each damage event, we update the damage, the stress, and the strength of each element. Thus, by recording the counts of failed elements, the damage evolution associated with the progressive failure can be captured in simulation.

In RFPA2D, the model is loaded in a quasi-static fashion. At each loading step, the stress field is first examined, and those elements that are stressed beyond the pre-defined strength threshold levels are assumed to be irreversibly damaged. The stiffness and strength of the damaged elements are degraded. The model, with its associated new parameters, is then reanalyzed. The next load increment is added only when no elements strained beyond the strength-threshold remain that correspond to the equilibrium stress field and a compatible strain field. The model iterates to follow the evolution of failure along a stress path in pseudo-time.

Mathematically speaking, RFPA is using the continuum damage principle to simulate the crack propagation. The fractured rock is assumed to be an equivalent continuous medium, and the crack corresponds to the damaged zone or flaw. Although the behaviors of the mesoscopic elements are governed by a simple constitutive law, the complex macroscopic behaviors, e.g., the fracturing process, can be represented by the interactions and coalescences of the damaged elements (with more and more isolated failure elements or cracks appearing, these elements or cracks may connect to each other to form the macroscopic fractures). Neither calculation of the stress intensity factor nor mesh rezoning is required in the model. One of the main features of this model is that there is no need for a pre-existing crack to simulate the crack initiation and propagation. This approach to simulate fracture is similar to a smeared crack model, i.e., the crack is smeared over the whole element and no special singular element is used for the finite element analysis of these mesoscopic elements, which has the advantage of leaving the mesh topology untouched and greatly simplifies the simulation of crack initiation, propagation, and coalescence ( Malvar and Fourney, 1990; Pietruszczak and Xu, 1995; Pearce et al., 2000). In recent years, with the advancement in the performance of computers, an increasing number of researchers have tried to use a similar principle to model rock/concrete damage behaviors ( Pietruszczak and Xu, 1995; Fang and Harrison, 2002; Ma et al., 2011).

Numerical model set up

To numerically investigate the fracture spacing behavior in layered rocks subjected to different possible driving forces, a typical three-layer model without pre-existing fractures is employed. As shown in Fig. 4, the geometry of the model is 1,800 mm×400 mm. Thickness of the central layer (T_c) is 66.7 mm. The overall thickness of the model (T) is 400 mm, which can be defined as T=T_c+T_b+T_t. The mesh for the model contains 540×120= 64,800 elements that are produced randomly by a computer, according to the Weibull distribution for a given homogeneity index. We fix the whole bottom boundary of the model in the vertical direction and plain strain is assumed for all calculations. In the layered model, we postulate that the two materials across the layer boundaries are welded together; that is, no slip is permitted along the layer boundaries. A compressive vertical stress (σ_v) of 5.0 MPa is imposed on the top boundary.

In this study, four common possible forces that drive fracturing infilling are considered, i.e., heat effect, fluid pressure, fracture-normal extension load, and pure compression load. Consequently, four cases are modeled based on the loading configurations:

1) The thermo-mechanical loading case: the three layers are subjected to a uniform temperature arising from 15°C, at a rate of 1°C per step. The coupled process between stress/strain and heat transfer in the deforming rock mass is based on thermoelastic consolidation theory ( Li et al., 2013);

2) The hydro-mechanical loading case: a uniform water pressure with an increasing rate of 0.1 MPa/step is applied to the left and right boundaries and inner elements of the model, while the top and bottom boundaries are assumed to be impermeable. The coupled process between stress/strain and water flow in the deforming rock mass is governed by Biot’s consolidation theory ( Tang et al., 2002; Li et al., 2013);

3) The fracture-normal extension loading case: a horizontal displacement with an increasing rate of 0.05€mm/step is uniformly imposed on the left and right boundaries of the model;

4) The pure compression loading case: a horizontal stress (σ_h) of 1.0 MPa is imposed on the right and left boundaries in a vertical direction. A compressive vertical stress (σ_ν) with an increasing rate of 0.3 MPa/step is uniformly imposed on the top surface of the model.

Simulations are performed based on the uniform model shown in Fig. 4, each loading configuration may trigger a different fracture-driving mechanism. The material parameters employed in the simulation are listed in Table 1.

Numerical results

Fracture infilling process in the case under thermo-mechanical loading

Figure 5 illustrates the development and formation of fractures in detail, with the relative magnitude of the elastic modulus of the mesoscopic elements indicated by the gray scale. A brighter shade of gray corresponds to a higher Young’s modulus. In Fig. 5, the dark elements represent the nucleated flaw. Fractures form through the connection of the isolated flaws. Figure 6 shows the corresponding evolution of stress field during the fracture infilling process. The gray scale indicates the relative magnitude of the stress within the elements. The elements with a lighter shade of gray have relatively higher stresses.

According to the order of the fracture formation, we can generally divide the whole process into three stages. In the first stage (Fig. 6(a)), randomly distributed flaws form in the central layer at its weakest points where the local tensile stress induced by the TMS reaches its local tensile strength. The tensile stress is mainly caused by the heating expansion of the central layer itself in addition to the contrasts in elastic properties of the involved stratigraphy. At this stage, the initial flaws are isolated from the seed locations and do not interact. During the following stage, the process of fracture infilling takes place, where an increased number of isolated flaws connect with each other, new fractures infill continuously between the earlier formed fractures, and the spacing to fractured layer thickness ratio decreases gradually. During the last stage, the fractures are spaced so closely that no more new fractures can infill, even with increasing temperature. As shown in Fig. 6(f), the openings of the infilled fractures increase which can partially accommodate for the increasing strain induced by TMS. Furthermore, more interface fractures between adjacent layers begin to initiate and propagate. As shown in the partially enlarged boxes in Fig. 6(f), many fractures nucleate and propagate along the top and bottom interfaces. Based on the results of the numerical simulation of the model, we can assume that the number of fractures fully cut across the central layer keeps 35 in the condition of saturation.

Before the macro-fractures formed in the central layer, no localized zone or crack nucleation is observed, and thus, it is very difficult to predict where the nucleation will begin. Immediately after a certain stress is reached, zones of intense fracture event activity are identified, and nucleation sites rapidly evolve into a macro-fracture that cuts across the central layer. Figure 7 shows the temporal distribution of fracture events that occurred during the infilling process. When a new fracture forms, the fracturing process is accompanied by a rapid increase in the fracture events count. While in saturation state, the magnitude of fracture events is small and the cumulative fracture event counts nearly remain constant. This indicates that no new fractures form in the saturation state.

In addition to the sample presented in Fig. 6, four more samples (referred to as Samples 2, 3, 4, and 5) are modeled numerically while holding all of the material parameters and boundary conditions constant to eliminate the influence of the local heterogeneity on the fracture infilling. The final stage of the modelling for the additional four samples after fracture saturation is listed in Fig. 8. The overall fracturing patterns in the four models are similar to those in Fig. 6. The failure modes are sensitive to the local disorder feature of the central layer material. In the thermo-mechanical loading case, the mean critical spacing to layer thickness ratio for the five samples stated above is 0.73.

Fracture infilling process in the case under hydro-mechanical loading

As shown in Fig. 9, three stages are observed in the case under thermo-mechanical loading, including short fractures initiating, fracture infilling, and fracture saturation and interface debonding. In this case, the tensile stress is primarily from the local pore pressure and the contrasts in elastic properties of the involved stratigraphy. After flaws initiate, the process of fluid-driven fracture infilling is then governed by a balance of forces which act to open or close a flaw. Flaw opening is induced by pore pressure within a crack whereas flaw closure occurs under a compressive remote stress. Until the fluid pressure in a flaw exceeds the remote compressive stress, the crack will remain closed with no infilling fracture. The ratio of the critical fracture spacing to layer thickness for this particular case is approximately 1.174.

Fracture infilling process in the case under horizontal extension

As shown in Fig. 10, the overall fracturing pattern in this case is similar to that in Fig. 6 and Fig. 9; the profiles of the fracture are characterized by numerous branchings, T-shaped propagating (interface debonding), and termination points. In this case, the tensile stress is primarily from the local horizontal extension and the contrasts in elastic properties of the involved stratigraphy. When the fracture spacing is very small, the required tensile displacement must be considerably larger to allow for production of a new fracture because most of the additional strain is accommodated for by further opening of existing fractures and interface debonding. When the fracture spacing reaches a certain level, no opening-mode fracture can be induced, even when the tensile loading approaches the infinite. In this case, the ratio of the critical fracture spacing to layer thickness for this particular case is approximately 0.66.

Fracture infilling process in the case under pure compressive loading

The infilling process in the case under pure compressive loading is shown in Fig. 11. The overall fracturing patterns are similar to those in Figs. 6, 9, and 10. The infilling process is characterized by numerous branchings and termination points, and finally, the fractures are spaced and no new fractures can infill, even with increased vertical compressive stress. The mean ratio of the critical fracture spacing to layer thickness is approximately 1.0. Only a few interface delamination and through-going fracturing occurrences are observed in this particular case. The tensile stress is mainly due to the contrasts in elastic properties of the involved stratigraphy. Even though applying a larger horizontal compressive stress to the models may result in unfavorable conditions for the development of areas of tensile stress between the fractures, a vertical compression (representing the overburden in nature) tends to facilitate the development of tensile stress and keeps the fractures open. Of course, if the horizontal compressive strain is at least twice as large as the vertical compressive strain, the end result may be complete closure of the fractures and the tensile stress around them would no longer exist. This corresponds to extremely compressive conditions in nature ( de Joussineau and Petit, 2007).

Discussion

The influence of heterogeneity on the fracturing process

Theoretically, a new fracture should form at a certain location in between two existing adjacent fractures. However, as shown in Figs. 6, 9, 10, and 11, this is not the case for all new, i.e., the variation in fracture mode is highly sensitive to the local disordered feature of the rock. As a result, in addition to the arbitrary location of an infilling fracture, its surface is flexural and rough, which is consistent with the field observations as shown in Fig. 1. In reality, there are two types of failure for different materials: high-stress and low-strength. In a homogeneous material, failure begins at the high-stress site, whereas in heterogeneous material, such as rock, failure may start at the weaker locations due to the presence of pores, micro-fractures, grain boundaries, etc. For this reason, Fairhurst (1964) introduced the notion of “stress severity”, which represents the ratio of the theoretical stress at the moment of failure to the stress that would theoretically be necessary for failure to occur at any given point. Heterogeneity is the main explanation for the failure that occurs at locations where stress is not necessarily at its highest level.

One of the advantages of the numerical simulation is that detailed information about the stress distribution can be obtained, including data on failure-induced stress redistribution. For example, Figure 12 illustrates the tensile stress variation across section A–A (a section in the central layer as shown in Fig. 4) in the stages both before and after fracture saturation in the thermo-mechanical loading case. The influences of material heterogeneity on the stress distribution are illustrated very clearly. We find that the maximum fluctuations of σ'_min are about 50% of the mean stress for this case. It is obvious that ignoring these stress fluctuations may result in wrong conclusions if fracturing process is considered. Although the stress distribution before fracture saturation is statistically heterogeneous on a macro-scale, all elements are under tensile stress. After fracture saturation, the stress between two infilled fractures is changed from tensile to compressive stress and obvious compressive stresses are observed. When large fracture zones develop, highly non-uniform stress distributions also develop, especially when the fracture zone is not immediately stress free. Due to the local heterogeneity and the strong interaction between the isolated flaws in this high-stress field, tensile stresses can still be observed; however, most tensile stress values cannot reach the mean tension strength of the central layer (5.0 MPa). In addition, such tensile areas are not likely to allow for the initiation and propagation of subsequent fracture infilling.

The influence of delamination on the fracturing process

As shown in Figs. 5, 6, 9, and 10, a very commonly observed mode of fracture pattern formation is that fractures may nucleate and propagate along the interface between layers, that is, delamination or decohesion ( Hu et al., 1988). Physically, delamination will partially disconnect the fractured layer from the neighboring layers. In other words, once delamination begins, the deformation in the fractured layer becomes less affected by the neighboring layers. In the model, we postulate that the two materials across the layer boundaries are welded together; that is, no slip is permitted along the layer boundaries (the interface between layers is assumed to be perfectly bonded and there is no interface element). Thus, it should be emphasized that in reality, the delamination elements are actually fractured elements which are located in the central layer and adjacent to the top and bottom layers, as shown in the partially enlarged boxes in Fig. 6. Based on the infilling process described above, it has been shown that the top and bottom layers apply a local traction at all points to the central embedded layer in all models. Although the delamination can be commonly observed in the cases under thermo-mechanical, hydro-mechanical, and direct tension load, few occurrences of interface delamination and through-going fracturing are observed in the case under pure compressive load (Fig. 11). It can thus be assumed that the tensile stress resulting from the contrasts in elastic properties of the involved stratigraphy in this particular model is not sufficient for the generation of closely spaced fractures and those additional driving forces (i.e., heat effect, fluid pressure, and additional extension) are needed.

By incorporating these additional driving forces, more evident delamination can be observed. Consequently, this effect will result in a longer length scale of fracture spacing and will make the critical spacing to layer thickness ratio much greater, as discussed by other researchers (e.g., He and Hutchinson, 1989; Thouless, 1989; Cherepanov, 1994; Bai et al., 2000). As the slip effect between layers (e.g., Ji et al., 1998; Jain et al., 2007) is not considered in our modeling, the delamination can be considered as one of the main reasons for significantly changing the local stress state and inhibiting additional layer-confined fracturing. The interface delamination and through-going fracturing (fracture nucleation and propagation along the interface between layers) reduce the stress concentration induced by the contrasts in elastic properties between layers. That is, the delamination will stop the transfer of stress from the neighboring layers to the central layer and, in turn, prevent further infilling of new fractures between existing fractures in the central layer.

The influence of overburden stress on the fracturing process

In Figs. 6, 9, and 10, some crack-like flaws (partially infilling fractures) between existing infilled fractures can still be observed, even in the stage of saturation. How might such flaws be propagated to form a complete infilling fracture? Although Bai et al. (2000) did not conduct a detailed analysis of fracture controlled by the fluid pressure and thermal expansion, they concluded that the sequential infilling by such crack-like flaws between existing fractures, if such flaws exist, could occur at depth in the Earth’s crust under the action of overburden stress. To investigate the influence of overburden stress on the fracturing process of partially infilled cracks, while also keeping material parameters and all other boundary conditions identical, the models with different values of overburden stress (σ_ν= 1.0, 3.0, 7.0, and 9.0 MPa) are considered.

Figure 13 shows the numerically obtained critical ratios of fracture spacing to layer thickness under different overburden loading conditions. As long as the additional driving forces, e.g., heat effect, fluid pressure, and additional extension, etc., are incorporated, the critical ratio of fracture spacing to layer thickness in a case under a particular overburden stress can be reduced below 1.0 (the most often cited value ( Bai and Pollard, 2000a)). All the simulations produce a similar relationship between spacing and thickness: spacing decreases rapidly at first, before slowing down and approaching a nearly constant value. One explanation for this decrease is that with increased compressive overburden stress, the tensile strength of the rock between the existing fractures is enhanced accordingly. Another explanation is that the interfacial debonding becomes more difficult when the interfacial normal stress increases due to an increased compressive overburden stress.

Numerical investigation on the stress field

Regardless of the type of driving force, the fracture infilling is determined by the stress evolution between two infilled fractures. To clearly understand the evolution of stress, a modified model based on Fig. 4 is employed as shown in Fig. 14. In this model, a homogeneous medium is used (i.e., the homogeneity index m is set large enough to represent a very homogeneous material). Four fractures are pre-assigned in the central layer and are equally-spaced along the central layer, perpendicular to the long axis and fully transect the layer height (T_c). With this defined model, the transition of the stress distribution is examined as a function of the fracture spacing to the fractured layer thickness ratio (S/T_c).

Figure 15 is the numerically obtained stress transition (σ_x) along the section B–B. It is shown that, for the cases under thermo-mechanical, hydro-mechanical, directly horizontal extension, and pure compression loading, the critical spacing to fractured layer thickness ratio is approximately 0.73, 1.17, 0.66, and 1.0, respectively. This critical value determines the stress state between the adjacent fractures: when the ratio is below this critical value, the stress along line B–B is compressive, while above the critical value the stress is tensile. The compressive stress results in unfavorable conditions for the development of areas of tensile stress between the infilled fractures.

The variation of fracture spacing with different driving mechanisms

Based on the previous studies, the fracture spacing to layer thickness ratios (S/T_c) in the cases under different driving forces are summarized and listed in Table 2, as are the numerically obtained ratios of S/T_c. Field observations actually provide a broad range of fracture spacing to layer thickness ratios from less than 0.1 to greater than 10, with the most commonly reported values between 0.3 and 1.2 ( Bai and Pollard, 2000a). The range of critical values obtained in this study encompasses the often-cited spacing to layer thickness ratios in the literature for well-developed fracture sets.

Numerical results show that a far-field fracture-normal extension stress, the most widely studied condition, can generate a lower ratio of S/T_c. Although absolute tensile stresses can occur under certain geological circumstances (e.g., above the neutral fiber in a buckle fold, and a long-term movement of remote boundaries), the predominant state of stress in the brittle crust is triaxial or polyaxial compression. Even so, Crosby (1882) argued that tensile stresses did not exist in deeply buried rocks because heat and the enormous pressure of the overlying strata caused lateral expansion and prevented the contraction necessary for crack propagation. A layer can acquire a tensile stress under a remote compressive stress if it is bonded to another significantly softer and thicker layer. The amount of stress transfer between layers depends on contrasts in layer thickness and elastic properties. Tensile stress and failure within the stiff layer were demonstrated to occur for properties well within ranges typical for sedimentary rocks ( Bourne, 2003). The formation of fractures in layered rocks in compressive environments is possible, but the presence of closely spaced open fractures at depth in the Earth’s crust cannot be perfectly addressed by the pure compressive loading condition. As a consequence, the coupled thermo-mechanical or hydro-mechanical loading may be more germane to in situ fracture propagation than the pure fracture-normal loading (pure extension or compression), particularly for the sedimentary rocks buried at depth. As shown in Fig. 13, given the heat effect or fluid pressure is incorporated, closely spaced, open fractures can generate (i.e., the critical ratio of S/T_c under a particular overburden stress can be reduced below 1.0). Due to a complex tectonic and burial history with the associated thermo-poro-elastic effects, the rock material is far from being intact at depth ( Narr, 1991). Concerning fracture aperture, pressurized fluids and/or residual stresses due to contrasted properties of the layers may help to maintain the open fractures at depth.

Concluding remarks

Based on a 2D numerical study, we examined the role of different physical processes on the development of tensile stress and fracture infilling behavior in heterogeneous and layered rocks. It is shown that each of the following: thermal effect, internal fluid pressure, fracture-normal extension, and even the pure compressive loading, can be considered as possible driving forces behind the propagation of infilling fractures. Nevertheless, coupled thermo-mechanical or hydro-mechanical loading may be more germane to in situ fracture propagation than pure fracture-normal loading, particularly for the sedimentary rocks buried at depth. The ratio of fracture spacing to layer thickness decreases with depth. By incorporating heat effect or fluid pressure, the critical ratio can be reduced to a smaller value far lower than 1.0.

Whatever the driving force, the fracture infilling is determined by the stress evolution between two infilled fractures. An important concept that we have numerically confirmed and quantified during these simulations is fracture saturation, which is attributed to stress relaxation primarily caused by: 1) the mechanism of strain accommodation: the fracture spacing initially decreases as extensional strain increases. At a certain ratio of fracture spacing to layer thickness, no new fractures form and additional strain is accommodated by further opening of existing fractures; and 2) additional fractures, such as interface delamination which reduces the stress concentration induced by the contrasts in elastic properties between layers, thus preventing further infilling of new fractures.

References

[1] Aveston J, Cooper G A, Kelly A (1971). The Properties of Fiber Composites. In: Conference Proceeding of National Physical Laboratory. Guildford: IPC Science and Technology Press, 15–26

[2] Bai T, Pollard D D (2000a). Fracture spacing in layered rocks: a new explanation based on the stress transition. J Struct Geol, 22(1): 43–57

[3] Bai T, Pollard D D (2000b). Closely spaced fractures in layered rocks: initiation mechanism and propagation kinematics. J Struct Geol, 22(10): 1409–1425

[4] Bai T, Pollard D D, Gao H (2000). Explanation for fracture spacing in layered materials. Nature, 403(6771): 753–756

[5] Becker A, Gross M R (1996). Mechanism for joint saturation in mechanically layered rocks: an example from southern Israel. Tectonophysics, 257(2–4): 223–237

[6] Bourne S J (2003). Contrast of elastic properties between rock layers as a mechanism for the initiation and orientation of tensile failure under uniform remote compression. J Geophys Res, 108(B8): 2395

[7] Cherepanov G P (1994). On the theory of thermal stresses in a thin film on a ceramic substrate. J Appl Phys, 75(2): 844–849

[8] Cooke M L, Simo J A, Underwood C A, Rijken P (2006). Mechanical stratigraphic controls on fracture patterns within carbonates and implications for groundwater flow. Sediment Geol, 184(3–4): 225–239

[9] Crosby W O (1882). On the classification and origin of joint structures. Proc Boston Soc Nat Hist, 22: 72–85

[10] De Joussineau G, Petit J P (2007). Can tensile stress develop in fractured multilayers under compressive strain conditions. Tectonophysics, 432(1–4): 51–62

[11] Engelder T, Fischer M P (1996). Loading configurations and driving mechanisms for joints based on the Griffith energy-balance concept. Tectonophysics, 256(1–4): 253–277

[12] Engelder T, Lacazette A (1990). Natural hydraulic fracturing. In: Proceeding of international symposium on rock joints. Norway: Loen, 35–43

[13] Engelder T, Peacock D C P (2001). Joint development normal to regional compression during flexural-flow folding: the Lilstock buttress anticline, Somerset, England. J Struct Geol, 23(2–3): 259–277

[14] Eyal Y, Gross M R, Engelder T, Becker A (2001). Joint development during fluctuation of the regional stress field in southern Israel. J Struct Geol, 23(2–3): 279–296

[15] Fairhurst C (1964). On the validity of the Brazilian test for brittle materials. Int J Rock Mech Min Sci, 1(4): 535–546

[16] Fang Z, Harrison J P (2002). Development of a local degradation approach to the modelling of brittle fracture in heterogeneous rocks. Int J Rock Mech Min Sci, 39(4): 443–457

[17] Fischer M P, Gross M R, Engelder T, Greenfield R J (1995). Finite-element analysis of the stress distribution around a pressurized crack in a layered elastic medium: implications for the spacing of fluid-driven joints in bedded sedimentary rocks. Tectonophysics, 247(1–4): 49–64

[18] Gross M R (1993). The origin and spacing of cross joints: examples from Monterey Formation, Santa Barbara coastline, California. J Struct Geol, 15(6): 737–751

[19] Gross M R, Bahat D, Becker A (1997). Relations between jointing and faulting based on fracture-spacing ratios and fault-slip profiles: a new method to estimate strain in layered rock. Geology, 25(10): 887–890

[20] Gross M R, Engelder T (1995). Strain accommodated by brittle failure in adjacent units of the Monterey Formation, U.S.A.: scale effects and evidence for uniform displacement boundary conditions. J Struct Geol, 17(9): 1303–1318

[21] He M Y, Hutchinson J W (1989). Kinking of a crack out of an interface. J Appl Mech, 56(2): 270–278

[22] Helgeson D E, Aydin A (1991). Characteristics of joint propagation across layer interfaces in sedimentary rocks. J Struct Geol, 13(8): 897–911

[23] Hu M S, Thouless M D, Evans A G (1988). Decohesion of thin films from brittle substrates. Acta Metall, 36(5): 1301–1307

[24] Iyer K, Podladchikov Y Y (2009). Transformation-induced jointing as a gauge for interfacial slip and rock strength. Earth Planet Sci Lett, 280(1–4): 159–166

[25] Jain A, Guzina B B, Voller V R (2007). Effects of overburden on joint spacing in layered rocks. J Struct Geol, 29(2): 288–297

[26] Ji S, Saruwatari K (1998). A revised model for the relationship between joint spacing and layer thickness. J Struct Geol, 20(11): 1495–1508

[27] Ji S, Zhu Z, Wang Z (1998). Relationship between joint spacing and bed thickness in sedimentary rocks: effects of interbed slip. Geol Mag, 135(5): 637–655

[28] Kelly A, Tyson W R (1965). Tensile properties of fiber-reinforced metals: copper/tungsten and copper/molybdenum. J Mech Phys Solids, 13(6): 329–350

[29] Kingery W D (1955). Factors affecting thermal stress resistance of ceramic materials. J Am Ceram Soc, 38(1): 3–15

[30] Lachenbruch A H (1961). Depth and spacing of tension cracks. J Geophys Res, 66(12): 4273–4292

[31] Ladeira F L, Price N J (1981). Relationship between fracture spacing and bed thickness. J Struct Geol, 3(2): 179–183

[32] Larsen B, Grunnaleite I, Gudmundsson A (2010). How fracture systems affect permeability development in shallow-water carbonate rocks: an example from the Gargano Peninsula, Italy. J Struct Geol, 32(9): 1212–1230

[33] Li L C, Tang C A, Fu Y F (2009). Influence of heterogeneity on fracture behavior in multi-layered materials subjected to thermo-mechanical loading. Comput Mater Sci, 46(3): 667–671

[34] Li L C, Tang C A, Wang S Y (2012). A numerical investigation of fracture infilling and spacing in layered rocks subjected to hydro-mechanical loading. Rock Mech Rock Eng, 45: 753–765

[35] Li L C, Tang C A, Wang S Y, Yu J (2013). A coupled thermo-hydrologic-mechanical damage model and associated application in a stability analysis on a rock pillar. Tunn Undergr Space Technol, 34: 38–53

[36] Ma G W, Wang X J, Ren F (2011). Numerical simulation of compressive failure of heterogeneous rock-like materials using SPH method. Int J Rock Mech Min Sci, 48(3): 353–363

[37] Malvar L J, Fourney M E (1990). A three dimensional application of the smeared crack approach. Eng Fract Mech, 35(1–3): 251–260

[38] Mandl G (2005). Rock Joints. The Mechanical Genesis. Heidelberg: Springer, 27–48

[39] McClintock F A, Argon A S (1966). Mechanical Behavior of Materials. Reading: Addison-Wesley, 1–770

[40] Narr W (1991). Fracture density in the deep subsurface: techniques with applications to Point Arguello oil field. Am Assoc Pet Geol Bull, 75: 1300–1323

[41] Olson J E (1993). Joint pattern development: effects of subcritical crack growth and mechanical crack interaction. J Geophys Res, 98(B7): 12251–12265

[42] Pascal C, Angelier J, Cacas M C, Hancock P (1997). Distribution of joints: probabilistic modelling and case study near Cardiff (Wales, U.K.). J Struct Geol, 19(10): 1273–1284

[43] Pearce C J, Thavalingam A, Liao Z, Bicanic N (2000). Computational aspects of the discontinuous deformation analysis framework for modeling concrete fracture. Eng Fract Mech, 65(2–3): 283–298

[44] Pietruszczak S, Xu G (1995). Brittle response of concrete as a localization problem. Int J Solids Struct, 32(11): 1517–1533

[45] Pollard D D, Segall P (1987). Theoretical displacements and stresses near fracture in rock: with applications to faults, joints, veins, dikes, and solution surfaces. In: Atkinson BK ed. Fracture Mechanics of Rock. London: Academic Press, 277–349

[46] Price N J (1966). Fault and Joint Development in Brittle and Semi-brittle Rocks. Oxford: Pergamon Press, 15–31

[47] Price N J, Cosgrove J W (1990). Analysis of geologic structures. Cambridge: Cambridge University Press, 120–145

[48] Ramsay J G, Huber M I (1987). The Techniques of Modern Structural Geology (Vol 2): Folds and Fractures. San Diego: Academic Press, 1–700

[49] Savalli L, Engelder T (2005). Mechanisms controlling rupture shape during subcritical growth of joints in layered rocks. Geol Soc Am Bull, 117(3): 436–449

[50] Secor D T (1965). Role of fluid pressure in jointing. Am J Sci, 263(8): 633–646

[51] Sibson R H (1996). Structural permeability of fluid-driven fault-fracture meshes. J Struct Geol, 18(8): 1031–1042

[52] Tang C A, Liang Z Z, Zhang Y B, Chang X, Xu T, Wang D G, Zhang J X, Liu J S, Zhu W C, Elsworth D (2008). Fracture spacing in layered materials: a new explanation based on two-dimensional failure process modeling. Am J Sci, 308(1): 49–72

[53] Tang C A, Liu H, Lee P K K, Tsui Y, Tham L G (2000). Numerical studies of the influence of microstructure on rock failure in uniaxial compression — Part I: effect of heterogeneity. Int J Rock Mech Min Sci, 37(4): 555–569

[54] Tang C A, Tham L G, Lee P K K, Yang T H, Li L C (2002). Coupling analysis of flow, stress and damage (FSD) in rock failure. Int J Rock Mech Min Sci, 39(4): 477–489

[55] Taylor D W (1981). Carbonate petrology and depositional environments of the limestone member of the Carmel Formation, near Carmel Junction, Kane County, Utah. Brigham Young University Geology Studies, 28: 117–133

[56] Thouless M D (1989). Some mechanics for the adhesion of thin films. Thin Solid Films, 181(1–2): 397–406

[57] Weibull W (1951). A statistical distribution function of wide applicability. J Appl Mech, 18: 293–297

[58] Wong T F, Wong R H C, Chau K T, Tang C A (2006). Microcrack statistics, Weibull distribution and micromechanical modeling of compressive failure in rock. Mech Mater, 38(7): 664–681

[59] Wu H, Pollard D D (1995). An experimental study of the relationship between joint spacing and layer thickness. J Struct Geol, 17: 887–905

[60] Yin H M (2010). Fracture saturation and critical thickness in layered materials. Int J Solids Struct, 47(7–8): 1007–1015

[61] Zhu W C, Tang C A (2004). Micromechanical model for simulating the fracture process of rock. Rock Mech Rock Eng, 37(1): 25–56

[62] Zuo J P, Xie H P, Zhou H W, Peng S P (2010). SEM in-situ investigation on thermal cracking behavior of Pingdingshan sandstone at elevated temperatures. Geophys J Int, 181: 593–603

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg