Interactive roles of geometrical distribution and geomechanical deformation of fracture networks in fluid flow through fractured geological media

Qinghu Lei, Xiogung Wng, Ki-Bok Min, Jonny Rutqvist

a Department of Earth Sciences, ETH Zürich, Zürich, Switzerland

b Laboratoire HydroSciences Montpellier, Université de Montpellier, Montpellier, France

c Department of Energy Systems Engineering, Seoul National University, Seoul, Republic of Korea

d Energy Geosciences Division, Lawrence Berkeley National Laboratory, Berkeley, USA

ABSTRACT In this study, the combined effects of geometrical distribution and geomechanical deformation of fracture networks on fluid flow through fractured geological media are investigated numerically. We consider a finite-sized model domain in which the geometry of fracture systems follows a power-law length scaling. The geomechanical response of the fractured rock is simulated using a hybrid finitediscrete element model, which can capture the deformation of intact rocks, the interaction of matrix blocks, the displacement of discrete fractures and the propagation of new cracks. Under far-field stress loading,the locally variable stress distribution in the fractured rock leads to a stress-dependent variable aperture field controlled by compression-induced closure and shear-induced dilatancy of rough fractures. The equivalent permeability of the deformed fractured rock is calculated by solving for the fracture-matrix flow considering the cubic relationship between fracture aperture and flow rate at each local fracture segment. We report that the geometrical connectivity of fracture networks plays a critical role in the hydromechanical processes in fractured rocks.A well-connected fracture system under a high stress ratio condition exhibits intense frictional sliding and large fracture dilation/opening, leading to greater rock mass permeability. However, a disconnected fracture network accommodates much less fracture shearing and opening,and has much lower bulk permeability.We further propose an analytical solution for the relationship between the equivalent permeability of fractured rocks and the connectivity metric (i.e. percolation parameter) of fracture networks, which yields an excellent match to the numerical results. We infer that fluid flow through a well-connected system is governed by traversing channels (forming an “in parallel” architecture) and thus equivalent permeability is sensitive to stress loading(due to stress-dependent fracture permeability),whilst fluid flow through a disconnected system is more ruled by matrix (linking isolated clusters “in series”) and has much less stress dependency.

Keywords:Fracture Stress Aperture Connectivity Permeability

1. Introduction

Natural fractures often form complex networks in the Earth’s crust and serve as important pathways for fluid migration in subsurface geological media (Tsang and Neretnieks, 1998). The geometrical distribution of these natural discontinuities controls the global connectivity of passage systems (Bonnet et al., 2001),while the geomechanical deformation of contacting fracture surfaces under in situ stress loading dominates the local transmissivity of individual channels(Rutqvist and Stephansson,2003).Thus,the understanding of fluid flow in fractured geological formations requires a comprehensive characterisation of both the geometrical properties and geomechanical responses of the embedded natural fracture networks (Zimmerman and Main, 2004; Lei et al., 2017a).

Fracture patterns are spatially organised by mechanical interactions that emerge during their growth (Pollard and Aydin,1988), which creates a hierarchical geometry exhibiting longrange correlations from macroscale frameworks to microscale fabrics (Barton, 1995). Such geometrical complexities are often depicted in terms of fracture size, density, orientation and aperture as well as connectivity(Bonnet et al.,2001).By mimicking the distributions of these properties following statistical/stochastic principles, synthetic discrete fracture networks are often generated to represent natural discontinuity systems (Dershowitz and Einstein,1988), based on which subsurface fluid flow properties and processes can be studied (Long et al.,1982;Long and Billaux,1987;Berkowitz,2002;Min et al.,2004a;Wang et al.,2017).Many theoretical and numerical investigations have suggested that the hydraulic connectivity of fracture networks, embedded in impervious or low permeability rocks, crucially controls the bulk permeability and fluid movements (Hestir and Long, 1990;Berkowitz, 1995; Bour and Davy, 1997, 1998; Renshaw, 1999; de Dreuzy et al., 2001a; Liu et al., 2016). In addition, the non-trivial impacts of variable fracture apertures on fluid flow have been analysed via making ad hoc hypotheses for their statistics (de Dreuzy et al., 2001b; Baghbanan and Jing, 2007; Klimczak et al.,2010), whereas the sensitive nature of fracture openings in response to geomechanical processes was often omitted.

