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.

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 [13].

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 [58]. 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 [1012]. The aperture heterogeneity may lead to a highly channelized flow pattern [7, 1315]. Extensive field observations have shown that the contribution of individual fractures in the network to the total flow is not equal [1618]. Instead, the majority of the fluid flow is restricted to only a small number of fractures [1922]. 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 [2329]. 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 difficulty to obtain, the sampled fracture patterns are limited to 2D. The fracture network covering an area of 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): where is the normal closure, is the effective normal compressive stress, is the initial normal stiffness, and 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 and values [39]. In this model, the peak shear stress is given by [40]: where is the proportion of total fracture area sheared through asperities, is the dilation angle, is the shear strength of the asperity (i.e., cohesion of the intact rock), and is the basic friction angle which is substituted using the residual friction angle [41]. If does not exceed the uniaxial compressive strength of the intact rock , the values of and are, respectively, given as [40]: where is the initial dilation angle when , and and are empirical parameters with suggested values of 1.5 and 4.0, respectively. The residual shear stress is given as [41]:

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

The fracture aperture under coupled normal and shear loadings is thus given by [10]: where is the initial aperture, and 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, where is the intrinsic permeability and 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 is the velocity vector field, is the local pressure at element nodes, is the dynamic fluid viscosity, and 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, 4548].

The network equivalent permeability is obtained using Darcy’s law: where is the volumetric flow rate through the inlet face of the fractured rock, is the length of the model domain, is the outlet cross-section area, and and are the hydraulic head at inlet and outlet, respectively.

The correlation dimension of the flow rate field is defined as [49, 50]: where is the proportion of flow within an elemental area to the cumulative amount of flow in the entire grid, and is the dimension of the grid elements. In this work, we use a 30 by 30 square grid in our calculation of , thus equals to 0.6 m.

The flow channeling density indicator , similar to the participation ratio in the physics literature [51, 52], is defined as [46, 47]: where is the fracture flow rate, is the fracture length, is the cumulative fracture length of the entire network. The inverse of defines the average spacing of main flow paths [46]. A small 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 different stress scenarios: (i)  MPa (reference case); (ii)  MPa; (iii)  MPa,  MPa; (iv)  MPa,  MPa; and (v)  MPa,  MPa. The initial aperture field is assumed constant as  mm, which is representative for nature fractures [5355].

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 is assumed constant in space with a value of  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 -direction; the results for -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 , 30°, 60°, 90°, 120°, and 150°. The probability density functions (PDFs) of these cases are given in Figure 4.

When  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  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  MPa and  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  MPa and  MPa, the shear-induced dilation of fracture apertures is more prevailing along some long fractures (Figure 3). When 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, . 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.,  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  MPa and  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  MPa and  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 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., 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  MPa and  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  MPa and  MPa, there is only a quantitative variation in the relative magnitude of fracture flow when is varied; however, when  MPa and  MPa or  MPa and  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 , the correlation dimension of flow rate , and the flow channeling density indicator , for various stress conditions. The equivalent permeability may quantify the bulk flow magnitude through the fractured rock, while and further indicate the flow distribution within the fracture network.

As shown in Figure 9(a), if the fractured rock is isotropically loaded ( or 5 MPa), 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 MPa and MPa), decreases by about several times due to compression-induced aperture closure (Figure 3). On the other hand, both and increase, indicating a less channelized flow distribution in the fractured rock. If the fractured rock is anisotropically stressed, varies with the far-field stress orientation. When MPa and MPa, the becomes larger than that of the isotropic stress loading cases of MPa and MPa but is still smaller than that of the zero-stress loading case (i.e., MPa and MPa) due to the dominant normal closure effect under confining stresses. The variation of with far-field stress orientation is quite small. A variation in and 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 MPa and MPa), the 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, 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, and also become smaller, indicating a higher level of flow channelization. The variation magnitude of and with stress becomes more prominent. The variation trends of and as the far-field stress field orientation appear similar (Figures 9(b) and 9(c)). As is further increased to 20 MPa, there is a significant overall increase in (Figure 9(a)). For all far-field stress orientations, the values are now larger than that of the fractured rock under the zero-stress loading. Both and 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 and values remain high and similar to those of the case of MPa and MPa. When , 130°, and 140°, an enhanced reduction in is observed (Figures 9(b) and 9(c)), while there is a slightly less decrease in . 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 permeability enhancement and flow structure alternation under high differential stresses (Figure 7), we show, in Figure 10, the derived , , and 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 and , 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).

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 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 ([email protected]).

Supplementary Materials

In this supplementary file, we provide flow simulation results in the -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 -direction (Figures S2, S3, S4, and S5) are similar to those in the -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 -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 in the cases of -direction, whereas it occurs when and 150° in -direction (Figure S3 and Figure 6). Note that is the angle between maximum stress and positive direction. The variation tendency of the quantitative flow indices (, , and ) under different stress conditions (Figures S4 and S5) is coincident with the cases in the -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 -direction flow. The difference between the results of - and -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)