# Hydraulic and Mechanical Relationship of Individual Fracture in Rock under Compression and Shearing: Theoretical Study

Authors: Zhiqiang Zhou; Yu Zhao; Chaolin Wang
Journal: Geofluids
 (2020)

---

## Abstract

In this paper, a new approach has been developed for predicting the hydraulic and mechanical relationship of individual fractures subjected to normal stress and compression-shear stress. Considering that the closure process of rough fracture subjected to normal stress can be divided into two phases (linear behavior and nonlinear behavior), a relationship between normal stress and fracture aperture is derived through the minimum potential energy principle. Then, a formulation for calculating fracture permeability during shearing and compression processes is developed. Furthermore, a formulation for determining fracture aperture during the crack growth process is obtained, which is further implanted into the permeability model to predict the hydraulic behavior of fractured rock during fracture propagation. This new model not only considers the normal deformation of the fracture but also, and more importantly, integrates the effect of fracture propagation and shear dilation. Theoretical studies demonstrate that fracture permeability increases nonlinearly during fracture propagation. At last, experimental results and analytic results are compared to demonstrate the usefulness of the proposed models, and satisfactory agreements are obtained.


---

## Body

## 1. Introduction

The coupling between hydraulic and mechanical behavior is of considerable interest in several areas of rock engineering, such as underground tunnels [1], naturally fractured petroleum reservoirs [2], rock slope stability, and underground disposal of radioactive waste. Fractures of rock masses are key factors governing the coupling between hydraulic and mechanical behavior [3]. For engineering applications, an improved understanding of the stress-flow coupling of individual fracture, which is fundamental to the hydromechanical analysis of fractured rock, is needed. However, to accurately describe the hydraulic conductivity of fractured rock masses under complex stress is an issue that is not completely resolved. This difficulty mainly comes from that the heterogeneous nature of rock materials contains various natural fractures of different scales [4–6]. Consequently, many investigators have been focusing on the hydraulic and mechanical properties of fractures by conducting experimental tests and analytical research [7–9].Previous studies of fluid flow through fractured rock have mainly focused on the effects of hydromechanical coupling, fracture characteristics (such as aperture distribution, surface roughness, and contact area), or a combination of both effects [10–12]. It is examined by a multitude of experiments that aperture distribution, surface roughness, and contact area have a significant influence on flow behaviors in fractures [13–16], and all of them are established with a function of normal or shear stress. Laboratory investigations on fractured rocks show that shear deformation can induce fracture opening, which is recognized as shear dilatancy [17]. After recognizing the dilation phenomenon in shear processes, shear dilatancy and related hydromechanical behavior of fractured rocks were studied extensively. Yeo et al. [18] investigated the effect of shear displacement on the aperture and permeability of rock fracture by using radial and unidirectional flow test and numerical simulation. Chen et al. [19] presented a strain-dependent formulation to address the change in hydraulic conductivity resulted from normal and shear loads. Wei and Liu [20] presented a fractal-based model for fracture deformation under shearing and compression. Although considerable progress has been made in the past four decades, models for describing the hydromechanic behavior with different conditions should be further studied [21–23].The existing fractures in rock mass may close, open, or grow, and new cracks may be induced under external loads [8, 19, 24–27]. During the process of fracture propagation, the length and aperture of fracture, which is widely regarded as factors influencing the flow properties of the rock mass, will change. In all the above models for describing fracture hydromechanical properties, however, only the obvious factors such as aperture distribution, surface geometry, and contact area are well investigated; the variation of fracture length induced by fracture propagation is not taken into account. In practice, the effect of fracture propagation on fluid flow was not ignored, especially when the rock mass suffers relatively large external loads. Specifically, a variation of fracture length will change the average fracture aperture and flow path. However, to the best knowledge of the authors, limited literature for investigating fracture hydromechanical properties has considered the effect of propagation on fracture characteristics.Therefore, in this paper, we extend the above work and develop a new approach to establish the hydromechanical behavior for fractured rock under shear and normal stress. A simple review of the continuously yielding fracture model proposed by Cundall and Hart [28] is first introduced, and then, an improved relationship for shear stress and displacement is built with the consideration of mechanical properties on fracture surface roughness. Based on the principle of minimum potential energy, the relationships between stress and deformation of open fracture are developed for the compression process and shearing-compression process, respectively. And further considering fracture propagation under compression-shear, a formulation for calculating fracture aperture in the crack growth region is obtained. Unlike most previous models, the new model gives particular emphasis to the effect of fracture propagation and shear dilation.


## 2. Materials and Methods

The fluid flow through fracture of rock masses has commonly been described by the cubic law. It is derived from the assumption that the incompressible and laminar fluid flow occurs between two smooth-parallel plates. It is often denoted as [29].
(1)Q=Cb3∇H,where Q is the flowrate, b is the fracture aperture, ∇H is the water head difference through fracture surface, and C is the constant related to gravitational acceleration, flow viscosity coefficient, and the shape and size of the fracture surface. Natural rock fractures, however, are featured with rough surfaces and variable aperture, which conflicts with the assumption. Many scholars have been focused on the studying of flow characteristics of rough fractures and proposed various modified cubic law. Neuzil and Tracy [29] put forward a probability density function to modify fracture width based on the assumption that width changes in the vertical direction and keeps a constant in a parallel direction. Walsh [15] built a modifier formula containing a ratio of area contact by imitating the heat conduction theory. To account for the deviation of the nature fractures from the ideally smooth plates, Chen et al. [19] presented an equivalent formulation for fracture permeability by adopting a dimensionless constant, ξ; it can be expressed by
(2)k=ξb212,where k is the fracture permeability. The aperture at any fracture deformation is given as [19]:
(3)b=b0+Δb,where b0 and Δb are the initial aperture and the change in aperture due to external loads.


## 3. Relationship between Shear Stress and Shear Displacement