The presence of natural fractures can induce strong local stress perturbations in geological media subjected to far-field stress loading (Lei and Gao, 2018, 2019). The locally varying stress distribution leads to variable normal and/or shear forces across fracture walls of widely-ranging sizes and orientations,producing a variety of fracture responses such as closure, sliding, dilatancy and propagation (Min et al.,2004b;Latham et al.,2013;Lei et al.,2016). Since the conductivity of fractures is critically dependent on the third power of fracture apertures (Witherspoon et al.,1980), these geomechanical processes accommodated in fractures can considerably affect the bulk hydrological properties of fractured porous media, especially if the matrix is much less permeable than fractures (Rutqvist, 2015). Overburden-induced confinement of fractured rocks tends to reduce fracture apertures and suppress fluid flow (Bandis et al., 1983; Barton et al.,1985), leading to the general trend of decreased rock mass permeability with increased formation depth (Rutqvist and Stephansson, 2003). On the other hand, large differential stresses imposed within the critically stressed crust can promote sliding and dilatancy along rough fractures, resulting in flow localisation and permeability enhancement, as revealed by many numerical results(Sanderson and Zhang,1999;Min et al.,2004b;Baghbanan and Jing, 2008; Latham et al., 2013; Lei et al., 2014,2015; 2017b; Lang et al., 2018; Kang et al., 2019; Jiang et al.,2019). However, these previous modelling studies mainly focused on sophisticated geomechanical processes while placing little emphasis on the role of geometrical properties, and therefore the analysis was usually only based on a specifically selected or generated fracture network.

In reviewing the extensive studies in the literature about hydromechanical processes in fractured rocks (Rutqvist and Stephansson, 2003; Lei et al., 2017a), we identify that little effort has been devoted to understanding the mutual effects of fracture network geometry and geomechanics on subsurface fluid flow as well as distinguishing the relative importance of each process.Thus, in this paper, we aim to use the state-of-the-art numerical simulation to gain insights into how these two aspects interactively affect the hydrological properties of fractured rocks.The remainder of the paper is organised as follows. Section 2 describes a set of approaches for modelling fracture network geometry, solid skeleton deformation, rough fracture behaviour and fluid flow field.Section 3 presents the model setup and boundary conditions for numerical experiments, with the numerical results further elucidated in Section 4.Finally,a discussion is given and conclusions are drawn in Section 5.

2. Methodology

2.1. Fracture network model

Natural fracture systems often exhibit a broad range of fracture lengths that can be described by a power-law statistical model with a density function given as (Bonnet et al., 2001; Lei and Wang,2016):

wherelis the fracture length,Lis the domain size,ais the powerlaw length exponent,Dis the fractal dimension,and α is the density term. The only intrinsic characteristic length scales in this model are the smallest and largest fracture lengths, i.e.lminandlmax,respectively. In numerical simulations,Lis the scale of the modelling domain, which usually meetslmin＜＜L＜＜lmax. Extensive outcrop data suggest thatDvalue generally varies between 1.5 and 2, andafalls between 1.5 and 3.5 (Bonnet et al., 2001).

The fracture intensity γ (also known as the mass density), i.e.total length of fractures per unit area,is related to the first moment of the density distribution of fracture lengths as

wherel′denotes the fracture length included in the domain of an areaAL=L2.The percolation parameterpas a connectivity metric of fracture networks is given by(Bour and Davy,1997):

The higher thepis, the more connected the system is. The network is statistically connected ifpis greater than the percolation thresholdpc,whose value is in general scale-invariant and within a narrow range between 5 and 7,i.e.pc=6±1 for two-dimensional(2D) random fracture networks withD= 2 (Bour and Davy,1997).

2.2. Geomechanical model

The geomechanical model is based on a hybrid finite-discrete element method (FDEM) (Munjiza, 2004), which can realistically capture the stress/strain evolution in intact rocks, interaction between matrix blocks, deformation of pre-existing fractures, and propagation of new cracks (Lei et al., 2017a). The FDEM model represents 2D fractured rock using an unstructured, fullydiscontinuous mesh of three-node triangular finite elements,which are linked by four-node joint elements(see Fig.1).There are two types of joint elements: unbroken joint elements inside the matrix and broken joint elements along fractures (Lei et al., 2016).The joint elements (either broken or unbroken) are created and embedded between the edges of triangular element pairs before the numerical simulation and no remeshing is performed during later computation.

