<figure><img src="image_1.png" /></figure>

HAL
open science

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

Chuanyin Jiang, Xiaoguang Wang, Zhixue Sun, Qinghua Lei

► To cite this version:

Chuanyin Jiang, Xiaoguang Wang, Zhixue Sun, Qinghua Lei. The Role of In Situ Stress in Organizing Flow Pathways in Natural Fracture Networks at the Percolation Threshold. Geofluids, 2019, 2019, pp.1-14. 10.1155/2019/3138972. hal-03168260

HAL Id: hal-03168260

https://hal.science/hal-03168260

Submitted on 14 Mar 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Research Article

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

Chuanyin Jiang,<fnref n="1" /> Xiaoguang Wang ID,<fnref n="2" /> Zhixue Sun,<fnref n="1" />,<fnref n="2" /> and Qinghua Lei<fnref n="3" />

Correspondence should be addressed to Xiaoguang Wang; x.wang0127@gmail.com

Received 24 July 2019; Accepted 13 November 2019; Published 16 December 2019

Guest Editor: Tanguy Robert

Copyright © 2019 Chuanyin Jiang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

## 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

<footnote n="1">¹School of Petroleum Engineering, China University of Petroleum (East China), Qingdao, China</footnote>

<footnote n="2">²Hydrosciences Montpellier, UMR 5569 CNRS-IRD-UM, Montpellier, France</footnote>

<footnote n="3">³Department of Earth Sciences, ETH Zurich, Switzerland</footnote>
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 the *in 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 Section 2, 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 (Figure 1(a)) [33]. Due to direct field measurements for a complete 3D characterization which is difficult to obtain, the sampled fracture patterns are limited to 2D. The fracture network covering an area of 18 m × 18 m 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.

### 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):

$$v_n = \frac{\sigma_n v_m}{k_{n0} v_m + \sigma_n}, \quad (1)$$

where $v_n$ is the normal closure, $\sigma_n$ is the effective normal compressive stress, $k_{n0}$ is the initial normal stiffness, and $v_m$ 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 $u_p$ and $u_r$ values [39]. In this model, the peak shear stress $\tau_p$ is given by [40]:

$$\tau_p = \sigma_n \tan (\phi_b + \phi_i)(1 - a_s) + a_s c, \quad (2)$$

where $a_s$ 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 $\phi_b$ is the basic friction angle which is substituted using the residual friction angle $\phi_r$ [41]. If $\sigma_n$ does not exceed the uniaxial compressive strength of the intact rock $\sigma_u$, the values of $a_s$ and $\phi_i$ are, respectively, given as [40]:

$$a_s = 1 - \left(1 - \frac{\sigma_n}{\sigma_u}\right)^{m_1}, \\ \tan \phi_i = \begin{cases} \left(1 - \frac{\sigma_n}{\sigma_u}\right)^{m_2} \tan \phi_{i0}, & \text{for } u \le u_r, \\ 0, & \text{for } u > u_r, \end{cases} \quad (3)$$

where $\phi_{i0}$ is the initial dilation angle when $\sigma_n = 0$, and $m_1$ and $m_2$ are empirical parameters with suggested values of 1.5 and 4.0, respectively. The residual shear stress $\tau_r$ is given as [41]:

$$\tau_r = \sigma_n \tan \phi_r. \quad (4)$$

The dilational displacement $v_s$ is related to the shear displacement $u$ in an incremental form as [40]:

$$dv_s = - \tan \phi_i du. \quad (5)$$
<figure><img src="image_3.png" /><figcaption>FIGURE 1: (a) The location of the outcrop. (b) A natural fracture network mapped at Hornelen Basin, Norway [33]. (b–d) Fracture statistics.</figcaption></figure>

The fracture aperture h under coupled normal and shear loadings is thus given by [10]:

$$h = \begin{cases} h_0 + w, & \sigma_n < 0, \\ h_0 - v_n - v_s, & \sigma_n \ge 0, \end{cases} \quad (6)$$

where $h_0$ 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,

$$\nabla \cdot (k \nabla p) = 0, \quad (7)$$

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:

$$\mathbf{u}^e = -\frac{k^e}{\mu} \nabla p^e, \quad (8)$$