Previous studies demonstrate that rock joints show high nonlinear shear properties under direct shear stress. It is reflected in the features that normal stiffness varies with normal stress, and internal friction angle decreases with increasing shear displacement. Cundall and Hart [28] proposed a continuously yielding joint model (referred to as the CY model). This model was intended to produce nonlinear behavior and shear dilation and could be applied to analyze the internal mechanism of progressive damage of rock joints under shear stress. Based on the CY model, a new relationship between shear stress and shear displacement is built in this paper.A demonstration of the CY model with a typical shear stress-displacement curve, normal displacement curve, and the bounding shear strength curve is shown in Figure1. The response to shear stress in the CY model is calculated incrementally as [28]
(4)Δτ=F0ksΔδs,where Δδs is the shear displacement increment, ks is shear stiffness, and F0 is the elastic recoverable ratio, which depends on the distance from the actual shear stress curve to bounding strength curve, τm (Figure 1). Theoretically, F0 is written as [31]
(5)F0=1−τ/τm1−r.Figure 1
A typical shear stress-displacement curve, normal displacement curve [30], and bounding shear strength curve [31].The factor,r, is used to reproduce the speed in approaching the bounding stress curve: a bigger r indicates a lager speed to the bounding stress. In practice, r is limited to 0.75 in order to avoid numerical noise. To simplify the equation, r is set to be zero in this paper.To analyze the effect of roughness in fractured rock, Oh [30] employed a formulation to substitute the bounding strength in the CY model with
(6)τm=σntanφr+ψ,where σn, φr, and ψ are normal stress, residual friction angle, and dilatancy angle, respectively.On the basis of test data, Barton et al. [14] proposed an empirical expression of dilatancy angle with the joint roughness coefficient (JRC), which was defined as
(7)ψmob=JRCmoblog10JCSσn,where JCS refers to the fracture surface compressive strength, subscript mob denotes that parameters change with the shear displacement during the process of shearing. Tse and Curden [32] proposed an empirical formulation to calculate JRC, and it would get the maximum value at the point of peak shear strength.Combining Eqs. (4) through (7) yields
(8)Δτ=ksΔδs1−τσntanφr+JRCmoblgJCS/σn.Combining with Eqs. (4) and (5) and using the initial condition “when δ=0, τ=0” gives
(9)τ=τm−τmexp−ksτmδs.Eq. (8) provides a new relationship between shear stress and displacement, which is highly nonlinear and is related to JRC. Previous studies focused on fracture shear properties are mainly suggested in an empirical manner and fail to explicitly describe the influence of mechanical properties on fracture surface roughness. Eq. (8) responds to our effort to consider the impact of fracture roughness during the process of shearing. Note that due to the failure of asperities, the dilation angle and JRC decrease with plastic shear displacement, which is aligned with fracture deformation.


## 4. Behavior of Fracture Deformation

In this section, we will discuss the hydromechanical behavior of a single fracture under different conditions, including normal stress, compression-shear stress, and fracture propagation process.


### 4.1. Fracture Closure under Normal Stress

As reviewed in the introduction section, for fractured rock, the relationship between normal displacement and normal stress is nonlinear. The deformation response of an open fracture subject to normal stress may include two phases: one is linear closure of the fracture, and the other is nonlinear closure due to contact areas. The observed nonlinear stress-deformation behavior of fracture could be attributed to the increasing contact areas and crushing of asperities. In this subsection, fracture closure due to the normal stress is investigated through the minimum potential energy principle.An infinite plate with rough fractured surfaces under far-field compression stress is shown in Figure2. With the working of normal stress, fracture deformation occurs, and new elastic energy is produced. Then, the energy of the fracture element is [33]
(10)U=U0+U1−W1,where U0 is the initial elastic energy stored in the fracture element. The new elastic energy, U1, induced by normal displacement, corresponds to the mode-I crack field. The work W1 applied by the normal stress σn can be expressed as [33]
(11)W1=12∑Fs=12∫s12HσnUnxds1,where F and s are the work and displacement, respectively, and Unx refers to displacement function. The region of integration, s1, depends on the fracture length 2a. The effective working force, Hσn, is composed of at least two major parts, which correspond to linear closure and nonlinear closure. For the first phase, there are no contact areas in the fracture; thus, σn is the effective stress doing work. When normal stress increases to closure stress, σc, a portion of stresses are counteracted due to surfaces contacting and asperities crushing. A reduction factor m is used to account for the effect of surfaces contacting and asperities crushing [34]. Then, Hσn can be determined by
(12)Hσn=σn,σn<σcσnm,σn≥σc,Figure 2
Sketch of mechanical analysis of fracture in compression state.Note that with normal stress increasing, the growth rate of effective working stress decreases because of the increasing contact areas and the number of contacts. Therefore, the factorm should be dependent on effective stress.The Griffith theory made it possible for the first time to like within the limits of the thermodynamic approach the process of failure of a solid with the energy transformations accompanying the formation of new fracture surfaces [35]. Griffith [36] argued that the value of deformation in elliptic fracture under compression is maximal in the middle and decreases to both ends. Therefore, the displacement function is assumed as
(13)Unx=±α1δna2−x21/2a,where a, α1, and δn are the half fracture length, correction coefficient of the fracture shape, and normal displacement, respectively. Substituting Eq. (13) into Eq. (11), and integrating, we obtain
(14)W1=π2aBα1Hσnδn,where B is the width of the material. According to Ashby and Hallam [33], new elastic energy U1 accumulated by σn is
(15)U1=α3EBδn2,where α3 and E are the constant and fracture elasticity modulus, respectively.The minimum potential energy principle for a conservative elastic system states that of all geometrically admissible displacement fields, that which brings the loaded structure into elastic equilibrium also minimizes the potential energy functional [37]. Minimizing the energy of element during the deformation process, hence
(16)dUdδn=0.Applying Eqs. (10)–(15) into Eq. (16), and ignoring the changes of U0, yields
(17)δn=Δb=πaα14α3EHσn.Hence, a theoretical model of fracture closure for rough fracture is formulated, which is totally determined by effective working stressHσn and a set of parameters. It is to be noted that for the first phase (σn<σc) Eq. (17) can be rewritten as Hook’s law by defining 8α3E/πα1 as the fracture stiffness.
(18)εn=δn2a=σnkn,where εn is the normal strain and kn is fracture stiffness.