The motion of finite elements is governed by the forces acting on the elemental nodes (Munjiza, 2004):

Fig. 1. Representation of a 2D fractured rock using an unstructured, fullydiscontinuous mesh of three-node triangular finite elements linked by four-node joint elements.

whereMis the lumped nodal mass matrix;xis the vector of nodal displacements;fintis the internal nodal force induced by the deformation of triangular elements; andfextis the external nodal force including external loadflcontributed by boundary conditions and body forces, cohesive bonding forcefbcaused by the deformation of unbroken joint elements, and contact forcefcgenerated by the contact interaction via broken joint elements. The deformation of intact rocks is captured by linear-elastic constant-strain finite elements with the continuity constrained by the bonding forces of unbroken joint elements (Munjiza et al., 1999). The interaction of discrete matrix bodies is calculated based on the penetration of finite elements via broken joint elements (Munjiza and Andrews, 2000). The elasto-plastic fracturing in formation rocks is modelled by a smeared crack (i.e. cohesive zone) method that can capture the nonlinear stress-strain behaviour of the plastic zone ahead of each fracture tip(Munjiza et al.,1999).The equations of motion of the FDEM model are solved through an explicit time integration scheme based on the forward Euler method.

2.3. Fracture constitutive model

The closure of rock fractures under compression is calculated based on a hyperbolic relation (Bandis et al.,1983) (Fig. 2a):

wherevnis the normal closure, σnis the effective normal compressive stress derived from the Cauchy stress tensor of adjacent finite elements,kn0is the initial normal stiffness,andvmis the maximum allowable closure.

The shear deformation of rock fractures is calculated based on an elasto-plastic constitutive model with strain-softening(Goodman,1976; Saeb and Amadei,1992) (Fig. 2b). In the elastic phase, the shear stress τ increases linearly with the shear displacementu,and the slope of the stress-displacement curve is given by the shear stiffnessks. During this stage, the opposing fracture walls ride over each other’s asperities, resulting in dilational displacement in the normal direction (Fig. 2c). The peak shear stress τpis eventually reached when the displacement arrives at the peak shear displacementup, beyond which the asperities begin to shear off and irreversible damage on the surfaces starts to occur. If the fracture continues to slide, the shear stress decreases linearly to the residual shear stress τr,during which the asperities are crushed and sheared off and the dilation continues.Finally,when the displacement exceeds the residual displacementur, the shear stress remains constant (i.e. τ = τr), and no further dilation develops. The dependency of the shear behaviour of fractures on normal stress loading is described using a constant displacement model characterised with fixedupandurvalues(Goodman,1976).

The peak shear stress τpis given by (Ladanyi and Archambault,1969; Saeb and Amadei,1992):

whereasis the proportion of total fracture area sheared through asperities; φiis the dilation angle;cis the shear strength of the asperity(i.e.cohesion of the intact rock);and φbis the basic friction angle which,for unweathered conditions,can be substituted using the residual friction angle φr(Barton and Choubey,1977).If σndoes not exceed the uniaxial compressive strength of the intact rock σu,the values ofasand φiare respectively given as (Ladanyi and Archambault,1969):

where φi0is the initial dilation angle when σn=0;andm1andm2are the empirical parameters with suggested values of 1.5 and 4,respectively. The residual shear stress τris given as (Barton and Choubey,1977):

Fig. 2. Fracture deformation model: (a) Normal closure vn as a function of effective normal stress σn; (b) Shear stress τ as a function of shear displacement u; and (c) Dilational displacement vs as a function of shear displacement u.

The dilational displacementvsunder a constant normal stress condition is related to the shear displacement in an incremental form as (Saeb and Amadei,1992):

The fracture aperturehunder coupled normal and shear loadings is thus given by(Lei et al., 2016):

whereh0is the initial aperture,andwis the separation of opposing fracture walls if the fracture is under tension. The local fracture permeability is then calculated based on the cubic law ash2/12(Witherspoon et al.,1980).

2.4. Fluid flow model

Fluid flow through fractured rock with multiple intersecting fractures and permeable matrix is further solved. Single-phase steady-state flow of incompressible fluid with constant viscosity through porous media,in absence of sources and sinks,is governed by the continuity equation and Darcy’s law, which is reduced to a Laplace equation as

