# The Role of In Situ Stress in Organizing Flow Pathways in Natural Fracture Networks at the Percolation Threshold

Authors: Chuanyin Jiang; Xiaoguang Wang; Zhixue Sun; Qinghua Lei
Journal: Geofluids
 (2019)

---

## Abstract

We investigated the effect of in situ stresses on fluid flow in a natural fracture network. The fracture network model is based on an actual critically connected (i.e., close to the percolation threshold) fracture pattern mapped from a field outcrop. We derive stress-dependent fracture aperture fields using a hybrid finite-discrete element method. We analyze the changes of aperture distribution and fluid flow field with variations of in situ stress orientation and magnitude. Our simulations show that an isotropic stress loading tends to reduce fracture apertures and suppress fluid flow, resulting in a decrease of equivalent permeability of the fractured rock. Anisotropic stresses may cause a significant amount of sliding of fracture walls accompanied with shear-induced dilation along some preferentially oriented fractures, resulting in enhanced flow heterogeneity and channelization. When the differential stress is further elevated, fracture propagation becomes prevailing and creates some new flow paths via linking preexisting natural fractures, which attempts to increase the bulk permeability but attenuates the flow channelization. Comparing to the shear-induced dilation effect, it appears that the propagation of new cracks leads to a more prominent permeability enhancement for the natural fracture system. The results have particularly important implications for predicting the hydraulic responses of fractured rocks to in situ stress fields and may provide useful guidance for the strategy design of geofluid production from naturally fractured reservoirs.


---

## Body

## 1. Introduction

Fractured rocks host a significant proportion of the world’s georesources, e.g., groundwater, hydrocarbon, and geothermal energy. Fractured reservoirs are known to be highly heterogeneous with the main flow paths dominated by intricate fracture networks. The spatial distribution and organization of natural fractures in the subsurface are highly complex, often exhibiting long-range correlation and spatial. Furthermore, the transmissivity of fractures shows strong variations in both magnitude and space. As a result, fractured reservoirs often accommodate strong flow heterogeneities across multiple length scales [1–3].A typical workflow for fractured reservoir characterization starts by constructing a flow model, which requires the geometrical and hydraulic properties of the studied fracture network to be defined [4]. Such a model is often difficult to build due to limited subsurface measurements. Therefore, fracture patterns from outcrop analogues are commonly used to enhance the understanding regarding the key characteristics of the fracture network [5–8]. With a thorough fracture characterization, the emergent geological model may provide a good geometrical representation of the studied fracture system; however, it may still fail completely to reproduce the system’s hydrodynamic behavior. This failure arises from the difficulties in relating surface characterization to subsurface fracture network properties (e.g., fracture aperture and connectivity), due to resolution limits of geophysical tools [9]. Moreover, since fractures are subject to in situ stresses in the subsurface environment, fracture apertures can be often strongly modified from their initial values and exhibit a highly heterogeneous distribution [10–12]. The aperture heterogeneity may lead to a highly channelized flow pattern [7, 13–15]. Extensive field observations have shown that the contribution of individual fractures in the network to the total flow is not equal [16–18]. Instead, the majority of the fluid flow is restricted to only a small number of fractures [19–22]. Consequently, in order to predict the hydrodynamic behavior of fractured rocks, it is crucial to include the impact of in situ stresses on fluid flow.A few recent studies have shown that thein situ stresses may exert a complex impact on the bulk hydraulic properties of fractured rocks due to a variety of fracture responses such as closure, sliding, dilatancy, and propagation [23–29]. Based on orthogonally imposed stress boundary conditions, numerical simulations have revealed that confining stresses tend to reduce fracture apertures and flow magnitude [23, 28], while large differential stresses often cause strong sliding along preferentially oriented rough fractures, which dilates fracture aperture and leads to permeability enhancement and flow channelization [7, 14, 25, 26, 28, 30]. On the other hand, the propagation of new cracks may also generate changes to the bulk flow properties through modifying the connectivity of fracture networks. This effect may be significant as natural fracture networks are often found to be close to the percolation threshold, i.e., critically connected [31, 32].In reviewing the literature, it is found that little effort has been devoted to quantifying the changes of flow structures in critically connected fracture systems due to stress loading and understanding the mechanisms of geomechanical effects that alter fluid flow. In particular, which geomechanical process (e.g., fracture dilation or propagation) dominates flow structure alteration remains an open question. In this work, we perform a generic study using high-fidelity numerical simulations with representative parameters to explore how the fluid flow properties change with in situ stresses for a natural fracture network whose connectivity state is close to the percolation threshold and further elucidate the mechanisms underpinning the changes. The remainder of the paper is organized as follows. In Section2, we present the numerical methods used in this work for modeling fracture network geometry, geomechanical deformation, and fluid flow of a naturally fractured rock. The set-up and boundary conditions of numerical experiments are presented in Section 3. The numerical simulation results are given in Section 4 with an emphasis on comparing the relative change in flow magnitude and organization caused by shear-induced dilation and fracture propagation. Finally, a brief discussion and conclusions are given in Section 5.