### 4.2. Fracture Aperture under Compression-Shear Stress

A fracture subject to compression-shear stress will slip and dilate to facilitate the opening of fracture. As illustrated in Figure3(a), the normal stress and shear stress of fracture surface are computed by
(19)σn=12σ1+σ3+12σ1−σ3cos2α,τs=12σ1−σ3sin2α,where α is the inclination measured in terms of the angle from the principal axes.Figure 3
Sketch of mechanical analysis of fracture in the compression-shear state: (a) stress state; (b) deformation feature.
(a)(b)


#### 4.2.1. Fracture Aperture before Propagation

Considering a fracture to be embedded into a rock subject to compression-shear stress, a similar method and notations in the previous subsection are used to investigate the fracture deformation. Then, the energy of fracture element subject to compression-shear stress is given as(20)U=U0+U1+U2−W2,where U2 is new elastic energy induced by shear displacement and can be calculated by [33]
(21)U2=α3EBδs2,where δs refers to shear displacement. One can observe from Figure 3(b) that the working applied by the shear load is positive and applied by the normal load is negative due to shear dilation. Then, the whole acting, W2, applied by τs and σn can be expressed as
(22)W2=12∑Fs=12∫s12τsUsxds1−12∫s12HσnUnxds1,where Usx is the shear displacement function and can be similarly given as follows
(23)Usx=±α1δsa2−x21/2a.Substituting Eqs. (13) and (23) into Eq. (22) yields
(24)W2=π2aBα1τsδs−Hσnδn.The relationship between the total fracture displacement (δ) and δs or δn can be written as
(25)δn=δsinψ,δs=δcosψ,where ψ is the dilation angle. Applying Eqs. (13), (15), and (20)–(25) into Eq. (16) yields
(26)δ=πaα14α3Eτscosψ−Hσnsinψ.As discussed before, the fractured rock under shear stress will dilate to facilitate the opening of fracture, which means few stresses are counteracted during the dilation process. Therefore, the function of working stressHσn is set to σn. Combining Eqs. (25) and (26), the change of aperture subject to compression-shear stress is
(27)Δb=δn=πaα1sinψ4α3Eτscosψ−σnsinψ.Combining Eqs. (2), (3), and (27), the stress-dependent formulation of fracture permeability is determined by
(28)k=ξb0+Δb212=ξ12b0+aα1sinψα3Eτscosψ−σnsinψ2.Eq. (9) presents a relation between shear stress τs and shear displacement δs. Accordingly, a formulation of strain-dependent permeability is developed by combining Eq. (9) and (28). It is important to note that our Eq. (28) is applicable when the rock dilates. Other mechanic deformations can be similarly analyzed with this method.


#### 4.2.2. Fracture Aperture during Propagation

Eq. (27) describes fracture aperture change during compression-shear stress without considering the effect of fracture propagation. In a general case, fractures under high stress involve aperture changes in the crack growth region. As illustrated in Figure 4(a), when stress (σ1orσ3) increases to the initial cracking value, stress-induced cracks are produced. Assuming that stress-induced cracks are straight lines and parallel to σ3 (Figure 4(a)), then the normal stress of wing crack, Fn, equals to σ3. A number of experimental results and theoretical calculations demonstrate that the initial cracking is in accordance with mode-I under compression-shear [24, 38, 39].Figure 4
Sketch of the mechanical analysis of wing cracks in the compression-shear state: (a) stress state; (b) deformation feature.
(a)(b)Similar to Eq. (22), the total work W3 done by τs, σn, and Fn is
(29)W3=12∑Fs=12∫s12τsUsxds1−12∫s12σnUnxds1+12∫s22FnUFnxds2,where UFn is the displacement of the wing cracks. The integral region, s2, depends on the wing crack length l. Martin [40] provided a general equation for estimating the length of the crack extension once the cracks initiate; it was written as:
(30)l=A+B0expCσ1σ3,where A, B0, and C are the curve-fitting parameters.In order to illuminate the fracture deformation, the system (original crack and wing crack pair) is replaced by an equivalent single straight crack [33]. Similar to Eqs. (13) and (23), the following relationships are formulated by substituting a+l for a(31)Usx=±α1δsa+l2−x21/2a+l,Unx=±α1δna+l2−x21/2a+l.Assuming that the displacement of the wing crack,UFnx, is proportional to the distance from the wing crack tip, then it can be defined by
(32)UFnx=±α2δ1−xa+l,where α2 is the constant. Integrating Eq. (29) with Eqs. (31) and (32), W3 is derived to be
(33)W3=aBα1τsδs−σnδn1−11+L21/2+1+Larcsin11+L+α2α3aLBδ,where L=l/a. By employing the Taylor series expansion, Eq. (33) can be adequately approximated by
(34)W3=2aBα1τsδs−σnδn1−161+L2+α2α3aLBδ.Note that the value in the bracket of the first item approximately equals to 1 forL>1; then, Eq. (34) is reduced to the following form:
(35)W3=aB2α1τsδs−2α1σnδn+α2α3Lδ.Combining Eqs. (15), (16), (20), (21), and (35) yields
(36)δ=aα1α3Eτscosψ−σnsinψ+α22α1σ3L.Therefore, the aperture and fracture permeability during the fracture propagation process can be expressed as(37)b=b0+aα1sinψα3Eτscosψ−σnsinψ+α22α1σ3L,(38)k=ξ12b0+aα1sinψα3Eτscosψ−σnsinψ+α22α1σ3L2.Eq. (38) made it possible for the first time to predict the permeability evolution of rock during fracture propagation.