wherekis the intrinsic and isotropic permeability of the porous media with local variability permitted, andPis the fluid pressure calculated at the nodes of unstructured finite element grids. The element-wise constant barycentric velocity is resolved based on the pressure gradient vector field by applying Darcy’s law given by

whereueis the vector field of element-wise constant velocities;Peis the local element pressure field;μ is the dynamic fluid viscosity;andkeis the local permeability of a matrix volumetric element with an assumed constant value or a lower dimensional fracture element having a stress-dependent value,i.e.h2/12.By applying a prescribed macroscopic pressure differential on each pair of opposite boundary surfaces with no-flow conditions on the remaining ones parallel to the flow direction,the pressure diffusion is computed for all the fracture and matrix elements of the entire domain. The equivalent permeability of the fractured media is then derived using element volume weighted averaging of pressure gradients and fluxes:

3. Model setup

In this study, we generate a series of 2D synthetic fracture networks in a square domain of sizeL=10 m(Fig.3).The location and orientation of fractures are assumed purely random, i.e.nominally homogeneous (i.e.D= 2) and isotropic. The fracture lengths follow the power-law scaling, with the bounds given bylmin=L/50 = 0.2 m andlmax= 50L= 500 m. We explore five different length exponent cases, i.e.a= 1.5, 2, 2.5, 3 and 3.5, and four mean fracture intensity scenarios, i.e. γ = 1.25 m-1, 2.5 m-1,3.75 m-1and 5 m-1. For each combination ofaand γ, ten realizations are generated,with theirpvalues also derived.As shown in Fig. 3, whena≤2, the system is dominated by domain-sized,very large fractures; whena≥3, the system mainly consists of small fractures;when 2 ＜a＜ 3,both large and small fractures tend to be important.It can also be noted that an increasedaleads to a reduction in the geometrical connectivity, i.e. decreasedp. The fracture networks generated here cover the scenarios of below(p＜ 5), around (5 ≤p≤ 7) and beyond (p> 7) the percolation threshold, representing the disconnected, transitional and connected regimes, respectively.

The assumed material properties of the fractured rocks are given in Table 1, based on the typical ranges for crystalline rocks in the literature (Lama and Vutukuri, 1978; Bandis et al., 1983; Zoback,2007). The energy release rates are estimated based on empirical correlations (Zhang, 2002; Jin et al., 2011). The problem domain containing distributed fractures is discretised using an unstructured mesh with an average element size of 0.05 m. The penalty term and damping coefficient are chosen to be 500 GPa and 2 × 105kg/(m s), respectively, based on the recommendations inthe literature (Munjiza, 2004; Mahabadi, 2012). Effective far-field stresses are loaded orthogonally to the model (Fig. 4a), and we consider three different scenarios: (i)Sx= 5 MPa,Sy= 5 MPa, (ii)Sx=10 MPa,Sy=5 MPa,and(iii)Sx=15 MPa,Sy=5 MPa,such that the effective far-field stress ratioSx/Syis 1, 2, and 3, respectively.Single-phase steady-state fluid flow through the deformed fractured rock having a stress-dependent aperture field is further modelled by imposing the classical permeameter boundary condition (Fig. 4b): two opposite boundary surfaces of the model domain have a fixed pressure drop (i.e. 10 kPa), while the two orthogonal boundaries parallel to the flow direction are impervious.Matrix permeabilitykmis assumed to be 1×10-18m2,which gives a high fracture-matrix permeability contrast so that the flow is predominated by fractures (Matth?i and Belayneh, 2004).

Table 1Material properties of the fractured rocks (Lama and Vutukuri,1978; Bandis et al.,1983; Zoback, 2007).

The constructed models are numerically solved using the methods described in Section 2. It is worth pointing out that the advantages of our approaches, compared to many other previous studies/methods(e.g.Min et al.,2004b;Baghbanan and Jing,2008),include simulation of crack growth in intact rocks and consideration of fluid flow in permeable matrix.Newly-formed cracks under stress loadings may cause coalescence of initially isolated fractures/clusters, capture what is essential for modelling the hydromechanical behaviour of fracture networks around the percolation threshold. The consideration of fluid flow in permeable matrix is important for studying fracture networks below the percolation threshold, where fluid cannot migrate purely through fractures.

4. Results

4.1. Geomechanical responses