## 2. Models and Methods

### 2.1. Natural Fracture Network

The natural fracture network used in our simulations was mapped on a field outcrop of the Devonian sandstone at the Hornelen Basin, Norway (Figure1(a)) [33]. Due to direct field measurements for a complete 3D characterization which is difficulty to obtain, the sampled fracture patterns are limited to 2D. The fracture network covering an area of 18m×18m consists three major fracture sets of mean set orientation of 5°, 50°, and 120° (Figures 1(b)–1(d)). This 2D natural fracture network has a connectivity close to the percolation threshold such that the system is in a critical state in between well-connected and disconnected [34, 35]. The fracture network has been studied extensively in the past to investigate the geometrical organization and distribution of fractures [2, 34, 36] and their influence on fluid flow [8, 37] and mass transport [5, 33] properties of fractured rocks. In this paper, we use this fracture network to further examine the impact of in situ stresses on the flow structure.Figure 1
(a) The location of the outcrop. (b) A natural fracture network mapped at Hornelen Basin, Norway [33]. (b–d) Fracture statistics.
(a)
(b)
(c)
(d)


### 2.2. Geomechanical Model

The geomechanical deformation of the fracture network in response to in situ stresses is simulated using a hybrid finite-discrete element method [38]. The geomechanical model can realistically capture the deformation of intact rock, interaction of matrix blocks, variability of local stresses, displacement of preexisting fractures, and propagation of new cracks [14, 27]. The closure of rock fractures under compression is calculated based on a hyperbolic relation (Bandis et al., 1983):
(1)vn=σnvmkn0vm+σn,where vn is the normal closure, σn is the effective normal compressive stress, kn0 is the initial normal stiffness, and vm is the maximum allowable closure. The dependency of the shear behaviour of fractures on the normal stress loading is described using a constant displacement model parameterised with fixed up and ur values [39]. In this model, the peak shear stress τp is given by [40]:
(2)τp=σntanϕb+ϕi1−as+asc,where as is the proportion of total fracture area sheared through asperities, i is the dilation angle, c is the shear strength of the asperity (i.e., cohesion of the intact rock), and ϕb is the basic friction angle which is substituted using the residual friction angle ϕr [41]. If σn does not exceed the uniaxial compressive strength of the intact rock σu, the values of as and ϕi are, respectively, given as [40]:
(3)as=1−1−σnσum1,tanϕi=1−σnσum2tanϕi0,foru≤ur,0,foru>ur,where ϕi0 is the initial dilation angle when σn=0, and m1 and m2 are empirical parameters with suggested values of 1.5 and 4.0, respectively. The residual shear stress τr is given as [41]:
(4)τr=σntanϕr.The dilational displacementvs is related to the shear displacement u in an incremental form as [40]:
(5)dvs=−tanϕidu.The fracture apertureh under coupled normal and shear loadings is thus given by [10]:
(6)h=h0+w,σn<0,h0−vn−vs,σn≥0,where h0 is the initial aperture, and w is the separation of opposing fracture walls if the fracture is under tension.