## 4.1. Fracture Closure under Normal Stress

As reviewed in the introduction section, for fractured rock, the relationship between normal displacement and normal stress is nonlinear. The deformation response of an open fracture subject to normal stress may include two phases: one is linear closure of the fracture, and the other is nonlinear closure due to contact areas. The observed nonlinear stress-deformation behavior of fracture could be attributed to the increasing contact areas and crushing of asperities. In this subsection, fracture closure due to the normal stress is investigated through the minimum potential energy principle.An infinite plate with rough fractured surfaces under far-field compression stress is shown in Figure2. With the working of normal stress, fracture deformation occurs, and new elastic energy is produced. Then, the energy of the fracture element is [33]
(10)U=U0+U1−W1,where U0 is the initial elastic energy stored in the fracture element. The new elastic energy, U1, induced by normal displacement, corresponds to the mode-I crack field. The work W1 applied by the normal stress σn can be expressed as [33]
(11)W1=12∑Fs=12∫s12HσnUnxds1,where F and s are the work and displacement, respectively, and Unx refers to displacement function. The region of integration, s1, depends on the fracture length 2a. The effective working force, Hσn, is composed of at least two major parts, which correspond to linear closure and nonlinear closure. For the first phase, there are no contact areas in the fracture; thus, σn is the effective stress doing work. When normal stress increases to closure stress, σc, a portion of stresses are counteracted due to surfaces contacting and asperities crushing. A reduction factor m is used to account for the effect of surfaces contacting and asperities crushing [34]. Then, Hσn can be determined by
(12)Hσn=σn,σn<σcσnm,σn≥σc,Figure 2
Sketch of mechanical analysis of fracture in compression state.Note that with normal stress increasing, the growth rate of effective working stress decreases because of the increasing contact areas and the number of contacts. Therefore, the factorm should be dependent on effective stress.The Griffith theory made it possible for the first time to like within the limits of the thermodynamic approach the process of failure of a solid with the energy transformations accompanying the formation of new fracture surfaces [35]. Griffith [36] argued that the value of deformation in elliptic fracture under compression is maximal in the middle and decreases to both ends. Therefore, the displacement function is assumed as
(13)Unx=±α1δna2−x21/2a,where a, α1, and δn are the half fracture length, correction coefficient of the fracture shape, and normal displacement, respectively. Substituting Eq. (13) into Eq. (11), and integrating, we obtain
(14)W1=π2aBα1Hσnδn,where B is the width of the material. According to Ashby and Hallam [33], new elastic energy U1 accumulated by σn is
(15)U1=α3EBδn2,where α3 and E are the constant and fracture elasticity modulus, respectively.The minimum potential energy principle for a conservative elastic system states that of all geometrically admissible displacement fields, that which brings the loaded structure into elastic equilibrium also minimizes the potential energy functional [37]. Minimizing the energy of element during the deformation process, hence
(16)dUdδn=0.Applying Eqs. (10)–(15) into Eq. (16), and ignoring the changes of U0, yields
(17)δn=Δb=πaα14α3EHσn.Hence, a theoretical model of fracture closure for rough fracture is formulated, which is totally determined by effective working stressHσn and a set of parameters. It is to be noted that for the first phase (σn<σc) Eq. (17) can be rewritten as Hook’s law by defining 8α3E/πα1 as the fracture stiffness.
(18)εn=δn2a=σnkn,where εn is the normal strain and kn is fracture stiffness.


## 4.2. Fracture Aperture under Compression-Shear Stress

A fracture subject to compression-shear stress will slip and dilate to facilitate the opening of fracture. As illustrated in Figure3(a), the normal stress and shear stress of fracture surface are computed by
(19)σn=12σ1+σ3+12σ1−σ3cos2α,τs=12σ1−σ3sin2α,where α is the inclination measured in terms of the angle from the principal axes.Figure 3
Sketch of mechanical analysis of fracture in the compression-shear state: (a) stress state; (b) deformation feature.
(a)(b)


### 4.2.1. Fracture Aperture before Propagation

Considering a fracture to be embedded into a rock subject to compression-shear stress, a similar method and notations in the previous subsection are used to investigate the fracture deformation. Then, the energy of fracture element subject to compression-shear stress is given as(20)U=U0+U1+U2−W2,where U2 is new elastic energy induced by shear displacement and can be calculated by [33]
(21)U2=α3EBδs2,where δs refers to shear displacement. One can observe from Figure 3(b) that the working applied by the shear load is positive and applied by the normal load is negative due to shear dilation. Then, the whole acting, W2, applied by τs and σn can be expressed as
(22)W2=12∑Fs=12∫s12τsUsxds1−12∫s12HσnUnxds1,where Usx is the shear displacement function and can be similarly given as follows
(23)Usx=±α1δsa2−x21/2a.Substituting Eqs. (13) and (23) into Eq. (22) yields
(24)W2=π2aBα1τsδs−Hσnδn.The relationship between the total fracture displacement (δ) and δs or δn can be written as
(25)δn=δsinψ,δs=δcosψ,where ψ is the dilation angle. Applying Eqs. (13), (15), and (20)–(25) into Eq. (16) yields
(26)δ=πaα14α3Eτscosψ−Hσnsinψ.As discussed before, the fractured rock under shear stress will dilate to facilitate the opening of fracture, which means few stresses are counteracted during the dilation process. Therefore, the function of working stressHσn is set to σn. Combining Eqs. (25) and (26), the change of aperture subject to compression-shear stress is
(27)Δb=δn=πaα1sinψ4α3Eτscosψ−σnsinψ.Combining Eqs. (2), (3), and (27), the stress-dependent formulation of fracture permeability is determined by
(28)k=ξb0+Δb212=ξ12b0+aα1sinψα3Eτscosψ−σnsinψ2.Eq. (9) presents a relation between shear stress τs and shear displacement δs. Accordingly, a formulation of strain-dependent permeability is developed by combining Eq. (9) and (28). It is important to note that our Eq. (28) is applicable when the rock dilates. Other mechanic deformations can be similarly analyzed with this method.