while $\mathbf{u}^e$ is the velocity vector field, $p^e$ is the local pressure at element nodes, $\mu$ is the dynamic fluid viscosity, and $k^e$ 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].
<figure><img src="image_4.png" /><figcaption>FIGURE 2: (a) Far-field stress boundary conditions for geomechanical simulation. (b) Hydraulic boundary conditions for fluid flow simulation.</figcaption></figure>

The network equivalent permeability is obtained using Darcy's law:

$$k_{\text{eff}} = \frac{\mu Q L}{A(p_{\text{in}} - p_{\text{out}})}, \qquad (9)$$

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 $p_{\text{in}}$ and $p_{\text{out}}$ are the hydraulic head at inlet and outlet, respectively.

The correlation dimension of the flow rate field $D_2$ is defined as [49, 50]:

$$D_2 = \lim_{r \to 0} \frac{\log \sum_{k=1}^{M} P_k^2}{\log r}, \qquad (10)$$

where $P_k$ 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 $D_2$, thus r equals to 0.6 m.

The flow channeling density indicator $d_Q$, similar to the participation ratio in the physics literature [51, 52], is defined as [46, 47]:

$$d_Q = \frac{1}{S} \cdot \frac{(\sum_f L_f \cdot Q_f)^2}{(\sum_f L_f \cdot Q_f^2)}, \qquad (11)$$

where $Q_f$ is the fracture flow rate, $L_f$ is the fracture length, S is the cumulative fracture length of the entire network. The inverse of $d_Q$ defines the average spacing of main flow paths [46]. A small $d_Q$ 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° (Figure 2(a)). We consider five

<table><caption>TABLE 1: Material properties for geomechanical computation.</caption><thead><tr><th>Properties</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Density</td><td>2500</td><td>kg·m<sup>-3</sup></td></tr><tr><td>Young's modulus</td><td>20</td><td>GPa</td></tr><tr><td>Poisson's ratio</td><td>0.25</td><td>—</td></tr><tr><td>Internal friction coefficient</td><td>0.85</td><td>—</td></tr><tr><td>Intact rock cohesion</td><td>40</td><td>MPa</td></tr><tr><td>Tensile strength</td><td>20</td><td>MPa</td></tr><tr><td>Mode I energy release rate</td><td>396</td><td>J/m<sup>2</sup></td></tr><tr><td>Mode II energy release rate</td><td>495</td><td>J/m<sup>2</sup></td></tr><tr><td>Fracture initial normal stiffness</td><td>50</td><td>GPa/m</td></tr><tr><td>Fracture residual friction angle</td><td>30</td><td>deg</td></tr><tr><td>Fracture initial dilation angle</td><td>10</td><td>deg</td></tr><tr><td>Fluid density</td><td>1000</td><td>kg/m<sup>3</sup></td></tr><tr><td>Dynamic viscosity</td><td>0.001</td><td>Pa·s</td></tr><tr><td>Rock density</td><td>2650</td><td>kg/m<sup>3</sup></td></tr><tr><td>Gravity acceleration</td><td>9.8</td><td>m/s<sup>2</sup></td></tr><tr><td>Initial fracture aperture</td><td>0.1</td><td>mm</td></tr><tr><td>Matrix permeability</td><td>1.0 × 10<sup>-15</sup></td><td>m<sup>2</sup></td></tr></tbody></table>

different stress scenarios: (i) $S_{\text{max}} = S_{\text{min}} = 0.0$ MPa (reference case); (ii) $S_{\text{max}} = S_{\text{min}} = 5.0$ MPa; (iii) $S_{\text{max}} = 10.0$ MPa, $S_{\text{min}} = 5.0$ MPa; (iv) $S_{\text{max}} = 15.0$ MPa, $S_{\text{min}} = 5.0$ MPa; and (v) $S_{\text{max}} = 20.0$ MPa, $S_{\text{min}} = 5.0$ MPa. The initial aperture field is assumed constant as $h_0 = 0.1$ mm, which is representative for nature fractures [53–55].