### 2.3. Fluid Flow Model

We model a single-phase steady-state flow of incompressible fluids through the fractured rock (formed by a fracture population and permeable matrix) by solving the Laplace equation,(7)∇∙k∇p=0,where k is the intrinsic permeability and p is the fluid pressure. The model domain is discretized by an unstructured mesh with 2D triangular matrix elements and 1D lower dimensional fracture elements. The pressure field is solved using the finite element method. The flow velocity at the barycenter of each element is computed based on the pressure gradient vector field by applying Darcy’s law:
(8)ue=−keμ∇pe,while ue is the velocity vector field, pe is the local pressure at element nodes, μ is the dynamic fluid viscosity, and ke is the local permeability. In the computational mesh, the local permeability is constant for all matrix elements while spatially variable for fracture elements. The fracture permeabilities are determined according to the cubic law [42] based on the stress-dependent fracture apertures.


### 2.4. Quantitative Indices for Flow Field Evaluation

To quantify the stress effects on fluid flow, we compute three network-scale indices: (i) the network equivalent permeability, (ii) the correlation dimension of the flow rate field, and (iii) the flow channeling density indicator. The network equivalent permeability provides a measure of the overall flow capacity [43, 44] of the fracture network, while the correlation dimension and flow channeling density indicator characterize how the flow is spatially organized in the system [30, 45–48].The network equivalent permeability is obtained using Darcy’s law:(9)keff=μQLApin−pout,where Q is the volumetric flow rate through the inlet face of the fractured rock, L is the length of the model domain, A is the outlet cross-section area, and pin and pout are the hydraulic head at inlet and outlet, respectively.The correlation dimension of the flow rate fieldD2 is defined as [49, 50]:
(10)D2=limr→0log∑k=1MPk2logr,where Pk is the proportion of flow within an elemental area to the cumulative amount of flow in the entire grid, and r is the dimension of the grid elements. In this work, we use a 30 by 30 square grid in our calculation of D2, thus r equals to 0.6 m.The flow channeling density indicatordQ, similar to the participation ratio in the physics literature [51, 52], is defined as [46, 47]:
(11)dQ=1S∙∑fLf∙Qf2∑fLf∙Qf2,where Qf is the fracture flow rate, Lf is the fracture length, S is the cumulative fracture length of the entire network. The inverse of dQ defines the average spacing of main flow paths [46]. A small dQ indicates a large distance between main flow paths, and thus a highly channeled flow pattern.


## 2.1. Natural Fracture Network

The natural fracture network used in our simulations was mapped on a field outcrop of the Devonian sandstone at the Hornelen Basin, Norway (Figure1(a)) [33]. Due to direct field measurements for a complete 3D characterization which is difficulty to obtain, the sampled fracture patterns are limited to 2D. The fracture network covering an area of 18m×18m consists three major fracture sets of mean set orientation of 5°, 50°, and 120° (Figures 1(b)–1(d)). This 2D natural fracture network has a connectivity close to the percolation threshold such that the system is in a critical state in between well-connected and disconnected [34, 35]. The fracture network has been studied extensively in the past to investigate the geometrical organization and distribution of fractures [2, 34, 36] and their influence on fluid flow [8, 37] and mass transport [5, 33] properties of fractured rocks. In this paper, we use this fracture network to further examine the impact of in situ stresses on the flow structure.Figure 1
(a) The location of the outcrop. (b) A natural fracture network mapped at Hornelen Basin, Norway [33]. (b–d) Fracture statistics.
(a)
(b)
(c)
(d)


## 2.2. Geomechanical Model