### 4.2.2. Fracture Aperture during Propagation

Eq. (27) describes fracture aperture change during compression-shear stress without considering the effect of fracture propagation. In a general case, fractures under high stress involve aperture changes in the crack growth region. As illustrated in Figure 4(a), when stress (σ1orσ3) increases to the initial cracking value, stress-induced cracks are produced. Assuming that stress-induced cracks are straight lines and parallel to σ3 (Figure 4(a)), then the normal stress of wing crack, Fn, equals to σ3. A number of experimental results and theoretical calculations demonstrate that the initial cracking is in accordance with mode-I under compression-shear [24, 38, 39].Figure 4
Sketch of the mechanical analysis of wing cracks in the compression-shear state: (a) stress state; (b) deformation feature.
(a)(b)Similar to Eq. (22), the total work W3 done by τs, σn, and Fn is
(29)W3=12∑Fs=12∫s12τsUsxds1−12∫s12σnUnxds1+12∫s22FnUFnxds2,where UFn is the displacement of the wing cracks. The integral region, s2, depends on the wing crack length l. Martin [40] provided a general equation for estimating the length of the crack extension once the cracks initiate; it was written as:
(30)l=A+B0expCσ1σ3,where A, B0, and C are the curve-fitting parameters.In order to illuminate the fracture deformation, the system (original crack and wing crack pair) is replaced by an equivalent single straight crack [33]. Similar to Eqs. (13) and (23), the following relationships are formulated by substituting a+l for a(31)Usx=±α1δsa+l2−x21/2a+l,Unx=±α1δna+l2−x21/2a+l.Assuming that the displacement of the wing crack,UFnx, is proportional to the distance from the wing crack tip, then it can be defined by
(32)UFnx=±α2δ1−xa+l,where α2 is the constant. Integrating Eq. (29) with Eqs. (31) and (32), W3 is derived to be
(33)W3=aBα1τsδs−σnδn1−11+L21/2+1+Larcsin11+L+α2α3aLBδ,where L=l/a. By employing the Taylor series expansion, Eq. (33) can be adequately approximated by
(34)W3=2aBα1τsδs−σnδn1−161+L2+α2α3aLBδ.Note that the value in the bracket of the first item approximately equals to 1 forL>1; then, Eq. (34) is reduced to the following form:
(35)W3=aB2α1τsδs−2α1σnδn+α2α3Lδ.Combining Eqs. (15), (16), (20), (21), and (35) yields
(36)δ=aα1α3Eτscosψ−σnsinψ+α22α1σ3L.Therefore, the aperture and fracture permeability during the fracture propagation process can be expressed as(37)b=b0+aα1sinψα3Eτscosψ−σnsinψ+α22α1σ3L,(38)k=ξ12b0+aα1sinψα3Eτscosψ−σnsinψ+α22α1σ3L2.Eq. (38) made it possible for the first time to predict the permeability evolution of rock during fracture propagation.


## 4.2.1. Fracture Aperture before Propagation

Considering a fracture to be embedded into a rock subject to compression-shear stress, a similar method and notations in the previous subsection are used to investigate the fracture deformation. Then, the energy of fracture element subject to compression-shear stress is given as(20)U=U0+U1+U2−W2,where U2 is new elastic energy induced by shear displacement and can be calculated by [33]
(21)U2=α3EBδs2,where δs refers to shear displacement. One can observe from Figure 3(b) that the working applied by the shear load is positive and applied by the normal load is negative due to shear dilation. Then, the whole acting, W2, applied by τs and σn can be expressed as
(22)W2=12∑Fs=12∫s12τsUsxds1−12∫s12HσnUnxds1,where Usx is the shear displacement function and can be similarly given as follows
(23)Usx=±α1δsa2−x21/2a.Substituting Eqs. (13) and (23) into Eq. (22) yields
(24)W2=π2aBα1τsδs−Hσnδn.The relationship between the total fracture displacement (δ) and δs or δn can be written as
(25)δn=δsinψ,δs=δcosψ,where ψ is the dilation angle. Applying Eqs. (13), (15), and (20)–(25) into Eq. (16) yields
(26)δ=πaα14α3Eτscosψ−Hσnsinψ.As discussed before, the fractured rock under shear stress will dilate to facilitate the opening of fracture, which means few stresses are counteracted during the dilation process. Therefore, the function of working stressHσn is set to σn. Combining Eqs. (25) and (26), the change of aperture subject to compression-shear stress is
(27)Δb=δn=πaα1sinψ4α3Eτscosψ−σnsinψ.Combining Eqs. (2), (3), and (27), the stress-dependent formulation of fracture permeability is determined by
(28)k=ξb0+Δb212=ξ12b0+aα1sinψα3Eτscosψ−σnsinψ2.Eq. (9) presents a relation between shear stress τs and shear displacement δs. Accordingly, a formulation of strain-dependent permeability is developed by combining Eq. (9) and (28). It is important to note that our Eq. (28) is applicable when the rock dilates. Other mechanic deformations can be similarly analyzed with this method.


## 4.2.2. Fracture Aperture during Propagation