We simulate a single-phase steady-state flow through the deformed fracture networks by applying the classical permeameter-type boundary condition (Figure 2(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 $k_m$ is assumed constant in space with a value of $1 \times 10^{-15}$ m$^2$. 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
<figure><img src="image_5.png" /><figcaption>FIGURE 3: Distributions of fracture apertures in the fracture network under different stress conditions.</figcaption></figure>

along the x-direction; the results for y-direction flow are given in the supplementary document (available here).

**4. Results**

We analyze the geomechanical responses of the fractured rocks under different far-field stress conditions. Figure 3 shows the distributions of fracture apertures in the fractured rocks under various stress conditions for representative cases of $\theta = 0°, 30°, 60°, 90°, 120°$, and $150°$. The probability density functions (PDFs) of these cases are given in Figure 4.

When $S_{max} = S_{min} = 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 $S_{max} = S_{min} = 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 $S_{max} = 10$ MPa and $S_{min} = 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 $S_{max} = 15$ MPa and $S_{min} = 5$ MPa, the shear-induced dilation of fracture apertures is more prevailing along some long fractures (Figure 3). When $S_{max}$ 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, $S_{max}$. 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., $S_{max} = S_{min} = 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 $S_{max} = 10$ MPa and $S_{min} = 5$ MPa, the main flow network still shows an overall reduction in flow velocity due to compression-induced aperture closure
<figure><img src="image_6.png" /><figcaption>FIGURE 4: Probability density function of fracture apertures in the fracture network under various stress conditions.</figcaption></figure>

(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 S_max = 15 MPa and S_min = 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 S_max 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
<figure><img src="image_7.png" /><figcaption>FIGURE 5: Distributions of flow velocity in the fracture network under different stress conditions.</figcaption></figure>

conditions (i.e., $S_{max}$ = 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 $S_{max}$ = 10 MPa and $S_{min}$ = 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 $S_{max}$ = 10 MPa and $S_{min}$ = 5 MPa, there is only a quantitative variation in the relative magnitude of fracture flow when $\theta$ is varied; however, when $S_{max}$ = 15 MPa and $S_{min}$ = 5 MPa or $S_{max}$ = 20 MPa and $S_{min}$ = 5 MPa, except for a more drastic variation in flow magnitude with stress, the locations of flowing fractures are also changed qualitatively as $\theta$ 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.

To quantify the impact of stress variation on flow magnitude and organization, we further derive three flow indicators, i.e., the equivalent network permeability $k_{eff}$, the correlation dimension of flow rate $D_2$, and the flow channeling density indicator $d_Q$, for various stress condi-

tions. The equivalent permeability may quantify the bulk flow magnitude through the fractured rock, while $D_2$ and $d_Q$ further indicate the flow distribution within the fracture network.

As shown in Figure 9(a), if the fractured rock is isotropically loaded ($S_{max}$ = $S_{min}$ = 0 or 5 MPa), $k_{eff}$ 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 $S_{max}$ = 5 MPa and $S_{min}$ = 5 MPa), $k_{eff}$ decreases by about several times due to compression-induced aperture closure (Figure 3). On the other hand, both $D_2$ and $d_Q$ increase, indicating a less channelized flow distribution in the fractured rock. If the fractured rock is anisotropically stressed, $k_{eff}$ varies with the far-field stress orientation. When $S_{max}$ = 10 MPa and $S_{min}$ = 5 MPa, the $k_{eff}$ becomes larger than that of the isotropic stress loading cases of $S_{max}$ = 5 MPa and $S_{min}$ = 5 MPa but is still smaller than that of the zero-stress loading case (i.e., $S_{max}$ = 0 MPa and $S_{min}$ = 0 MPa) due to the dominant normal closure effect under confining stresses. The variation of $k_{eff}$ with far-field stress orientation is quite small. A variation in $D_2$ and $d_Q$ 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 $S_{max}$ = 15 MPa and $S_{min}$ = 5 MPa), the $k_{eff}$ 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, $k_{eff}$ increases as the increase of differential stress so far but
<figure><img src="image_8.png" /><figcaption>FIGURE 6: Probability density function of the flow velocity in the fracture network under various stress conditions.</figcaption></figure>

is still smaller than the that of the fractured rocks under the zero-stress loading. Moreover, $D_2$ and $d_Q$ also become smaller, indicating a higher level of flow channelization. The variation magnitude of $D_2$ and $d_Q$ with stress becomes more prominent. The variation trends of $D_2$ and $d_Q$ as the far-field stress field orientation appear similar (Figures 9(b) and 9(c)). As $S_{max}$ is further increased to 20 MPa, there is a significant overall increase in $k_{eff}$ (Figure 9(a)). For all far-field stress orientations, the $k_{eff}$ values are now larger than that of the fractured rock under the zero-stress loading. Both $D_2$ and $d_Q$ exhibit an increased variability with far-field stress

orientation. Moreover, their variation trends are no longer similar. When $\theta$ falls in the range between 70° and 90°, the $D_2$ and $d_Q$ values remain high and similar to those of the case of $S_{max}$ = 10 MPa and $S_{min}$ = 5 MPa. When $\theta$ = 0°, 130°, and 140°, an enhanced reduction in $d_Q$ is observed (Figures 9(b) and 9(c)), while there is a slightly less decrease in $D_2$. It seems that the cases with high bulk permeability have less localized flow, although this correlation is associated with large uncertainties.

To further investigate the relative contribution of shear displacement and new fracture propagation to the
<figure><img src="image_9.png" /><figcaption>FIGURE 7: Formation of new cracks under high differential stress conditions. (a) All fractures; (b) location and aperture of new cracks.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 8: Flow rate distributions at the outlet of fracture networks under different stress condition.</figcaption></figure>

permeability enhancement and flow structure alternation under high differential stresses (Figure 7), we show, in Figure 10, the derived $k_{\text{eff}}$, $D_2$, and $d_Q$ 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 system-

atic decrease in both $D_2$ and $d_Q$, 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
<figure><img src="image_11.png" /><figcaption>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.</figcaption></figure>

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).

## 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 connec-

tivity 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
<figure><img src="image_12.png" /><figcaption>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.</figcaption></figure>

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 with *in 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.

## Data Availability

The simulation data used to support the findings of this study are available from the corresponding author upon request.

## Additional Points

*Highlights.* (1) Stress effect on flow organization in a natural fracture network is investigated. (2) The role of shear-induced dilation and stress-driven fracture on fluid flow is evaluated. (3) New cracks may serve as red links to critically enhance the network connectivity. (4) Shear-induced dilation tends to be a second-order factor on altering permeability

## Conflicts of Interest

The authors declare that they have no conflicts of interest.

## Acknowledgments

ZS acknowledges funding from the National Natural Science Foundation of China (No. 51774317), the Fundamental Research Funds for the Central Universities (No. 18CX2100A), and the National Science and Technology Major Project (No. 2016ZX05011004-004). All the data of numerical simulation results in this paper are available upon request to the correspondence author X. Wang (xiaoguang.wang@umontpellier.fr).

## Supplementary Materials

In this supplementary file, we provide flow simulation results in the *y*-direction for various stress cases. As shown in Figure S1, when calculating the flow field, a vertical hydraulic gradient is applied to the fracture networks, where the fracture apertures are derived from geomechanical simulations with zero, isotropic, and anisotropic far-field stress loadings. The material properties of the fracture networks and fluids used in these simulations are the same as listed in Table 1. Figure S2 displays the flow fields for various stress loading cases. Overall, the results in the *y*-direction (Figures S2, S3, S4, and S5) are similar to those in the *x*-direction (Figures 5, 6, 9, and 10), revealing that the effect of geomechanical deformation depends on the angle between maximum stress and flow direction. The *y*-direction flow is mainly restricted to the 50° and 120° fracture sets. The large extent of PDFs' spreading towards high velocities tends to occur when $\theta = 90°$ in the cases of *y*-direction, whereas it occurs when $\theta = 30°$ and 150° in *x*-direction (Figure S3 and Figure 6). Note that $\theta$ is the angle between maximum stress and *x* positive direction. The variation tendency of the quantitative flow indices ($k_{\text{eff}}/k_{\text{m}}$, $D_2$, and $d_Q$) under different stress conditions (Figures S4 and S5) is coincident with the cases in the *x*-direction (Figures 9 and 10), except that the corresponding transition points are shifted and the shape of the curves is changed compared to the respective case for *x*-direction flow. The difference between the results of *x*- and *y*-directions may be caused by the interplay among the maximum stress direction, mean flow direction, and fracture network geometry (e.g., the direction of the main fracture sets and the organization of fractures of different sets). (*Supplementary Materials*)

## References

[1] P. M. Adler and J.-F. Thovert, *Fractures and Fracture Networks, Theory and Applications of Transport in Porous Media*, Kluwer Academic Publishers, Dordrecht, Netherlands, 1999.

[2] B. Berkowitz, "Characterizing flow and transport in fractured geological media: a review," *Advances in Water Resources*, vol. 25, no. 8-12, pp. 861-884, 2002.

[3] S. P. Neuman, "Trends, prospects and challenges in quantifying flow and transport through fractured rocks," *Hydrogeology Journal*, vol. 13, no. 1, pp. 124-147, 2005.

[4] I. Miller, G. Lee, and W. Dershowitz, "MAFIC-matrix/fracture interaction code with heat and solute transport-user
documentation, version 1.6," Golder Associates Inc., Red-
mond, WA, USA, 1999.

[5] S. Geiger and S. Emmanuel, "Non-Fourier thermal transport in fractured geological media," *Water Resources Research*, vol. 46, no. 7, 2010.

[6] N. J. Hardebol, C. Maier, H. Nick, S. Geiger, G. Bertotti, and H. Boro, "Multiscale fracture network characterization and impact on flow: a case study on the Latemar carbonate platform," *Journal of Geophysical Research Solid Earth*, vol. 120, no. 12, pp. 8197-8222, 2015.

[7] Q. Lei, X. Wang, J. Xiang, and J.-P. Latham, "Polyaxial stress-dependent permeability of a three-dimensional fractured rock layer," *Hydrogeology Journal*, vol. 25, no. 8, pp. 2251-2262, 2017.

[8] N. E. Odling, P. Gillespie, B. Bourgine et al., "Variations in fracture system geometry and their implications for fluid flow in fractures hydrocarbon reservoirs," *Petroleum Geoscience*, vol. 5, no. 4, pp. 373-384, 1999.

[9] F. D. Day-Lewis, K. Singha, and A. M. Binley, "Applying petrophysical models to radar travel time and electrical resistivity tomograms: resolution-dependent limitations," *Journal of Geophysical Research: Solid Earth*, vol. 110, no. B8, 2005.

[10] Q. Lei, J.-P. Latham, and J. Xiang, "Implementation of an empirical joint constitutive model into finite-discrete element analysis of the geomechanical behaviour of fractured rocks," *Rock Mechanics and Rock Engineering*, vol. 49, no. 12, pp. 4799-4816, 2016.

[11] Q. Lei, J.-P. Latham, J. Xiang, C.-F. Tsang, P. Lang, and L. Guo, "Effects of geomechanical changes on the validity of a discrete fracture network representation of a realistic two-dimensional fractured rock," *International Journal of Rock Mechanics and Mining Sciences*, vol. 70, pp. 507-523, 2014.

[12] X. Zhang and D. J. Sanderson, "Numerical study of critical behaviour of deformation and permeability of fractured rock masses," *Marine and Petroleum Geology*, vol. 15, no. 6, pp. 535-548, 1998.

[13] S. Follin, L. Hartley, I. Rhén et al., "A methodology to constrain the parameters of a hydrogeological discrete fracture network model for sparsely fractured crystalline rock, exemplified by data from the proposed high-level nuclear waste repository site at Forsmark, Sweden," *Hydrogeology Journal*, vol. 22, no. 2, pp. 313-331, 2014.

[14] Q. Lei, J.-P. Latham, J. Xiang, and C.-F. Tsang, "Polyaxial stress-induced variable aperture model for persistent 3D fracture networks," *Geomechanics for Energy and the Environment*, vol. 1, pp. 34-47, 2015.

[15] C.-F. Tsang and I. Neretnieks, "Flow channeling in heterogeneous fractured rocks," *Reviews of Geophysics*, vol. 36, no. 2, pp. 275-298, 1998.

[16] H. Abelin, L. Birgersson, H. Widen, T. Ågren, L. Moreno, and I. Neretnieks, "Channeling experiment," Swedish Nuclear Fuel and Waste Management Co., 1990.

[17] H. Abelin, I. Neretnieks, S. Tunbrant, and L. Moreno, "Final report of the migration in a single fracture-experimental results and evaluation," Swedish Nuclear Fuel and Waste Management Co., 1985.

[18] P. A. Hsieh, S. P. Neuman, G. K. Stiles, and E. S. Simpson, "Field determination of the three-dimensional hydraulic conductivity tensor of anisotropic media: 2. Methodology and

application to fractured rocks," *Water Resources Research*, vol. 21, no. 11, pp. 1667-1676, 1985.

[19] H. Abelin, L. Birgersson, T. AAgren, and I. Neretnieks, "A channeling experiment to study flow and transport in natural fractures," in *MRS Online Proceedings Library Archive*, 1988.

[20] H. Abelin, L. Birgersson, J. Gidlund et al., "3--D migration experiment in sparsely fractured crystalline rock," *MRS Proceedings*, vol. 112, p. 199, 1987.

[21] Y. W. Tsang and C. F. Tsang, "Channel model of flow through fractured media," *Water Resources Research*, vol. 23, no. 3, pp. 467-479, 1987.

[22] Y. W. Tsang and C. F. Tsang, "Flow channeling in a single fracture as a two-dimensional strongly heterogeneous permeable medium," *Water Resources Research*, vol. 25, no. 9, pp. 2076-2080, 1989.

[23] A. Baghbanan and L. Jing, "Stress effects on permeability in a fractured rock mass with correlated fracture length and aperture," *International Journal of Rock Mechanics and Mining Sciences*, vol. 45, no. 8, pp. 1320-1334, 2008.

[24] W. Dang, W. Wu, H. Konietzky, and J. Qian, "Effect of shear-induced aperture evolution on fluid flow in rock fractures," *Computers and Geotechnics*, vol. 114, article 103152, 2019.

[25] P. K. Kang, Q. Lei, M. Dentz, and R. Juanes, "Stress-induced anomalous transport in natural fracture networks," *Water Resources Research*, vol. 55, no. 5, pp. 4163-4185, 2019.

[26] J.-P. Latham, J. Xiang, M. Belayneh, H. M. Nick, C.-F. Tsang, and M. J. Blunt, "Modelling stress-dependent permeability in fractured rock including effects of propagating and bending fractures," *International Journal of Rock Mechanics and Mining Sciences*, vol. 57, pp. 100-112, 2013.

[27] Q. Lei, J.-P. Latham, and C.-F. Tsang, "The use of discrete fracture networks for modelling coupled geomechanical and hydrological behaviour of fractured rocks," *Computers and Geotechnics*, vol. 85, pp. 151-176, 2017.

[28] K.-B. Min, J. Rutqvist, C.-F. Tsang, and L. Jing, "Stress-dependent permeability of fractured rock masses: a numerical study," *International Journal of Rock Mechanics and Mining Sciences*, vol. 41, no. 7, pp. 1191-1210, 2004.

[29] J. Rutqvist and O. Stephansson, "The role of hydromechanical coupling in fractured rock engineering," *Hydrogeology Journal*, vol. 11, no. 1, pp. 7-40, 2003.

[30] D. J. Sanderson and X. Zhang, "Critical stress localization of flow associated with deformation of well-fractured rock masses, with implications for mineral deposits," *The Geological Society*, vol. 155, no. 1, pp. 69-81, 1999.

[31] Q. Lei and X. Wang, "Tectonic interpretation of the connectivity of a multiscale fracture system in limestone," *Geophysical Research Letters*, vol. 43, no. 4, pp. 1551-1558, 2016.

[32] C. E. Renshaw, "Mechanical controls on the spatial density of opening-mode fracture networks," *Geology*, vol. 25, no. 10, pp. 923-926, 1997.

[33] N. E. Odling, "Scaling and connectivity of joint systems in sandstones from western Norway," *Journal of Structural Geology*, vol. 19, no. 10, pp. 1257-1271, 1997.

[34] O. Bour, "A statistical scaling model for fracture network geometry, with validation on a multiscale mapping of a joint network (Hornelen Basin, Norway)," *Journal of Geophysical Research*, vol. 107, no. B6, article 2113, 2002.
[35] Q. Lei and K. Gao, "Correlation between fracture network properties and stress variability in geological media," *Geophysical Research Letters*, vol. 45, no. 9, pp. 3994-4006, 2018.

[36] A. Roy, E. Perfect, W. M. Dunne, N. Odling, and J.-W. Kim, "Lacunarity analysis of fracture networks: evidence for scale-dependent clustering," *Journal of Structural Geology*, vol. 32, no. 10, pp. 1444-1449, 2010.

[37] N. E. Odling, "The scaling of hydraulic conductivity in rock fracture zones," *Geophysical Research Letters*, vol. 28, no. 15, pp. 3019-3022, 2001.

[38] A. Munjiza, *The Combined Finite-Discrete Element Method: Munjiza/Discrete Element Method*, John Wiley & Sons, Ltd, Chichester, UK, 2004.

[39] R. Goodman, *Methods of Geological Engineering*, St Paul. West Publishing Company, 1976.

[40] S. Saeb and B. Amadei, "Modelling rock joints under shear and normal loading," *International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts*, vol. 29, no. 3, pp. 267-278, 1992.

[41] N. Barton and V. Choubey, "The shear strength of rock joints in theory and practice," *Rock Mechanics*, vol. 10, no. 1-2, pp. 1-54, 1977.

[42] P. A. Witherspoon, J. S. Y. Wang, K. Iwai, and J. E. Gale, "Validity of cubic law for fluid flow in a deformable rock fracture," *Water Resources Research*, vol. 16, no. 6, pp. 1016-1024, 1980.

[43] P. S. Lang, A. Paluszny, and R. W. Zimmerman, "Permeability tensor of three-dimensional fractured porous rock and a comparison to trace map predictions," *Journal of Geophysical Research: Solid Earth*, vol. 119, no. 8, pp. 6288-6307, 2014.

[44] P. Renard and G. de Marsily, "Calculating equivalent permeability: a review," *Advances in Water Resources*, vol. 20, no. 5-6, pp. 253-278, 1997.

[45] C. Bruderer-Weng, P. Cowie, Y. Bernabé, and I. Main, "Relating flow channelling to tracer dispersion in heterogeneous networks," *Advances in Water Resources*, vol. 27, no. 8, pp. 843-855, 2004.

[46] J. Maillot, P. Davy, R. Le Goc, C. Darcel, and J. R. de Dreuzy, "Connectivity, permeability, and channeling in randomly distributed and kinematically defined discrete fracture network models," *Water Resources Research*, vol. 52, no. 11, pp. 8526-8545, 2016.

[47] A. Sornette, P. Davy, and D. Sornette, "Fault growth in brittle-ductile experiments and the mechanics of continental collisions," *Journal of Geophysical Research*, vol. 98, no. B7, pp. 12111-12139, 1993.

[48] X. Wang, Q. Lei, L. Lonergan, H. Jourde, O. Gosselin, and J. Cosgrove, "Heterogeneous fluid flow in fractured layered carbonates and its implication for generation of incipient karst," *Advances in Water Resources*, vol. 107, pp. 502-516, 2017.

[49] E. Bonnet, O. Bour, N. E. Odling et al., "Scaling of fracture systems in geological media," *Reviews of Geophysics*, vol. 39, no. 3, pp. 347-383, 2001.

[50] P. A. Cowie, D. Sornette, and C. Vanneste, "Multifractal scaling properties of a growing fault population," *Geophysical Journal International*, vol. 122, no. 2, pp. 457-469, 1995.

[51] R. J. Bell and P. Dean, "Atomic vibrations in vitreous silica," *Discussions of the Faraday Society*, vol. 50, p. 55, 1970.

[52] J. T. Edwards and D. J. Thouless, "Numerical studies of localization in disordered systems," *Journal of Physics C: Solid State Physics*, vol. 5, no. 8, pp. 807-820, 1972.

[53] W. C. Belfield, "Multifractal characteristics of natural fracture apertures," *Geophysical Research Letters*, vol. 21, no. 24, pp. 2641-2644, 1994.

[54] K. Bisdom, G. Bertotti, and H. M. Nick, "The impact of different aperture distribution models and critical stress criteria on equivalent permeability in fractured rocks," *Journal of Geophysical Research: Solid Earth*, vol. 121, no. 5, pp. 4045-4063, 2016.

[55] P. Szymczak and A. J. C. Ladd, "The initial stages of cave formation: beyond the one-dimensional paradigm," *Earth and Planetary Science Letters*, vol. 301, no. 3-4, pp. 424-432, 2011.

[56] C. A. Barton, M. D. Zoback, and D. Moos, "Fluid flow along potentially active faults in crystalline rock," *Geology*, vol. 23, no. 8, p. 683, 1995.
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

Hindawi
Submit your manuscripts at
www.hindawi.com