We analyse the geomechanical responses of the fractured rocks associated with a range of combinations of length exponentaand fracture intensity γ values under different far-field stress conditions.The numerical results for the case of γ=2.5 m-1are shown in Fig.5,while the results for other cases of γ=1.25 m-1,3.75 m-1and 5 m-1are given in Figs.A1,A2 and A3,respectively(see Appendix).Below, we take Fig. 5 as an example to elucidate the key geomechanical processes in fracture networks and similar phenomena or trends can also be observed in Figs. A1-A3.

Fig. 5a shows the distributions of local maximum principal stresses in the fractured rocks with γ = 2.5 m-1and differentavalues. When the far-field stress condition is isotropic, i.e.Sx=Sy= 5 MPa, the local stress distribution is very uniform. As the far-field stress ratioSx/Syincreases,stress fluctuations emerge in the system, especially whenSx/Syreaches 3. It can be seen that the high stress bands tend to align with the orientation of the applied far-field maximum principal stress but are significantly distorted due to the presence of pre-existing discontinuities. It seems that with the increase ofa(i.e. the system is more dominated by small fractures), the stress patterns become less heterogeneous.

Fig. 5b shows the distribution of shear displacements in the stressed fracture networks. When the far-field stresses are isotropically loaded, almost no shear displacement occurs in the system, irrespective of the geometrical distribution of fracture populations. AsSx/Syvalue increases to 2, noticeable shear displacements are accommodated along some of the large fractures that are oriented to favour frictional sliding. WhenSx/Syequals 3,preferentially-oriented,large fractures in the networks ofa≤2 are highly reactivated for shear slip, whereas the remaining relatively smaller fractures experience much less shearing. In the fracture networks ofa≥3,which are dominated by small cracks,frictional sliding is strongly suppressed, although some intermediate fractures might exhibit slight shear.In the fracture networks of 2 ＜a＜3, which consist of both large and small fractures, some large structures tend to be moderately sheared while the small cracks are mostly restrained for any sliding.

Fig.4. Model setup for(a)geomechanical simulation and(b)fluid flow simulation.Sx and Sy denote the effective far-field stresses loaded orthogonally to the model.P1 and P2 gives the macroscopic pressure difference imposed across the domain.

In Fig.5c, we show the distribution of fracture apertures in the fractured networks under the combined effects of compressioninduced closure and shear-induced dilatancy. WhenSx/Sy= 1, all fractures are evenly compressed and exhibit an aperture much lower than the initial value,i.e.0.1 mm.AsSx/Syincreases to 2,some traversing fractures in the networks ofa≤2 exhibit large apertures because of shear-induced dilation along these dominant structures that are preferentially-oriented for sliding. However, some small fractures sub-parallel toSyseem to be more closed due to the increasedSx, compared to those under isotropic compression.WhenSx/Syreaches 3,in the networks ofa≤2,drastic enlargement of fracture apertures occurs along those highly-sheared large discontinuities due to the dilational behaviour of dislocated rough fracture walls. However, the fracture networks ofa≥3 exhibit much less fracture opening, although some intermediate fractures still experience small amount of shear-induced aperture increase.In the fracture networks of 2 ＜a＜ 3,large fractures are associated with shear-induced wider apertures while small fractures only exhibit closure manner.

We further calculate the length-averaged mean shear displacementand mean fracture apertureof each fracture network and examine their variation as a function of various geometrical properties, i.e.a, γ andp, under different far-field stress loading conditions (Fig. 6). WhenSx/Sy= 1,is almost zero and insensitive to the change ofa, γ orp(Fig. 6a, c and e),whileis also independent ofa, γ orp(Fig. 6b, d and f) and exhibits a value much lower than the initial aperture ofh0=0.1 mm.AsSx/Syincreases, more frictional sliding is accommodated in the system due to the enhanced differential stress load. WhenSx/Sy> 1,increases significantly with the decrease ofa(Fig. 6a) or increase ofp(Fig. 6e), suggesting that the fracture networks associated with larger fractures or better connectivity tend to accommodate more frictional sliding driven by deviatoric stress loading. Similarly,increases considerably with the decrease ofa(Fig.6b)or increase ofp(Fig.6f),as a result of fracture dislocation and dilation under differential stresses. Especially, in the wellconnected fracture systems (p>pc) subjected to high stress ratio loading(Sx/Sy=3),the resulting aperturecan become even larger than the initial apertureh0(Fig.6f)due to the combined effects of geometrical properties and geomechanical processes. In addition,it is noticed that bothseem to be almost uncorrelated with γ (Fig. 6c and d). A further interpretation and discussion can be found in in Section 5.