Eq. (27) describes fracture aperture change during compression-shear stress without considering the effect of fracture propagation. In a general case, fractures under high stress involve aperture changes in the crack growth region. As illustrated in Figure 4(a), when stress (σ1orσ3) increases to the initial cracking value, stress-induced cracks are produced. Assuming that stress-induced cracks are straight lines and parallel to σ3 (Figure 4(a)), then the normal stress of wing crack, Fn, equals to σ3. A number of experimental results and theoretical calculations demonstrate that the initial cracking is in accordance with mode-I under compression-shear [24, 38, 39].Figure 4
Sketch of the mechanical analysis of wing cracks in the compression-shear state: (a) stress state; (b) deformation feature.
(a)(b)Similar to Eq. (22), the total work W3 done by τs, σn, and Fn is
(29)W3=12∑Fs=12∫s12τsUsxds1−12∫s12σnUnxds1+12∫s22FnUFnxds2,where UFn is the displacement of the wing cracks. The integral region, s2, depends on the wing crack length l. Martin [40] provided a general equation for estimating the length of the crack extension once the cracks initiate; it was written as:
(30)l=A+B0expCσ1σ3,where A, B0, and C are the curve-fitting parameters.In order to illuminate the fracture deformation, the system (original crack and wing crack pair) is replaced by an equivalent single straight crack [33]. Similar to Eqs. (13) and (23), the following relationships are formulated by substituting a+l for a(31)Usx=±α1δsa+l2−x21/2a+l,Unx=±α1δna+l2−x21/2a+l.Assuming that the displacement of the wing crack,UFnx, is proportional to the distance from the wing crack tip, then it can be defined by
(32)UFnx=±α2δ1−xa+l,where α2 is the constant. Integrating Eq. (29) with Eqs. (31) and (32), W3 is derived to be
(33)W3=aBα1τsδs−σnδn1−11+L21/2+1+Larcsin11+L+α2α3aLBδ,where L=l/a. By employing the Taylor series expansion, Eq. (33) can be adequately approximated by
(34)W3=2aBα1τsδs−σnδn1−161+L2+α2α3aLBδ.Note that the value in the bracket of the first item approximately equals to 1 forL>1; then, Eq. (34) is reduced to the following form:
(35)W3=aB2α1τsδs−2α1σnδn+α2α3Lδ.Combining Eqs. (15), (16), (20), (21), and (35) yields
(36)δ=aα1α3Eτscosψ−σnsinψ+α22α1σ3L.Therefore, the aperture and fracture permeability during the fracture propagation process can be expressed as(37)b=b0+aα1sinψα3Eτscosψ−σnsinψ+α22α1σ3L,(38)k=ξ12b0+aα1sinψα3Eτscosψ−σnsinψ+α22α1σ3L2.Eq. (38) made it possible for the first time to predict the permeability evolution of rock during fracture propagation.


## 5. Model Validation

In this section, we mainly focused on evaluating the usefulness of the proposed models against the experimental observations from the literature.


### 5.1. Validation of the Proposed Model under Normal Load

Many investigators have explored the relationship between fracture aperture and normal stress by conducting experiments. Two experimental datasets from Ref. [41, 42] are selected to validate the model proposed in Section 4.1.The information on the datasets adopted in this paper is shown in Table1. The values for fracture length and initial aperture are directly determined from laboratory tests. Fracture elastic modulus can be determined from the stress-strain curve in the compression test; in this paper, E is set to 4 GPa for both fractured rocks. The determination of threshold stress σc, for intact rock, may utilize the method of fracture volume strain proposed by Martin and Chandler [43]. For separated rock, σc is estimated from the closure curves by identifying the turning point from linear behavior to nonlinear behavior. Moreover, an empirical function of the factor m=σn/nlnσn is proposed based on the experimental data in Figure 5, in which n is a parameter related to contact areas and asperity degradation. The remaining parameter values (α1/α3 and n) are identified by fitting Eq. (17) with the measured data in Figure 5.Table 1
Parameter values used for matching the data obtained under normal stress.
Sample no.Initial aperture (b0)Fracture elastic modulus (E)Fracture half length (a)Threshold stress (σc)Fitting parametersRef.α1/α3nDr-10.38 mm4 GPa47 mm1.9 MPa7.730.59Malama et al. [42]Grd-10.505 mm4 GPa47 mm1.5 MPa17.620.23IG100.526 mm4 GPa75 mm0.5 MPa1.880.75Matsuki et al. [41]IG150.589 mm4 GPa75 mm0.63 MPa4.360.27IG200.39 mm4 GPa75 mm1.09 MPa7.990.16IG250.697 mm4 GPa75 mm0.4 MPa8.780.13Figure 5
Comparisons between the observed closure curves (from (a) Ref. [41] and (b) Ref. [42]) and theoretical results from Eq. (17).
(a)(b)The comparisons between the observed closure curves and theoretical results from Eq. (17) are depicted in Figure 5, which manifests that the mechanical deformation of the fractured rock under normal stress is well reproduced by the proposed model. It should be noted that for the two experimental datasets fitting, the values of the correlation coefficient (R2) are above 0.95. This confirms that the proposed model is capable of reproducing the change of fracture aperture during the normal loading process.


### 5.2. Validation of the Proposed Model under Compression-Shear Load