The geomechanical deformation of the fracture network in response to in situ stresses is simulated using a hybrid finite-discrete element method [38]. The geomechanical model can realistically capture the deformation of intact rock, interaction of matrix blocks, variability of local stresses, displacement of preexisting fractures, and propagation of new cracks [14, 27]. The closure of rock fractures under compression is calculated based on a hyperbolic relation (Bandis et al., 1983):
(1)vn=σnvmkn0vm+σn,where vn is the normal closure, σn is the effective normal compressive stress, kn0 is the initial normal stiffness, and vm is the maximum allowable closure. The dependency of the shear behaviour of fractures on the normal stress loading is described using a constant displacement model parameterised with fixed up and ur values [39]. In this model, the peak shear stress τp is given by [40]:
(2)τp=σntanϕb+ϕi1−as+asc,where as is the proportion of total fracture area sheared through asperities, i is the dilation angle, c is the shear strength of the asperity (i.e., cohesion of the intact rock), and ϕb is the basic friction angle which is substituted using the residual friction angle ϕr [41]. If σn does not exceed the uniaxial compressive strength of the intact rock σu, the values of as and ϕi are, respectively, given as [40]:
(3)as=1−1−σnσum1,tanϕi=1−σnσum2tanϕi0,foru≤ur,0,foru>ur,where ϕi0 is the initial dilation angle when σn=0, and m1 and m2 are empirical parameters with suggested values of 1.5 and 4.0, respectively. The residual shear stress τr is given as [41]:
(4)τr=σntanϕr.The dilational displacementvs is related to the shear displacement u in an incremental form as [40]:
(5)dvs=−tanϕidu.The fracture apertureh under coupled normal and shear loadings is thus given by [10]:
(6)h=h0+w,σn<0,h0−vn−vs,σn≥0,where h0 is the initial aperture, and w is the separation of opposing fracture walls if the fracture is under tension.


## 2.3. Fluid Flow Model

We model a single-phase steady-state flow of incompressible fluids through the fractured rock (formed by a fracture population and permeable matrix) by solving the Laplace equation,(7)∇∙k∇p=0,where k is the intrinsic permeability and p is the fluid pressure. The model domain is discretized by an unstructured mesh with 2D triangular matrix elements and 1D lower dimensional fracture elements. The pressure field is solved using the finite element method. The flow velocity at the barycenter of each element is computed based on the pressure gradient vector field by applying Darcy’s law:
(8)ue=−keμ∇pe,while ue is the velocity vector field, pe is the local pressure at element nodes, μ is the dynamic fluid viscosity, and ke is the local permeability. In the computational mesh, the local permeability is constant for all matrix elements while spatially variable for fracture elements. The fracture permeabilities are determined according to the cubic law [42] based on the stress-dependent fracture apertures.


## 2.4. Quantitative Indices for Flow Field Evaluation

To quantify the stress effects on fluid flow, we compute three network-scale indices: (i) the network equivalent permeability, (ii) the correlation dimension of the flow rate field, and (iii) the flow channeling density indicator. The network equivalent permeability provides a measure of the overall flow capacity [43, 44] of the fracture network, while the correlation dimension and flow channeling density indicator characterize how the flow is spatially organized in the system [30, 45–48].The network equivalent permeability is obtained using Darcy’s law:(9)keff=μQLApin−pout,where Q is the volumetric flow rate through the inlet face of the fractured rock, L is the length of the model domain, A is the outlet cross-section area, and pin and pout are the hydraulic head at inlet and outlet, respectively.The correlation dimension of the flow rate fieldD2 is defined as [49, 50]:
(10)D2=limr→0log∑k=1MPk2logr,where Pk is the proportion of flow within an elemental area to the cumulative amount of flow in the entire grid, and r is the dimension of the grid elements. In this work, we use a 30 by 30 square grid in our calculation of D2, thus r equals to 0.6 m.The flow channeling density indicatordQ, similar to the participation ratio in the physics literature [51, 52], is defined as [46, 47]:
(11)dQ=1S∙∑fLf∙Qf2∑fLf∙Qf2,where Qf is the fracture flow rate, Lf is the fracture length, S is the cumulative fracture length of the entire network. The inverse of dQ defines the average spacing of main flow paths [46]. A small dQ indicates a large distance between main flow paths, and thus a highly channeled flow pattern.