4.2. Hydrological properties

We derive the hydrological properties of the deformed fractured porous rocks from the single-phase steady-state flow simulation.Fig.7 shows the distribution of local fluid flow velocity in different fracture systems subjected to the far-field stress condition ofSx= 15 MPa andSy= 5 MPa and with the macroscopic pressure drop imposed along thex-(Fig. 7a) ory-direction (Fig. 7b). The numerical results for other far-field stress conditions are given in the Appendix(see Figs.A4 and A5).In the networks that are below the percolation threshold (p＜pc), the flow velocity is extremely low,since fluid has to migrate via intact rocks of low permeability that bridge disconnected fractures. On the other hand, the wellconnected networks (p>pc) accommodate significantly high flow velocities along through-going discontinuities, whereas the“background” small cracks with “dead-ends” provide much lower velocities. If the system is around the percolation threshold(p≈pc),e.g.the fracture network ofa=2.5 and γ =5 m-1,only one or two globally-connected pathways exist(s) for fast fluid migration,whereas other locally-connected clusters permit much slower flow. In such a critically-connected system (p≈pc), multiple clusters bounded by large fractures seem to exhibit distinct regimes of velocity magnitudes, implying the dominant roles of large discontinuities on the flow field.

Fig. 8 shows the variation of the equivalent permeabilitykeq(eitherkxxorkyy) of the fractured rocks as a function of various geometrical properties, i.e.a, γ andp, under different far-field stress loading conditions. It can be seen that, with the increase ofa,the permeability generally decreases(Fig.8a and b),because the system becomes more controlled by small fractures, which geometrically have less probability to form connected pathways(i.e. poorer connectivity) and geomechanically have less opportunity to experience frictional sliding (i.e. smaller dilation). Due to these geometrical and geomechanical effects,keqis very sensitive to stress loading ifais small (e.g.a= 1.5), whilst almost independent of the far-field stress state ifais large(e.g.a=3.5).The rock mass permeability seems to increase with γ, but the trend exhibits large uncertainties(Fig.8c and d).On the contrary,the correlation betweenkeqandpis very significant (Fig. 8e and f). If the fracture network is disconnected (p＜pc),keqis very small, due to the flow-restriction caused by low-permeability rocks that isolate discrete fractures or fracture clusters from connecting with each other across the entire domain.The stressdependent fracture deformational behaviour therefore exerts very minor impacts onkeq. However, as the fracture network becomes gradually connected (p≥pc), an abrupt increase inkeqoccurs, since the interconnected backbones of high transmissivity(compared to matrix rocks)start to play the controlling role in carrying fluid migration, which further leads to the variation ofkeqin response to the change of boundary stress loading as a result of the important geomechanical processes illustrated in Fig. 5, A1 and A2.

Fig. 6. The mean shear displacement (left panel) and mean fracture aperture (right panel) of the fractured network as a function of different geometrical properties,i.e. (a, b)length exponent a, (c, d) fracture intensity γ, and (e, f) percolation parameter p, under different far-field stress conditions. The error bars denote ±1 standard deviation of the ten realisations. The zone bounded by dotted lines in (e, f) corresponds to the percolation threshold pc = 6 ± 1.

We further derive analytical solutions to investigate the relationship betweenkeqandp.We characterise their correlation based on the three different regimes: disconnected (p＜pc), transitional(p≈pc) and connected (p>pc). If the system is connected, fluid flow is controlled by interconnected fractures that form the major pathways and the equivalent permeability of the percolated fracture system may be predicted using the percolation theory(Stauffer and Aharony,1985):

where β is a universal exponent and equals 1.1 for 2D systems(Hestir and Long,1990;Berkowitz and Balberg,1993),andkcis the equivalent permeability whenp-pcequals unity.kcis expected to be controlled by single fracture permeability,because the system is close to the percolation threshold and the connectivity is ruled by a limited number of“red links”(Davy et al.,2006).Thus,we postulate thatkc= λh3/(12L), where λ is the number of red links. We takepc=6,L=10 m and λ=3(papproachespcfrom above).Then,for the three different stress cases ofSx/Sy= 1, 2 and 3, we chooseh=0.05 mm,0.07 mm and 0.11 mm(Fig.6),respectively,and thus derivek0to be 3 × 10-15m2, 8.6 × 10-15m2and 2.5 × 10-14m2,respectively. The analytical solutions based on these parameters show an excellent match to the numerical results for the connected regime(see the dashed lines in the region ofp>pcin Fig.8e and f).For the disconnected regime (p＜pc),keqis related to lacunarity which characterises the gaps (i.e. rock bridges) between isolated fractures or clusters. Thus,keqis insensitive to stress loading that mainly affects fractures. We suspect the equivalent permeability follows