Esaki et al. [44] conducted a comprehensive study regarding shear deformation and dilatancy dependence of hydraulic conductivity for rock fractures through laboratory experiments. The coupled shear-flow tests were conducted on granite samples with a size of 120 mm in length, 100 mm in width, and 80 mm in height. Note that the tested sample has been split into two parts to create an artificial fracture. Eq. (27) should be applied for calculating fracture deformation rather than Eq. (37). The values of the parameters used for matching the data of Esaki et al. [44] are listed in Table 2. The initial aperture and fracture length were directly obtained from Ref. [44]. He also provided values for JRC and JRS, which are 9 and 162 MPa, respectively. The dilation angle, ψ, is identified as a constant for simplicity. In this case, its value can be determined from Eq. (7). The bounding strength can be calculated using Eq. (6), but the better alternative is to directly identify it from the tests. The rest of the parameter values are determined by fitting the experimental curves using relevant equations. It should be noted that the coupled shear-flow tests were conducted for one cycle of forwarding and reverse shearing, and only the forward experimental results are adopted in this paper. The results from the reverse shearing process can be also matched in the same method.Table 2
Parameter values used for matching the data obtained under compression-shear stress.
Normal stressInitial aperture (b0)Fracture elastic modulus (E)Dilation angle (ψ)Fracture lengthBounding strength (τm)Fitting parametersξksα1/α31 MPa0.15 mm4 GPa19.89°60 mm2.1 MPa0.040.32114.85 MPa0.15 mm4 GPa13.59°60 mm5 MPa0.061.2529.3710 MPa0.15 mm4 GPa10.89°60 mm11 MPa0.142.037.8420 MPa0.15 mm4 GPa8.18°60 mm22.1 MPa0.083.393.01The comparisons between the aperture analytically predicted by using the proposed model given in Eq. (27) and those measured from coupled shear-flow tests are shown in Figure 6. As can be seen in Figure 6, for different normal stresses, the observed shear dilatancy behavior is fairly well captured by Eq. (27) with a single group of parameter values. Slight discrepancies between the experimental curves and analytical curves are observed when the normal stress is small. This can be explained by the fact that the dilation angle undergoes negative decay under small normal stress.Figure 6
Comparisons between the observed curves (from Ref. [44]) and theoretical results from Eq. (27): (a) normal stress 1 MPa; (b) normal stress 5 MPa; (c) normal stress 10 MPa; and (d) normal stress 20 MPa.
(a)(b)(c)(d)


### 5.3. Hydraulic Behavior during Crack Propagation Process

Chen et al. [45] performed several triaxial compression tests on granite and found that the evolution of rock permeability experienced a clear decrease phase in the initial microcrack closure region and a dramatic permeability increase phase in the crack growth region. In this subsection, the observed data in the crack growth region are compared with the predicted results. However, the hydraulic behavior with multiple sets of fractures is complex and obviously distinguished from that of a single fracture due to varied aperture values and sizes in each fracture. Min et al. [46] proposed the concepts of “equivalent frequency” and “equivalent aperture” to approximate the complex fracture deformation. Using this method, the fracture permeability with multiple sets of fractures may be formulated in a simple form as follows:
(39)k=ξeqbeq212=ξeq12b0+aeqα1sinψα3Eτscosψ−σnsinψ+α22α1σ3L2,where the subscript eq refers to demote the equivalent concept, L=l/aeq, and the stress-induced crack length, l, can be obtained using Eq. (30).For simplicity, two-dimensional analysis is employed here, and the associative flow rule is used for modeling the dilatancy of fracture. Assuming that the equivalent fracture is orientated at the angle of failure surface (i.e.,α=π/4+ϕ), the stresses (τs and σn) of fracture surface can be determined by Eq. (19). Since there is no model of plasticity, the dilatancy angle is assumed equal to the friction angle. The parameter values identified in the model are summarized in Table 3, in which the initial aperture, elastic modulus, and friction angle are directly obtained from Ref. [45]. The remaining parameters are obtained by fitting the measured data in Figure 7 using Eq. (39). Note that α1, α2, and α3 in Eq. (39) are set to 1 to reduce the number of fitting parameters. The curves illustrated in Figure 7 depict that the evolution of permeability in the crack growth region is fairly well captured by the analytic model, supporting the usefulness of our approach.Table 3
Parameter values used for matching the data in the crack extension phase.
Initial apertureElastic modulus (E)Friction angle (ϕ)Fracture length (aeq)ξeqAB0C0.05μm45 GPa35°1000μm0.040.60.2340.277Figure 7
Comparisons between the observed curves (from Ref. [45]) and theoretical results from Eq. (39).


## 5.1. Validation of the Proposed Model under Normal Load

Many investigators have explored the relationship between fracture aperture and normal stress by conducting experiments. Two experimental datasets from Ref. [41, 42] are selected to validate the model proposed in Section 4.1.The information on the datasets adopted in this paper is shown in Table1. The values for fracture length and initial aperture are directly determined from laboratory tests. Fracture elastic modulus can be determined from the stress-strain curve in the compression test; in this paper, E is set to 4 GPa for both fractured rocks. The determination of threshold stress σc, for intact rock, may utilize the method of fracture volume strain proposed by Martin and Chandler [43]. For separated rock, σc is estimated from the closure curves by identifying the turning point from linear behavior to nonlinear behavior. Moreover, an empirical function of the factor m=σn/nlnσn is proposed based on the experimental data in Figure 5, in which n is a parameter related to contact areas and asperity degradation. The remaining parameter values (α1/α3 and n) are identified by fitting Eq. (17) with the measured data in Figure 5.Table 1
Parameter values used for matching the data obtained under normal stress.
Sample no.Initial aperture (b0)Fracture elastic modulus (E)Fracture half length (a)Threshold stress (σc)Fitting parametersRef.α1/α3nDr-10.38 mm4 GPa47 mm1.9 MPa7.730.59Malama et al. [42]Grd-10.505 mm4 GPa47 mm1.5 MPa17.620.23IG100.526 mm4 GPa75 mm0.5 MPa1.880.75Matsuki et al. [41]IG150.589 mm4 GPa75 mm0.63 MPa4.360.27IG200.39 mm4 GPa75 mm1.09 MPa7.990.16IG250.697 mm4 GPa75 mm0.4 MPa8.780.13Figure 5
Comparisons between the observed closure curves (from (a) Ref. [41] and (b) Ref. [42]) and theoretical results from Eq. (17).
(a)(b)The comparisons between the observed closure curves and theoretical results from Eq. (17) are depicted in Figure 5, which manifests that the mechanical deformation of the fractured rock under normal stress is well reproduced by the proposed model. It should be noted that for the two experimental datasets fitting, the values of the correlation coefficient (R2) are above 0.95. This confirms that the proposed model is capable of reproducing the change of fracture aperture during the normal loading process.