## 3. Model Set-Up

In this study, the stress-dependent fracture aperture fields are obtained by imposing effective far-field stresses to the fracture network at various angles ranging from 0° to 170° with an increment of 10° (Figure2(a)). We consider five different stress scenarios: (i) Smax=Smin=0.0 MPa (reference case); (ii) Smax=Smin=5.0 MPa; (iii) Smax=10.0 MPa, Smin=5.0 MPa; (iv) Smax=15.0 MPa, Smin=5.0 MPa; and (v) Smax=20.0 MPa, Smin=5.0 MPa. The initial aperture field is assumed constant as h0=0.1 mm, which is representative for nature fractures [53–55].Figure 2
(a) Far-field stress boundary conditions for geomechanical simulation. (b) Hydraulic boundary conditions for fluid flow simulation.
(a)
(b)We simulate a single-phase steady-state flow through the deformed fracture networks by applying the classical permeameter-type boundary condition (Figure2(b)): two opposite model boundaries have constant hydraulic heads, which create a fixed pressure gradient (i.e., 10 kPa), while the two orthogonal boundaries parallel to the flow direction are impervious. The matrix permeability km is assumed constant in space with a value of 1×10−15 m2. The material properties of the fractured rocks and fluids used in the simulations (Table 1) are typical for sandstone formations [5]. In this work, we focus on discussing the cases where the flow is along the x-direction; the results for y-direction flow are given in the supplementary document (available here).Table 1
Material properties for geomechanical computation.
Properties
Value
Unit
Density
2500
kg·m-3
Young’s modulus
20
GPa
Poisson’s ratio
0.25
—
Internal friction coefficient
0.85
—
Intact rock cohesion
40
MPa
Tensile strength
20
MPa
Mode I energy release rate
396
J/m2
Mode II energy release rate
495
J/m2
Fracture initial normal stiffness
50
GPa/m
Fracture residual friction angle
30
deg
Fracture initial dilation angle
10
deg
Fluid density
1000
kg/m3
Dynamic viscosity
0.001
Pa·s
Rock density
2650
kg/m3
Gravity acceleration
9.8
m/s2
Initial fracture aperture
0.1
mm
Matrix permeability
1.0×10−15
m2


## 4. Results