Fig.7. Distribution of flow velocity in the fracture networks associated with various length exponents a and fracture intensities γ under the far-field stress condition of Sx=15 MPa and Sy = 5 MPa. The macroscopic pressure drop is imposed along either the (a) x or (b) y direction.

wherekeq=kmifp= 0, and ζ is an exponent. The value of ζ is derived as follows: ifpc-pis unity,which gives ζ=log10(kc/km)/log10pc.By takingpc=6 and λ=1(papproachespcfrom below), we also obtain an excellent match to the simulation data (see the dashed lines in the region ofp＜pcin Fig. 8e and f).The rationale of determining the number of red links λ is given as follows.Whenpapproachespcfrom below,i.e.papproaches 5,one red link is expected and therefore λ=1.Whenpapproachespcfrom above,i.e.papproaches 7,two more traversing flow channels may be added and thus λ = 3 (see the definition of the percolation parameter in Eq. (3), which indicates that ifpincreases by 1, one more domain-sized percolating channel is added into the system).For the narrow transitional regime (p≈pc), we only use a simple interpolation from the solutions of the two neighbouring regimes(see the dashed lines in the region ofp≈pcin Fig. 8e and f). We have further tested the analytical formulation for a different matrix permeability value (km= 1 ×10-21m2) in the Appendix (Fig. A6),and a good match to numerical results is also obtained.

Fig.8. The equivalent permeability kxx(left panel)and kyy(right panel)of fractured rocks as a function of different geometrical properties,i.e.(a,b)length exponent a,(c,d)fracture intensity γ, and (e, f) percolation parameter p, under different far-field stress conditions. The error bars denote ±1 standard deviation of the ten realisations. The insets present permeability on the logarithmic scale. In (e, f), the zone bounded by dotted lines corresponds to the percolation threshold pc = 6 ± 1; and the dashed lines correspond to the analytical solutions, i.e. Eqs. (15) and (16).

5. Discussion and conclusions

In this paper, we presented a systematic investigation into the hydromechanical processes (e.g. stress distribution, fracture closure, shear dilation and fluid flow) in fracture networks of a broad range of geometrical distributions subjected to different stress loading conditions. We examined the variation of hydromechanical properties (e.g. shear displacement, fracture aperture and equivalent permeability) as a function of various geometrical properties (e.g. length exponenta, fracture intensity γ and percolation parameterp).We found that fracture intensity γ is not a good proxy for the hydromechanical behaviour of rock masses,because it poorly (or not at all) indicates the connectivity state of fracture networks (as manifested in Fig. 3 where the networks with the same γ value can have very different connectivity conditions).This observation is consistent with our field characterisation data at the Grimsel Test Site, Switzerland, where a very poor correlation was found between the interval transmissivity (derived from pulse tests) and fracture spacing (fracture intensity measured in one dimension based on borehole logs)(Brixel et al.,2020a,b).Thus,it is essential to characterise the length distribution of fracture populations in addition to the measurement of fracture intensity (or fracture spacing), as suggested by many previous studies (Bonnet et al., 2001; Davy et al., 2010). The combination of the information of fracture length and intensity permits the calculation of the percolation parameter, which is metric of the fracture network connectivity and serves as an excellent proxy to the hydromechanical behaviour of fractured rocks. This is consistent with past research findings in the literature about the connectivity control on the solid deformation(Zhang and Sanderson,1998;Harthong et al.,2012;Lei and Gao,2018)and fluid flow(Hestir and Long,1990;Bour and Davy, 1997; Renshaw, 1999; de Dreuzy et al., 2001a) of fractured geological media, while our work further advanced the understanding on the mutually-existing, interactively-acting impacts of geometrical connectivity and geomechanical condition on hydrological performance of the system.

