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

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 stressdependent 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 shearinduced 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.


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][2][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][6][7][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][11][12]. The aperture heterogeneity may lead to a highly channelized flow pattern [7,[13][14][15]. Extensive field observations have shown that the contribution of individual fractures in the network to the total flow is not equal [16][17][18]. Instead, the majority of the fluid flow is restricted to only a small number of fractures [19][20][21][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][24][25][26][27][28][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.

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

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): where v n is the normal closure, σ 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 τ p is given by [40]: 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 ϕ 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 a s and ϕ i are, respectively, given as [40]: where ϕ i0 is the initial dilation angle when σ 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 τ r is given as [41]: The dilational displacement v s is related to the shear displacement u in an incremental form as [40]: The fracture aperture h under coupled normal and shear loadings is thus given by [10]: where h 0 is the initial aperture, and w is the separation of opposing fracture walls if the fracture is under tension.

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, 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: while u e is the velocity vector field, p e is the local pressure at element nodes, μ 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.

Geofluids
The network equivalent permeability is obtained using Darcy's law: 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 in and p 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]: 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]: 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.

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°w ith an increment of 10° (Figure 2(a)). We consider five 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 × 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

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 θ = 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 compressioninduced 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 farfield 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

Geofluids
( 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 farfield 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  6 Geofluids 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 θ 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 θ 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 compressioninduced 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 farfield stress orientation becomes very complex under the combined effects of compression-induced closure, shearinduced dilation, and formation of new cracks. In general, k eff increases as the increase of differential stress so far but  7 Geofluids 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 farfield 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 θ 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 θ = 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   Geofluids permeability enhancement and flow structure alternation under high differential stresses (Figure 7), we show, in Figure 10, the derived k 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 9 Geofluids 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).

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  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 secondorder 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  11 Geofluids 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) Shearinduced dilation tends to be a second-order factor on altering permeability

Conflicts of Interest
The authors declare that they have no conflicts of interest.

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 θ = 90°in the cases of y-direction, whereas it occurs when θ = 30 and 150°in x-direction ( Figure S3 and Figure 6). Note that θ is the angle between maximum stress and x positive direction. The variation tendency of the quantitative flow indices (k eff /k 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 xand 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).