We analyze the geomechanical responses of the fractured rocks under different far-field stress conditions. Figure3 shows the distributions of fracture apertures in the fractured rocks under various stress conditions for representative cases of θ=0°, 30°, 60°, 90°, 120°, and 150°. The probability density functions (PDFs) of these cases are given in Figure 4.Figure 3
Distributions of fracture apertures in the fracture network under different stress conditions.Figure 4
Probability density function of fracture apertures in the fracture network under various stress conditions.WhenSmax=Smin=0 MPa, all fractures have an identical aperture, i.e., the prescribed initial value of 0.1 mm (Figures 3 and 4). If the fractured rock is isotropically stressed with Smax=Smin=5.0 MPa, all fractures tend to be uniformly compressed such that the deformed apertures are lower than the initial value. If the fractured rock is anisotropically stressed, the fracture network becomes to accommodate heterogeneous aperture distribution as a result of the superimposed effects of shear-induced dilation and compression-induced closure (Figures 3 and 4). When Smax=10 MPa and Smin=5 MPa, most fractures are closed by compression, while only a few small fractures exhibit apertures larger than the initial value of 0.1 mm (Figures 3 and 4). When Smax=15 MPa and Smin=5 MPa, the shear-induced dilation of fracture apertures is more prevailing along some long fractures (Figure 3). When Smax is further increased to 20 MPa, the number of fractures with large apertures greatly increases due to enhanced shear displacements driven by elevated differential stresses (Figure 3). The fractures’ shear behavior is strongly controlled by the relative angle between the orientation of far-field stresses and the mean orientation of fracture sets. For example, when the rotation angle of far-field stresses equals to 30°, the 0° and 50° fracture sets exhibit enhanced shear displacements, while the 120° set is more suppressed for shearing because it is oriented almost perpendicularly to the maximum principal stress, Smax. The anisotropic geomechanical response of the fractured rock is further revealed by the different shape of aperture PDFs for different orientations of far-field stress field (Figure 4).We obtain the flow velocity fields of the deformed fractured rocks by solving the single-phase steady-state flow equations (i.e., Equations (7) and (8)). Figure 5 displays the flow fields for various stress loading cases. The PDFs of fracture flow velocities are provided in Figure 6. As shown in Figure 5, the flow is highly channelized in a subnetwork consisting mainly of fractures from the 5° and the 30° sets, if zero far-field stresses are applied. When the fractured rock is isotropically loaded (i.e., Smax=Smin=5 MPa), the flow velocity in the fracture network is systematically reduced due to the compression-induced normal closure of fractures (Figures 5 and 6). If the fractured rock is loaded with the anisotropic stress field of Smax=10 MPa and Smin=5 MPa, the main flow network still shows an overall reduction in flow velocity due to compression-induced aperture closure (Figures 5 and 6), while locally some fractures exhibit increased velocity values due to moderate shear-induced aperture dilation (Figure 5). The flow velocity fields for various far-field stress orientations are similar (Figure 5). This similarity of flow distribution is also revealed by the resembled shape of PDFs for various θ cases (Figure 6). When Smax=15 MPa and Smin=5 MPa, the flow fields become more heterogeneous and sensitive to the variation in θ (Figure 5). Visible variations of flow field, reflected by differences in PDF shape and range, with far-field stress orientation occurs (Figure 6). At specific orientations of far-field stress field, e.g., 30° and 60°, we observe a conspicuous increase in flow velocity in some fracture clusters in the lower-right part of the fracture network, which is not a part of the original main flow network under zero stress loadings (Figure 5). When Smax is further increased to 20 MPa, a stronger dependence of the flow field on the far-field stress orientation is observed (Figure 5). The extent of PDFs’ spreading towards high velocities varies with stress (Figure 6). It seems that the flow field alteration under the two high stress ratio conditions (i.e., Smax=15 or 20 MPa) is driven by a different mechanism compared with the isotropic stress cases and anisotropic cases with a low stress contrast (i.e., the case of Smax=10 MPa and Smin=5 MPa). Indeed, we observe numerous new cracks formed in the cases of high differential stresses (Figure 7). These propagated fractures, although small in size, may create new flow paths connecting preexisting disconnected clusters in the fractured rock. To further elucidate this effect, we plot flow distributions in fractures intersecting the outlet, i.e., the right model boundary (Figure 8). We note that a change in the relative flow magnitude of effluent fractures may indicate an alteration of internal flow organization. When Smax=10 MPa and Smin=5 MPa, there is only a quantitative variation in the relative magnitude of fracture flow when θ is varied; however, when Smax=15 MPa and Smin=5 MPa or Smax=20 MPa and Smin=5 MPa, except for a more drastic variation in flow magnitude with stress, the locations of flowing fractures are also changed qualitatively as θ varies (Figure 8). We will further analyze which geomechanical effect, i.e., shear-induced dilation or new crack propagation, dominates this enhanced sensitivity to stress towards the end of this section.Figure 5
Distributions of flow velocity in the fracture network under different stress conditions.Figure 6
Probability density function of the flow velocity in the fracture network under various stress conditions.Figure 7
Formation of new cracks under high differential stress conditions. (a) All fractures; (b) location and aperture of new cracks.Figure 8
Flow rate distributions at the outlet of fracture networks under different stress condition.To quantify the impact of stress variation on flow magnitude and organization, we further derive three flow indicators, i.e., the equivalent network permeabilitykeff, the correlation dimension of flow rate D2, and the flow channeling density indicator dQ, for various stress conditions. The equivalent permeability may quantify the bulk flow magnitude through the fractured rock, while D2 and dQ further indicate the flow distribution within the fracture network.As shown in Figure9(a), if the fractured rock is isotropically loaded (Smax=Smin=0 or 5 MPa), keff does not vary with the stress field rotation, and neither does the flow distribution (Figures 9(b) and 9(c)). As the stress level increases (i.e., in the case of Smax=5 MPa and Smin=5 MPa), keff decreases by about several times due to compression-induced aperture closure (Figure 3). On the other hand, both D2 and dQ increase, indicating a less channelized flow distribution in the fractured rock. If the fractured rock is anisotropically stressed, keff varies with the far-field stress orientation. When Smax=10 MPa and Smin=5 MPa, the keff becomes larger than that of the isotropic stress loading cases of Smax=5 MPa and Smin=5 MPa but is still smaller than that of the zero-stress loading case (i.e., Smax=0 MPa and Smin=0 MPa) due to the dominant normal closure effect under confining stresses. The variation of keff with far-field stress orientation is quite small. A variation in D2 and dQ also occurs. This indicates that the flow structure becomes dependent on the far-field stress orientation. With the increase of the differential stress magnitude (i.e., when Smax=15 MPa and Smin=5 MPa), the keff variation as a function of far-field stress orientation becomes very complex under the combined effects of compression-induced closure, shear-induced dilation, and formation of new cracks. In general, keff increases as the increase of differential stress so far but is still smaller than the that of the fractured rocks under the zero-stress loading. Moreover, D2 and dQ also become smaller, indicating a higher level of flow channelization. The variation magnitude of D2 and dQ with stress becomes more prominent. The variation trends of D2 and dQ as the far-field stress field orientation appear similar (Figures 9(b) and 9(c)). As Smax is further increased to 20 MPa, there is a significant overall increase in keff (Figure 9(a)). For all far-field stress orientations, the keff values are now larger than that of the fractured rock under the zero-stress loading. Both D2 and dQ exhibit an increased variability with far-field stress orientation. Moreover, their variation trends are no longer similar. When θ falls in the range between 70° and 90°, the D2 and dQ values remain high and similar to those of the case of Smax=10 MPa and Smin=5 MPa. When θ=0°, 130°, and 140°, an enhanced reduction in dQ is observed (Figures 9(b) and 9(c)), while there is a slightly less decrease in D2. It seems that the cases with high bulk permeability have less localized flow, although this correlation is associated with large uncertainties.Figure 9
Quantitative indices for evaluating flow magnitude and distribution in fracture networks under different stress conditions: (a) effective permeability, (b) correlation dimension of flow rate, and (c) flow channeling density indicator.
(a)
(b)
(c)To further investigate the relative contribution of shear displacement and new fracture propagation to the permeability enhancement and flow structure alternation under high differential stresses (Figure7), we show, in Figure 10, the derived keff, D2, and dQ values after removing all newly propagated cracks from the stressed fracture network. By comparing Figure 10 against Figure 9, it can be seen that the bulk permeability decreases by several times with the removal of new cracks. In addition, there is also a systematic decrease in both D2 and dQ, implying a higher level of flow channelization without the new cracks. These results are consistent with our observations of new flow paths created by the propagated cracks (Figures 5 and 8). The comparison between Figures 9 and 10 further suggests that the permeability enhancement may be more dominated by the propagation of new fractures while shear-induced dilation plays a second-order role. It is also evident that the elevated flow channelization due to shear displacement has been attenuated by the presence of new flow paths created by propagated cracks. This phenomenon is considered to be related to the fact that the natural fracture network is close to the percolation threshold such that new cracks can critically transition the system from disconnected to connected (see more discussions in the following section).Figure 10
Quantitative indices for evaluating flow magnitude and distribution in fracture networks under different stress condition after removing the new cracks: (a) effective permeability, (b) correlation dimension of flow rate, and (c) flow channeling density indicator.
(a)
(b)
(c)