Based on our numerical results in Section 4, it seems that the“connectedness” of the fracture system is a prerequisite for notable hydromechanical processes taking effect in fractured rocks. Thus, the system is almost insensitive to the stress loading conditions if it is disconnected.This is because when the fracture network is not percolated, fluid needs to migrate via isolated clusters and the matrix in between them. Hence, the system has an “in series” flow structure (i.e. isolated clusters are linked with each other by matrix), such that the bulk permeability is approximated by the harmonic mean of the permeabilities of different components (de Marsily,1986) and thus more ruled by the matrix (because it has a low permeability and dominates the harmonic mean). As the system gradually approaches the percolation threshold from below, the size of matrix blocks that form the gaps between isolated clusters is reduced and, therefore, the equivalent permeability increases although the system is still disconnected(Fig.8e and f).For the connected regime,the system has an “in parallel” flow structure consisting of multiple traversing channels, such that the bulk permeability is approximated by the arithmetic mean of the permeabilities of different components (de Marsily, 1986) and thus controlled by fractures(because they have high permeability and dominate the arithmetic mean). Since fracture aperture (and thus fracture permeability) is dependent on the stress state, the equivalent permeability of fractured rocks also shows a strong dependency on the far-field stress loading. However, it is worth pointing out that the shear-induced aperture enlargement under high stress ratio conditions tends to enlarge the bulk permeability by several times but still within an order of magnitude. We inferred that geomechanical processes tend to exert a secondary-order influence (compared to the first-order role of geometrical connectivity) causing stress-dependent variation of permeability within about one order of magnitude, for which similar permeability variation ranges were also reported in previous studies based on 2D fracture networks (Min et al., 2004b; Baghbanan and Jing,2008; Latham et al., 2013; Lei et al., 2014). However, more pronounced stress effects might be expected in three-dimensional(3D) fracture systems (Lei et al., 2015, 2017b), which may be of interest for further investigations based on computationally expensive 3D simulations. Further work may also be needed to explore more complicated scenarios of fracture networks,such as anisotropic discontinuity orientations (Lei et al., 2014), multiple domain sizes (Min et al., 2004a), variable initial apertures (Kang et al., 2019), fractal spatial organisations (Davy et al., 2006), and nonlinear mechanical-to-hydraulic aperture correlations(Renshaw, 1995) for testing the generality and sensitivity of the observed phenomena in the current paper.

To conclude, we observed that the geometrical connectivity of fracture networks plays a critical role in the hydromechanical processes in fractured rocks. A well-connected fracture system under a high stress ratio exhibits intense frictional sliding and considerable fracture opening, eventually leading to fast fluid migration and large bulk permeability.Such a connected network is more suppressed for shearing activities under a more isotropic compression, exhibiting a lower permeability compared to that under a high stress ratio loading.A disconnected fracture network is composed of multiple clusters isolated from each other by geomechanically stiffer and hydraulically less permeable matrix rocks,producing mainly fracture closure behaviour and slow fluid flow,with low equivalent permeability insensitive to the far-field stress state. We derived an analytical solution for the relationship between the equivalent permeability of fractured rocks and the percolation parameter of fracture networks, which showed an excellent match to the numerical results. We suggested that the flow through a well-connected system is governed by traversing fractures or clusters “in parallel” and thus the equivalent permeability is sensitive to stress loading(due to the stress dependency of fracture permeability), whilst the flow through a disconnected system is more ruled by matrix, which links isolated clusters “in series”, and therefore very insensitive to stress loading. The observation and analysis of interactively superimposed geometrical and geomechanical effects on hydrological properties of fractured geological media as presented in this paper have important implications for understanding heterogeneous subsurface fluid flow and upscaling rock mass permeability.

Declaration of Competing Interest

The authors wish to confirm that there are no known conflicts of interests associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgments

Qinghua Lei is grateful for the support from Swiss National Science Foundation (Grant No. IZLCZ0_189882).Xiaoguang Wang was partially funded by PRC-CNRS Joint Research Project (Grant No. 5181101856). Ki-Bok Min was partially supported by the Korea-EU Joint Research Support Program of the National Research Foundation of Korea through a grant funded by the Korean Government’s Ministry of Science, ICT and Future Planning (Grant No. NRF-2015K1A3A7A03074226) and Research Institute of Energy & Resources, Seoul National University. Support for Jonny Rutqvist was provided by the U.S. Department of Energy under contract No. DE-AC02-05CH11231 to the Lawrence Berkeley National Laboratory.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2019.12.014.