## 5.2. Validation of the Proposed Model under Compression-Shear Load

Esaki et al. [44] conducted a comprehensive study regarding shear deformation and dilatancy dependence of hydraulic conductivity for rock fractures through laboratory experiments. The coupled shear-flow tests were conducted on granite samples with a size of 120 mm in length, 100 mm in width, and 80 mm in height. Note that the tested sample has been split into two parts to create an artificial fracture. Eq. (27) should be applied for calculating fracture deformation rather than Eq. (37). The values of the parameters used for matching the data of Esaki et al. [44] are listed in Table 2. The initial aperture and fracture length were directly obtained from Ref. [44]. He also provided values for JRC and JRS, which are 9 and 162 MPa, respectively. The dilation angle, ψ, is identified as a constant for simplicity. In this case, its value can be determined from Eq. (7). The bounding strength can be calculated using Eq. (6), but the better alternative is to directly identify it from the tests. The rest of the parameter values are determined by fitting the experimental curves using relevant equations. It should be noted that the coupled shear-flow tests were conducted for one cycle of forwarding and reverse shearing, and only the forward experimental results are adopted in this paper. The results from the reverse shearing process can be also matched in the same method.Table 2
Parameter values used for matching the data obtained under compression-shear stress.
Normal stressInitial aperture (b0)Fracture elastic modulus (E)Dilation angle (ψ)Fracture lengthBounding strength (τm)Fitting parametersξksα1/α31 MPa0.15 mm4 GPa19.89°60 mm2.1 MPa0.040.32114.85 MPa0.15 mm4 GPa13.59°60 mm5 MPa0.061.2529.3710 MPa0.15 mm4 GPa10.89°60 mm11 MPa0.142.037.8420 MPa0.15 mm4 GPa8.18°60 mm22.1 MPa0.083.393.01The comparisons between the aperture analytically predicted by using the proposed model given in Eq. (27) and those measured from coupled shear-flow tests are shown in Figure 6. As can be seen in Figure 6, for different normal stresses, the observed shear dilatancy behavior is fairly well captured by Eq. (27) with a single group of parameter values. Slight discrepancies between the experimental curves and analytical curves are observed when the normal stress is small. This can be explained by the fact that the dilation angle undergoes negative decay under small normal stress.Figure 6
Comparisons between the observed curves (from Ref. [44]) and theoretical results from Eq. (27): (a) normal stress 1 MPa; (b) normal stress 5 MPa; (c) normal stress 10 MPa; and (d) normal stress 20 MPa.
(a)(b)(c)(d)


## 5.3. Hydraulic Behavior during Crack Propagation Process

Chen et al. [45] performed several triaxial compression tests on granite and found that the evolution of rock permeability experienced a clear decrease phase in the initial microcrack closure region and a dramatic permeability increase phase in the crack growth region. In this subsection, the observed data in the crack growth region are compared with the predicted results. However, the hydraulic behavior with multiple sets of fractures is complex and obviously distinguished from that of a single fracture due to varied aperture values and sizes in each fracture. Min et al. [46] proposed the concepts of “equivalent frequency” and “equivalent aperture” to approximate the complex fracture deformation. Using this method, the fracture permeability with multiple sets of fractures may be formulated in a simple form as follows:
(39)k=ξeqbeq212=ξeq12b0+aeqα1sinψα3Eτscosψ−σnsinψ+α22α1σ3L2,where the subscript eq refers to demote the equivalent concept, L=l/aeq, and the stress-induced crack length, l, can be obtained using Eq. (30).For simplicity, two-dimensional analysis is employed here, and the associative flow rule is used for modeling the dilatancy of fracture. Assuming that the equivalent fracture is orientated at the angle of failure surface (i.e.,α=π/4+ϕ), the stresses (τs and σn) of fracture surface can be determined by Eq. (19). Since there is no model of plasticity, the dilatancy angle is assumed equal to the friction angle. The parameter values identified in the model are summarized in Table 3, in which the initial aperture, elastic modulus, and friction angle are directly obtained from Ref. [45]. The remaining parameters are obtained by fitting the measured data in Figure 7 using Eq. (39). Note that α1, α2, and α3 in Eq. (39) are set to 1 to reduce the number of fitting parameters. The curves illustrated in Figure 7 depict that the evolution of permeability in the crack growth region is fairly well captured by the analytic model, supporting the usefulness of our approach.Table 3
Parameter values used for matching the data in the crack extension phase.
Initial apertureElastic modulus (E)Friction angle (ϕ)Fracture length (aeq)ξeqAB0C0.05μm45 GPa35°1000μm0.040.60.2340.277Figure 7
Comparisons between the observed curves (from Ref. [45]) and theoretical results from Eq. (39).


## 6. Conclusions

The hydromechanical behavior of fractured rock is an important issue for practical applications. This paper develops the relationships between hydraulic and mechanical behavior of individual fracture subjects to normal and shear loads. Based on this study, the following conclusions can be drawn:(1)
The formulas for describing the hydraulic and mechanical relationship in the fractured rock are derived through the minimum potential energy principle, including the cases of fractured rock under compression, under compression, and shearing. The proposed approach enables easier fracture deformation analysis to include the effect of shear dilation and is validated by favorable comparisons between measured and analytical results(2)
On the basis of the minimum potential energy principle and cubic law, a permeability model with consideration of the fracture propagation is derived. Theoretical results show that the permeability increases nonlinearly during fracture propagation, and the proposed model is capable of predicting the permeability evolution during fracture propagation



---

*Source: 8848116-2020-11-19.xml*