## 5. Discussions and Conclusions

In this study, we have demonstrated how in situ stress variations impact fluid flow in a 2D natural fracture network that is critically connected (i.e., the network connectivity is close to the percolation threshold). We have shown that the stress variation may induce changes in the bulk permeability and flow organization. As the stress magnitude increases, fracture apertures decrease systematically, which in turn reduces fracture flow velocities and the bulk permeability. For anisotropic stress conditions, significant shear-induced dilation may occur along preferentially oriented fractures with respect to the anisotropic stress field. Although the shear displacement may cause a significant change in the internal flow organization (generally enhances the flow channelization intensity), its impact on the network’s macroscopic hydraulic properties is mild. These observations are consistent with previous studies where similar results about the stress-dependent permeability of 2D fracture networks were reported [7, 11, 14, 23, 26, 28]. We further show that under high differential stress conditions, the formation of new cracks is also prevailing in the fracture network, which may change the number and location of preferential flow paths. As a result, the bulk permeability may increase significantly. The new cracks in general reduce the level of flow channelization within the network as they may bridge unconnected fractures or clusters to form new paths branching the flow. This effect is particularly prominent when the fractured rock is critically stressed (i.e., with a stress ratio larger than 3).These observations indicate that in natural fracture networks close to the percolation threshold, the network connectivity has a primary control on fluid flow, while the shear-related geomechanical effects only pose a second-order impact on the hydraulic properties of the fractured rock. However, the effect of fracture propagation still reserves the power to regulate the hydraulic connectivity of the fracture network. These research findings highlight the importance of integrating geomechanical analyses in the practice of understanding and prediction of fluid flow in natural fracture systems as they are mostly critically connected [32] and under critical stress condition [56]. Our numerical simulations clearly show that without taking into account the stress effects, we may largely underestimate the level of flow channelization and overestimate the bulk permeability of naturally fractured rocks subjecting to anisotropic loadings with intermediate stress magnitude. However, when the fractured rock is loaded by anisotropic stress of high magnitudes, we tend to underestimate the bulk permeability and overestimate the flow channeling intensity without considering the fracture propagation process. Moreover, the reduction in hydraulic connection caused by stress orientation variation may also be overlooked, which may cause complexities or even the failure of the production scheme under implementation, since stress redistributions often occur during the production of geofluids.Although the present work provides a better understanding of hydromechanical behavior of a nature fracture network under different stress conditions, particularly the evolution of flow paths within situ stresses, a number of issues remain to be addressed in future work. First, the present work needs to be extended to 3D to further verify the consistency and generality of the main results discovered based on 2D simulations. The geometry and connectivity of 3D fracture networks are intricate (Bour and Davy, 1998). The stress effects are thus expected to be much more complex. Based on the findings from our recent 3D modeling work (e.g., [27]), the 2D analysis presented here may provide some indicative approximations, but certainly cannot reveal fully the polyaxial stress-dependent behavior of 3D fracture networks. However, due to the very expensive run time of 3D FEMDEM simulations, the number of fractures in the 3D model has to be very limited. The challenges in performing large-scale simulations need to be addressed before they can be used to drive pertinent 3D characterizations of stress effects at the scale of natural fracture systems. Moreover, our work only examined the hydromechanical behavior of a specific natural fracture network. The validity and universality of the findings for other natural fracture systems require further investigations. Furthermore, spatial variability and correlation of initial fracture aperture may also need to be integrated into the model to more realistically simulate the actual systems.



---

*Source: 3138972-2019-12-16.xml*