Polyaxial stress-induced variable aperture model for persistent 3D fracture networks

Geomechanics for Energy and the Environment 1 (2015) 34–47

Contents lists available at ScienceDirect

Geomechanics for Energy and the Environment journal homepage: www.elsevier.com/locate/gete

Polyaxial stress-induced variable aperture model for persistent 3D fracture networks Qinghua Lei a,∗ , John-Paul Latham a , Jiansheng Xiang a , Chin-Fu Tsang b,c a

Department of Earth Science and Engineering, Imperial College London, London, UK

b

Department of Earth Sciences, Uppsala University, Uppsala, Sweden

c

Earth Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA, USA

highlights • • • • •

A stress-induced aperture model is developed for 3D persistent fracture networks. Aperture variations due to fracture-scale and network-scale effects are modelled. Local hydraulic apertures can vary greatly in a single persistent fracture. Large aperture channels are formed in critically stressed rock masses. Flow localisation occurs when the stress ratio exceeds the critical threshold.

article

info

Article history: Received 24 November 2014 Received in revised form 12 February 2015 Accepted 12 March 2015 Available online 20 March 2015 Keywords: FEMDEM Fracture apertures Stress Permeability Flow localisation

∗

Corresponding author. E-mail address: [email protected] (Q. Lei).

http://dx.doi.org/10.1016/j.gete.2015.03.003 2352-3808/© 2015 Elsevier Ltd. All rights reserved.

abstract This paper presents a stress-induced variable aperture model to characterise the effect of polyaxial stress conditions on the fluid flow in three-dimensional (3D) persistent fracture networks. Geomechanical modelling of the fractured rock is achieved by the finite-discrete element method (FEMDEM), which can capture deformability of matrix blocks, heterogeneity of stress fields as well as sliding and opening of pre-existing fractures. Propagation of new cracks is not required for this study of persistent fracture systems. The deformed fracture network topologies include details of dilation, opening and closing of fracture apertures, from which the local variations in hydraulic apertures are derived. Stress-controlled distribution of fracture apertures is modelled with both fracture-scale and network-scale effects considered. Under a geomechanical condition with low differential stress ratio, fracture porosity is dominated by the fracture-scale roughness. However, with the increase of stress ratio, some favourably oriented fractures are reactivated for shearing, and matrix blocks are promoted to rotate and generate large openings along their boundaries, which tend to be the key contributors to the aperture field in such persistent systems. The flow behaviour is then considered for these stressed but static solid skeletons and is investigated using a finite element solution to the Laplace problem of single-phase fluid flow. The equivalent permeability tensor of each cubeshaped rock mass is computed based on a series of flow simulations under a macroscopic pressure differential applied at opposite model boundaries with no-flow conditions on the remaining boundaries. Components of the permeability tensor are found to vary more than three orders of magnitude with respect to the change of stress ratio. Large aperture channels formed under a critical stress state accommodate significant localisation features in the flow structure of the network. The results of this study have important

Q. Lei et al. / Geomechanics for Energy and the Environment 1 (2015) 34–47

35

implications for upscaling permeability to grid block properties for reservoir flow simulation. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction Fractured rock is a naturally occurring solid material embedded with various discontinuities, such as faults, bedding planes, joints and veins. Such geological structures, along which rupture has caused mechanical weakness in the rock, often dominate the hydro-mechanical process of the host media.1 Understanding the nontrivial influence of fractures on the overall behaviour of such highly disordered geological media has important implications for many engineering applications including geothermal energy, nuclear repository safety and petroleum recovery. Discrete fracture networks are often used to mimic naturally faulted or jointed geological formations. Compared to the conventional dual porosity model2 and analytical solution for mathematically idealised discontinuity networks,3–5 the discrete fracture approach possesses the advantage of explicit representation of fracture geometries together with specific description of their intersections.6,7 Flow properties, such as block or equivalent permeability tensor, of a finite-sized fracture system can be studied from steady state fluid flow modelling.8,9 Permeability in this paper refers to the equivalent permeability that is defined as a constant tensor in Darcy’s law to represent flow in a heterogeneous medium.10 It is different from the notion of effective permeability that is considered as an intrinsic material property based on the existence of representative elementary volume (REV) at a large homogenisation scale. Effect of stress on the permeability of fractured rocks has been widely investigated using two-dimensional (2D) fracture network models. For example, Zhang and Sanderson11 analysed the stress effects on the 2D permeability tensor of three sampled natural fracture networks. The level of differential stress was found to have a significant influence on both the magnitude and direction of rock mass permeability. Min et al.12 studied the stress-dependency of rock mass permeability with the effects of non-linear joint normal deformation and shear dilation considered and they observed significant stressinduced flow enhancement along connected shear fractures. Latham et al.13 investigated the influence of in-situ stress on the permeability of an outcrop-based fracture system with consideration of bending features and crack propagation. Lei et al.14 examined the stress effect on the validity of synthetic fracture networks for representing a naturally fractured rock in terms of geomechanical and hydraulic properties. These previous studies are mostly based on 2D fracture network models, where the threedimensional (3D) nature of fluid flow in fractured rock masses under differential stress conditions remains poorly understood. The elementary object of a fracture network model is the single fracture. The commonly used approximation

for laminar flow through a single fracture is the parallel plate model, for which the hydraulic transmissivity is proportional to the cube of the distance between separated plates.3 In fact, the distribution of local apertures in a geological fracture is not uniform and strongly affected by the roughness condition of the two facing walls. Laboratory measurements of flow in single fractures revealed deviations from the results predicted by the idealised parallel plate concept due to the effect of surface asperities.15 Mismatched rough wall surfaces can result in tortuous flow through void space while bypassing asperity regions that are in contact.16 Numerical simulation of fluid flow through rough-walled joints with surface morphology characterised by fractal distributions also demonstrated the existence of channels formed by large aperture areas and barriers caused by low aperture zones.17 The laminar flow through a rough fracture may be considered following an equivalent ‘cubic’ law with the constant aperture value replaced by an appropriately weighted average, i.e. hydraulic aperture.18,19 Effects of normal compression and shearing processes on the relation between the mechanical aperture and the equivalent hydraulic aperture for rough fractures have been combined into an empirical formulation that is based on extensive results of laboratory experiments.20,21 In reviewing the literature on fracture modelling, there appears to be two distinct research focuses which depend on the chosen scale of study. The first scale is at the level of the individual fracture in which the surface roughness is described in detail, and the second scale is at the level of the fracture network with emphasis on the overall properties. Each aspect needs methods adapted to mechanisms for the given scale and appropriate for their analysis and interpretation. Up to now, there are very few attempts to bridge these two scales in 3D numerical modelling, with an exception of recent work by de Dreuzy et al.22 that combined the effects of fracture-scale heterogeneity and the network-scale topology in fluid flow modelling of 3D discontinuity systems. However, mechanical stress that has a vital impact on the variability of aperture fields were assumed uniform and isotropic across their model, regardless of the effects of fracture orientation and interaction which are known to be highly significant. In this paper, we will integrate the stress-dependent aperture model for single fractures into the hydromechanical modelling of highly interconnected fracture networks. The stress-induced variable aperture model that we propose is then applied to investigate the flow heterogeneity caused by both fracture-scale roughness and network-scale interaction effects in an idealised 3D persistent fracture network under various polyaxial stress conditions. The paper will mainly focus on the stress effect, whereas the complexity of scale effect and the possible

36

Q. Lei et al. / Geomechanics for Energy and the Environment 1 (2015) 34–47

Fig. 1. Representation of (a) a persistent fracture system by (b) an unstructured grid with (c) continuous nodal configuration inside matrix bodies, where two finite elements share identical nodes on their connecting facet, and (d) discontinuous nodal configuration on fracture interfaces, where two finite elements have different nodes for their overlapping facets.

existence of an REV are beyond the scope of this study. Persistent fracture sets are used in this paper to remove the intricacy associated with fracture propagation — a topic which is discussed briefly at the end of this paper.

by an ellipsoid visualisation technique. A brief discussion is presented on stress-induced variable apertures and localised fluid flow, after which conclusions are drawn. 2. Numerical methods

The paper is organised as follows. Section 2 describes the numerical approach for solid modelling of multiple deformable block systems, the hydraulic aperture model that computes fluid-filled equivalent fracture space from deformed solids, and the computational solution to the discrete fracture and matrix flow. In Sections 3 and 4, a 3D discontinuity network involving three orthogonal sets of persistent fractures is subjected to an isotropic lithostatic effective stress condition for initial consolidation and further to various deviatoric stress conditions. Geomechanical response of the fractured rock including stress heterogeneity, matrix block rotation, and fracture shearing and opening is investigated with the resulting stress-dependent hydraulic aperture distribution further obtained. The equivalent permeability of the stressed fracture network is calculated from single-phase steady state flow simulation, with its anisotropic property illustrated

2.1. Solid modelling The persistent discrete fracture network is integrated with the finite-discrete element method (FEMDEM)23–25 to model the geomechanical behaviour of 3D fractured rocks under polyaxial (true-triaxial) in-situ stresses. The rock mass dissected by a persistent fracture population (Fig. 1(a)) is represented by an unstructured grid system (Fig. 1(b)) involving a continuous discretisation of matrix domain using four-noded tetrahedral elements (Fig. 1(c)) and a discontinuous configuration of fracture interfaces using six-noded joint elements (Fig. 1(d)). A joint element is formed by two triangular faces that belong to opposite volumetric finite elements and are associated with separate nodes but having coincident initial coordinates (Fig. 1(d)).

Q. Lei et al. / Geomechanics for Energy and the Environment 1 (2015) 34–47

37

Fig. 2. Connectivity analysis of fracture joint elements.

Kinetics of the multi-block geological system is governed by the Cauchy linear momentum equation:

  Ω

ρ

dv dt

−∇ ·σ



dΩ =

 Ω

bdΩ +

 Γ

tdΓ

(1)

where Ω is a discrete matrix domain, Γ is the domain surface, ρ is the material density, v is the nodal velocity, σ is the Cauchy stress tensor derived using a finite strain formulation, b is the body force, and t is the surface traction force by external loads and contact interactions. The FEMDEM solid model is capable of modelling both the deformation and interaction of dissected matrix bodies under various prescribed boundary conditions. Several important geological phenomena can be simulated in the mechanical experiments, such as stress heterogeneity in matrix blocks, re-activation of shear on pre-existing fracture walls, interaction between individual fractures as well as variability of aperture distribution in the fractures.13,14 Fracture propagation is not to be modelled since only persistent fracture networks are considered in this study. 2.2. Fracture aperture model Fracture space represented by separated interfaces of deformed solids in the mechanical model is transformed to lower dimensional surfaces associated with variable equivalent apertures for fluid flow. The aperture model presented here is aimed to capture the change in fluid conduits caused by the applied in-situ stresses to the rock mass skeleton. The aperture characterisation procedure includes topological identification of the lower dimensional fracture system and calculation of variable hydraulic apertures, as described below. 2.2.1. Characterisation of fracture system topology A generic algorithm has been developed in this research for the topological diagnosis of discrete fracture systems

involving complicated interconnections, where a fracture is dissected into several block facets (polygonal shape) bounded by the intersections with many other fractures. Each block facet is further discretised into a number of connected joint elements in the FEMDEM grid system. A connectivity analysis (Fig. 2) is first implemented for each joint element to recognise its three continuously connected neighbours (i.e. sharing the same edge with identical nodes). If the edge of a joint element is located on model boundaries or fracture intersections, it is considered having no neighbour through that edge and a value of −1 is assigned numerically. Identification of isolated block facets is achieved based on a ternary-tree data structure (Fig. 3), in which a joint element is represented by a tree-node that has one parent tree-node (except the 1st level tree-node) and three child tree-nodes corresponding to its three neighbours. A breadth-first search (BFS) is conducted to recognise connected components (i.e. joint elements belonging to the same facet) by scanning the built ternary-tree structure, where previously visited tree-nodes or unreal neighbour tree-nodes are marked to be dead (i.e. empty nodes in Fig. 3) and will not grow in further searching loops. Isolated block facets represented by multiple ternarytrees are further combined based on their connectivity and coplanarity state to form corresponding discrete fractures (Fig. 4). The 3D fracture space bounded by opposite fracture walls in the solid model is transformed into a lower dimensional system represented by the median surfaces between deformed facing walls with variable apertures calculated using the approaches described in the next section. 2.2.2. Characterisation of fracture aperture distribution Stress effect on variability of fracture attributes is characterised in two different respects: (i) opening and shearing caused by network-scale fracture and matrix interaction under applied in-situ stresses (named as

38

Q. Lei et al. / Geomechanics for Energy and the Environment 1 (2015) 34–47

Fig. 3. Identification of a block facet by breadth-first search (BFS) based on a ternary-tree data structure representing the topological connectivity of joint elements.

ment (Fig. 5) given as ameso = s · n

δs =



(2)

∥s∥2 − a2meso

(3)

where s is the vector defined by the nodal difference between the barycentres of the two facing triangles of a joint element, n is the unit normal vector of the joint element. Fracture-scale joint behaviour is characterised implicitly for each joint element based on two empirical parameters: joint roughness coefficient (JRC) and joint compressive strength (JCS).20 The roughness-induced a priori initial aperture value is calculated using an empirical relation given by26 a0 =

JRC



5

0.2

σc JCS

 JRC − 0.1 = 50

(4)

where a0 is the a priori initial microscopic aperture (mm),

Fig. 4. Identified discrete fractures formed by combined block facets.

mesoscopic effect), and (ii) closure and dilation governed by fracture-scale roughness under local compressive stress and shearing movement (named as microscopic effect). Network-scale fracture opening ameso and shear displacement δs are calculated explicitly from the FEMDEM mechanical model as the normal and tangential translation components between the opposite triangles of a joint ele-

σc is the uniaxial compressive strength (MPa) and is equal

to JCS (MPa) given an assumption that the effect of weathering can be ignored. The estimated initial microscopic aperture value is assigned to all fractures equally before the loading of local normal compressive stress. Closure of microscopic apertures is modelled by a non-linear hyperbolic function (Fig. 6(a))26 amicro = a0 −

Fig. 5. A deformed joint element.

1 1

vm

+

Kni σ′ n

(5)

Q. Lei et al. / Geomechanics for Energy and the Environment 1 (2015) 34–47

39

Fig. 6. Microscopic non-linear behaviour of (a) normal closure, and (b) shear dilation of a fracture with length of 0.5 m, JRC = 15 and JCS = 120 MPa under an assumed compressive normal stress of 10 MPa.

where amicro is the current microscopic aperture (mm), Kni is the initial normal stiffness (MPa/mm) given by Kni = −7.15 + 1.75 JRC + 0.02 ×

JCS a0

(6)

and σn′ is the effective normal stress derived from the FEMDEM solid model based on the Mohr circle transformation law

σn′ = nT · σ · n

(7)

where σ is the Cauchy stress tensor of the volumetric finite element located on the fracture walls and n is the unit normal of the joint element in a column vector format, and vm in Eq. (5) is the maximum allowable closure (mm) given by

vm = −0.1032 − 0.0074 JRC + 1.1350  −0.2510 JCS × a0

(8)

δs ≤ δpeak

(9) dδ + as,peak ,

δs > δpeak

where δs is the current shear displacement, δpeak is the peak shear displacement, σn′ is the effective normal compressive stress, as,peak is the dilational displacement corresponding to the peak shear displacement. The peak shear displacement δpeak , as the threshold for dilation, is a stress- and scale-dependent parameter and can be estimated using an empirical equation27 given by

δpeak = 0.0077 × L

0.45



 × cos JRC log10

0.34

σn′ JCS



JCS

σn′

 ah =

in which the coefficients derived from experimental measurements based on numerous joint samples of five different rock types under a third loading cycle26 are adopted, because in-situ fractures are considered more likely to behave in a manner similar to the third or fourth cycle.20 The coupled effect of fracture normal and shear displacements is modelled by an improved dilation formulation27 (Fig. 6(b)) to compute shear-induced dilational displacement as as       2δs JCS 1   δ − 1 tan JRC log10 ,   3 s δpeak σn′   as =  δ     0 . 381  JCS δpeak  s tan JRC log   10 σn′ δ δpeak

where L is the size of a block facet (m)28 defined as the length of the bounded facet in the local shearing direction. The integral part in Eq. (9) is numerically approximated by the quadratic solution of the three-point Simpson’s rule. It can be noted from Fig. 6(b) that the pre-peak contraction and post-peak dilation features are well captured by this dilation model. The microscopic mechanical aperture am is calculated as the summation of roughness-induced opening amicro and shear-induced dilational displacement as . The microscopic hydraulic aperture ah defined as an equivalent aperture for laminar flow is derived based on an empirical relation with the mechanical aperture:21

δs /δpeak ≤ 0.75

,

δs /δpeak ≥ 1.0

a1/2 m JRCmob

(11)

where JRCmob is the mobilised JRC due to the roughness degradation and can be estimated using a power-base empirical relation27 as

 JRCmob = JRC

δpeak δs

0.381

.

(12)

A linear interpolation is used to determine the value of hydraulic aperture in the transition phase, i.e. 0.75 < δs /δpeak < 1.0, of Eq. (11).21 As shown in Fig. 7, in the pre-peak phase, asperities of rough walls contract first with closed small voids and increased contact areas, which leads to a slight decline in the mechanical aperture and the ratio of hydraulic aperture to mechanical aperture, i.e. ah /am . Thereafter, the fracture walls begin to dilate with asperities not destroyed yet and both mechanical and hydraulic apertures exhibit an increasing trend. In the post-peak stage, where asperities get worn and damaged, the mechanical aperture continues to increase. However, the reduction of joint porosity associated with gouge production results in a decreased ratio of ah /am , and the hydraulic aperture seems to reach a plateau under further shear displacement. The total fracture hydraulic aperture a for fluid flow is calculated as



ameso + ah0 , ah ,

 (10)

a2m /JRC2.5 ,

a=

ameso > 0 ameso ≤ 0

(13)

40

Q. Lei et al. / Geomechanics for Energy and the Environment 1 (2015) 34–47

Fig. 7. Variation of mechanical aperture, hydraulic aperture and their ratio for the fracture with length of 0.5 m, JRC = 15 and JCS = 120 MPa during a shearing process with an assumed compressive normal stress of 10 MPa.

where ah0 is the hydraulic aperture corresponding to the roughness-induced initial aperture by substituting a0 into the first part of Eq. (11). The hydraulic aperture model presented here can capture the scenario where fractures are either opened (ameso > 0) or in contact (ameso ≤ 0). By piecewisely computing the local hydraulic aperture for each joint element of the topologically identified interconnected fractures, an equivalent fracture field for fluid laminar flow is obtained. 2.3. Fluid modelling 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 are reduced to a Laplace equation as

∇ · (k ∇ p) = 0

(14)

where k is the intrinsic and isotropic permeability of the porous media with local variability permitted, and p is the fluid pressure solved at nodes of unstructured finite element grids by employing the standard Galerkin method. The element-wise constant barycentric velocity is resolved based on the pressure gradient vector field by applying Darcy’s law given by ue = −

ke

µ

∇ pe

(15)

where ue is the vector field of element-wise constant velocities, pe is the local element pressure field, µ is the constant fluid viscosity, and ke is the local permeability of a matrix volumetric element with an assumed constant value or a lower dimensional fracture element having a variable value related to the local hydraulic aperture obeying the cubic law for laminar flow between parallel plates. 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, pressure diffusion is solved for all fracture and matrix elements of the entire domain. The equivalent permeability tensor of the fractured media is

computed using element volume weighted averaging of pressure gradients and fluxes for elements e within a restricted subvolume V of the flow region away from the borders to eliminate boundary effects9

 1  V

Ve

e

uej dV e =

 kij 1  µV

e

Ve

∂ pe e dV ∂ xi

(16)

where uej is the element-wise barycentric velocity in the j direction, ∂ pe /∂ xi is the element pressure gradient along xi , and kij is the components of the symmetric second-rank permeability tensor k: kxx kyx kzx

 k=

kxy kyy kzy

kxz kyz kzz

 (17)

whose eigenvectors give the maximum, medium and minimum principal equivalent permeability, i.e. kmax , kmed , and kmin , respectively. 3. Virtual experiment setup 3.1. Persistent fracture network The discontinuity system of a periodically fractured limestone involves three orthogonal sets of persistent fractures with their geological data given by Table 1. The two vertical sets are oblique at 45° to the model boundaries where far-field horizontal stresses are to be applied. In this study, dispersion of fracture orientation is ignored to avoid treating finite elements with extremely high aspect ratios caused by intersection between subparallel fractures from the same set. All fractures are assumed through-going (i.e. only persistent fractures are modelled), tending to provide an upper limit for rock deformability and permeability. In reality, such idealised persistent networks might still be representative of some special scenarios involving highly fractured ‘non strata bound’ sedimentary rock. Assumed material properties for this fractured limestone are given in Table 2.26,29 Due to the limits of current processing power, the numerical computation is technically constrained to consider only a

Q. Lei et al. / Geomechanics for Energy and the Environment 1 (2015) 34–47

41

Fig. 8. A 0.5 m × 0.5 m × 0.5 m fracture network consisting of three orthogonal sets of persistent fractures. Table 1 Geological data of the discontinuity system consisting of three orthogonal sets of persistent fractures. Fracture sets

Dip (°)

Dip direction (°)

Spacing (m)

Set 1 Set 2 Set 3

90 90 0

45 315 0

0.050 0.075 0.100

Table 2 Material properties of the fractured limestone. Properties Rock matrix: Density Young’s modulus Poisson’s ratio Fractures: Friction coefficient Tensile strength Cohesion JCS JRC Initial aperture a0

Value

Unit

2700 30 0.27

kgm−3 GPa –

0.6 0 0 120 15 0.3

– MPa MPa MPa – mm

Table 3 Loading scheme for the geomechanical experiment.

σx′

σy′

Phase I (lithostatic stress condition): 5 5 Phase II (deviatoric stress conditions): 5 5 5 10 5 11 5 12 5 13 5 14 5 15 5 20

σz′

σx′ /σy′

5

1.0

10 10 10 10 10 10 10 10

1.0 2.0 2.2 2.4 2.6 2.8 3.0 4.0

relatively small scale virtual experiment and a 0.5 m × 0.5 m × 0.5 m cube-shaped rock sample is extracted for analysis in this paper (Fig. 8).

3.2. Procedure for numerical experiment The fractured rock is considered to be at a depth of ∼350 m with a pore fluid pressure ratio (i.e. the ratio of pore fluid pressure to lithostatic stress) equal to 0.45, producing an overburden effective stress of 5 MPa. The rock sample is designed to be surrounded by a hollow-box shaped buffer zone having a width of 0.025 m and a reduced Young’s modulus of 0.3 GPa. The buffer material has no physically corresponding substance in a realistic rock mass. It is introduced purely as a means to create boundary conditions that have a less distorting effect in the corner regions of the main volume domain of interest. The effect of the buffer zone is to provide a semi-free displacement boundary constraint to accommodate potential large slipping in such persistent system. The bottom of the model is fixed in the vertical direction, to accommodate the body force effect, but has no constraint for movements in the horizontal plane (i.e. ‘‘roller’’ boundary condition). The solid model is loaded in two consecutive phases (Fig. 9(a) and Table 3). First, an isotropic stress field (σx′ = σy′ = σz′ = 5 MPa) is imposed to consolidate the rock sample under the effective lithostatic stress. Second, a series of deviatoric stress conditions is further loaded with a fixed σx′ = 5 MPa, various σy′ = 5–20 MPa, and an increased σz′ = 10 MPa to consider the evolution of corresponding strike-slip tectonic regimes under an enhanced overburden stress (Table 3). More stress conditions are explored for the horizontal stress ratio between 2.0 and 3.0 where the state is approaching the theoretical value for frictional sliding on ideally oriented pre-existing fracture walls (i.e. a ratio of 3.1) given that the friction coefficient equals 0.6 (see page 132 in Ref. [30]). Though, in the field, observed stress ratios are generally less than 2.0, we have used values larger than this for the sake of studying the effect from typical to extreme conditions to bring out clearly the system behaviour. A larger ratio may also represent conditions close to

42

Q. Lei et al. / Geomechanics for Energy and the Environment 1 (2015) 34–47

Fig. 9. Procedure for numerical experiment: (a) mechanical modelling with polyaxial stress conditions loaded by two phases, and (b) calculation of equivalent permeability based on single-phase steady state flow tests on stressed sample under a prescribed macroscopic pressure differential imposed on each pair of opposite boundary surfaces while the remaining boundaries are impervious.

an excavation or fluid injection point. In any case, the simulations may correspond to laboratory measurements where such stress ratios may be imposed. Single-phase steady state fluid flow through the deformed fracture network with stress-induced variable apertures is further modelled by imposing classical permeameter boundary conditions: two opposite boundary surfaces of the cube have fixed heads while the four orthogonal boundaries parallel to the flow direction are impervious (Fig. 9(b)). 4. Results 4.1. Fracture apertures The model that has arrived at equilibrium under the initial isotropic stress condition further adjusts to a new deformed state when various deviatoric stress fields are loaded. The stress ratio of σy′ to σx′ triggers stress heterogeneity in the matrix blocks (Fig. 10(a)), shear displacements along the two vertical sets (Fig. 10(b)), and even mesoscopic fracture openings caused by block rotations if the stress ratio is high enough (Fig. 10(c)). Hydraulic aperture of the stressed fracture networks is calculated as the summation of mesoscopic opening caused by fracture interaction and block rotation, and microscopic aperture governed by the surface roughness nature (Eq. (13)). Effect of stress generates significant fracture-scale heterogeneity for the distribution of hydraulic apertures in single fractures. Fig. 11 shows the heterogeneous aperture contour of a vertical fracture extracted from the network under the polyaxial stress condition with σy′ /σx′ = 3. Very large apertures are clustered in some local areas, which seem to be connected and form a slightly diverted vertical channel from the top to the bottom of the model.

Fig. 12 shows the network-scale distribution of hydraulic apertures in log scale under different polyaxial stress conditions. In the stress case of σy′ /σx′ = 1, hydraulic apertures are uniformly distributed and exhibit quite low magnitude in such an isotropic stress field. With the increase of the far-field stress ratio, heterogeneity of fracture apertures begins to emerge and develop. Especially in cases of σy′ /σx′ ≥ 3, very large hydraulic apertures are localised in some fractures of the two vertical sets that are favourably oriented for shearing. Fracture porosity is calculated as the proportion of fracture hydraulic aperture space to the total rock mass volume. The contributions from mesoscopic and microscopic effects are distinguished to isolate the sources of hydraulic apertures under different stress conditions (Fig. 13). In the case with a low stress ratio, e.g. σy′ /σx′ < 2.5, fracture porosity is mainly dominated by the microscopic roughness effect. As the stress ratio increases, the microscopic component exhibits moderate increase due to shear dilatancy, while the mesoscopic counterpart begins to manifest itself by a dramatic growth. As a result, the total porosity shows a continuous increasing trend under the increased differential stress ratio. It seems that the microscopic and mesoscopic porosity components as well as the total porosity display a positive linear relation with the stress ratio when σy′ /σx′ > 2.5. 4.2. Equivalent permeability Matrix permeability km is assumed to have a low value, i.e. 1 × 10−15 m2 , to produce a high fracturematrix permeability contrast and impose a condition close to fracture-only flow. Poroelastic effect of the Biot-type coupling of pore fluid pressure and solid elastic stress31 is only modelled for a particular scenario with the Biot

Q. Lei et al. / Geomechanics for Energy and the Environment 1 (2015) 34–47

43

Fig. 10. (a) Distribution of differential stress in the matrix blocks, (b) distribution of fracture shear displacement in log scale, and (c) fracture openings caused by block rotations (observed from the top of the model) of the fractured rock under the polyaxial stress condition of σx′ = 5 MPa, σy′ = 15 MPa, and σz′ = 10 MPa.

Fig. 11. Distribution of hydraulic aperture within a single fracture under the polyaxial stress condition of σx′ = 5 MPa, σy′ = 15 MPa, σz′ = 10 MPa.

Fig. 12. Distribution of hydraulic apertures in the fracture network under various polyaxial effective stress conditions: (a) σy′ /σx′ = 1, (b) σy′ /σx′ = 2, (c) σy′ / σx′ = 3, (d) σy′ /σx′ = 4, given that σx′ = 5 MPa and σz′ = 10 MPa.

44

Q. Lei et al. / Geomechanics for Energy and the Environment 1 (2015) 34–47

Fig. 13. Fracture porosity of the fractured rock under various polyaxial stress conditions: three curves represent the porosity induced by mesoscopic effect, microscopic effect and the value of total porosity, respectively.

Fig. 14. Equivalent permeability of the fractured rock under various polyaxial in-situ stress conditions.

coefficient for the solid skeleton compressibility equal to 1.0. The equivalent permeability of the fractured rock under various polyaxial stress conditions is derived from the steady state flow simulation, where a subvolume is conservatively chosen with a distance of 10% of the model size away from the nearest domain boundaries. As shown in Fig. 14, the increased stress ratio of σy′ to σx′ leads to considerable increase over several orders of magnitude in the diagonal of the permeability tensor, i.e. components, kxx , kyy , and kzz . A transition regime with steep permeability increase occurs when the far-field stress ratio is approaching the critical threshold, i.e. 3.1. The permeability tensor is visualised as a triaxial ellipsoid with three semi-principal axes indicating the magnitudes of maximum, medium, and minimum principal permeability, i.e. kmax , kmed , and kmin , respectively (Fig. 15). Normalisation is performed with respect to corresponding kmax since the absolute values span several orders of magnitude. In the case of σy′ /σx′ = 1, the permeability tensor ellipsoid is quite isotropic, despite of the intrinsic anisotropy

in fracture geometries. In the case with higher stress ratios, e.g. the one of σy′ /σx′ = 3, significant permeability anisotropy is induced by the deviatoric stress acting with respect to the favourably oriented vertical fractures, resulting in a very high permeability in the subvertical direction. The increased far-field stress ratio also leads to considerable change in flow patterns as illustrated by Fig. 16. In the case of σy′ /σx′ = 1, fluid spreads through the whole network due to the quite uniformly distributed apertures. However, in the case of σy′ /σx′ = 3, fluid flow is localised in some zigzag-shaped pathways corresponding to the large aperture channels formed by parts of some fractures of the two vertical sets. 5. Discussion Stress-induced heterogeneity of hydraulic apertures of a 3D persistent fracture network has been modelled with consideration of both fracture-scale and network-scale effects. In cases with lower stress ratio, fracture porosity

Q. Lei et al. / Geomechanics for Energy and the Environment 1 (2015) 34–47

(a) σy′ /σx′ = 1.

45

(b) σy′ /σx′ = 3.

Fig. 15. Ellipsoid visualisation, after normalisation with respect to kmax , of the permeability tensor of the fractured rock under different polyaxial stress conditions: (a) σx′ = 5 MPa, σy′ = 5 MPa, σz′ = 10 MPa, (b) σx′ = 5 MPa, σy′ = 15 MPa, σz′ = 10 MPa. Note kmax in (b) is >1000 times kmax in (a).

Fig. 16. Flow pathways of the fractured rock models under different polyaxial stress conditions: (a) σx′ = 5 MPa, σy′ = 5 MPa, σz′ = 10 MPa, (b) σx′ = 5 MPa, σy′ = 15 MPa, σz′ = 10 MPa. (Note the flow arrow sizes representing local flux magnitudes in the flow test of the case σy′ /σx′ = 3 are scaled down by a factor that is 20, 50, and 100 times the one of the case σy′ /σx′ = 1 for east-to-west, north-to-south, and top-to-bottom pattern, respectively.)

is mainly controlled by the fracture-scale microscopic roughness effect. With the increase of stress ratio, preexisting fractures were reactivated for shearing and matrix blocks were mobilised into rotation and sliding at the mesoscale, which created some large openings along block boundaries. As a result, even in the persistent fracture, local hydraulic apertures can vary greatly, as shown in Fig. 11. The formation of large aperture channels due to such network-scale mechanical interactions leads to

significant flow localisation and dramatic increase of overall hydraulic conductivity. The transition stage of permeability with steep growth that occurred when the far-field stress ratio is approaching the critical threshold (Fig. 14) shows consistency with the results of 2D fracture network modelling.12 The results of the case under a critically stressed state, e.g. σy′ /σx′ = 3, are of particular interest. First, the shear displacement is extremely heterogeneous, in spite of

46

Q. Lei et al. / Geomechanics for Energy and the Environment 1 (2015) 34–47

being given such regular geometrical configurations of fracture sets. The system finds equilibrium by activating sliding with local extremes of shear displacement as highlighted in Fig. 10(b). Locally, the sliding on the two vertical sets has created large aperture channels parallel to the active fractures (Fig. 11), which shows consistency with the field observation from boreholes that critically stressed faults with favourable orientations appear to have much higher hydraulic conductivity.30 The result supports what is already known of the strike-slip faulting regime, that significantly higher permeability can be anticipated in the vertical direction associated with localised flow along displacing and dilating fractures.32,33 This raises the question of whether the imposed boundary conditions with orthogonally applied stresses and semi-free displacement constraints are the most appropriate for modelling mechanical behaviour of the rock sample with such persistent fractures and whether the localisation effect is exaggerated by considering a domain with such few idealised fractures. However, some fundamental mechanisms captured in this idealised fractured rock model, e.g. stress-induced fracture dilation, block rotation and flow localisation, would probably exist in more complicated systems having arbitrarily shaped and oriented fractures, which has been proven in many 2D models.11–14 The high sensitivity of equivalent permeability (Figs. 14 and 15) and flow structure (Fig. 16) to the polyaxial stress condition indicates that special attention is required when the in-situ stress state of rock masses is significantly perturbed resulting from natural or human activities. For example, unloading effects during the excavation of underground infrastructures may cause significant stress redistribution surrounding the openings; injections and extractions of fluids during oil/gas reservoir production can significantly change the pore fluid pressure level and further vary the effective stress state of rock masses; multiple complex factors (e.g. underground excavation, radioactivity-induced heat transfer, and glaciation loading) can engender remarkable changes on the geomechanical condition of nuclear waste repositories. Such irreversible perturbations may lead to intensive fault reactivation, dramatic flow enhancement, and severe construction risk. Needless to say, the modelling methods employed in this study can also be applied to investigate more permeable matrix rock scenarios. The extreme nature of the flow anisotropy would be somewhat ameliorated by modelling a rock system with a more permeable matrix material. However, the realism of the ‘persistent-only’ fracture model with no new fracturing adopted here may come into question as such a rock mass with weak sedimentary rock properties may be weak enough to locally propagate new fractures before the exceptionally high in-situ stress ratios considered here could be generated. Unlike some other conventional 3D upscaling methods in the literature3–5,10 that do not require explicit mechanical and flow simulations to compute the equivalent hydraulic properties, the proposed approach here may not be a practical solution under the limits of current processing power. This is a particular problem for applications to real reservoirs with domains spanning hundreds or thousands of metres and consisting of millions of fractures. However,

this study still has important implications for upscaling permeability to grid block properties for 3D reservoir flow simulation. For example, the results obtained in this study imply that determination of an REV size, if it exists, may be a sophisticated process that requires many coupled effects to be considered in the model including not only the description of geometrical features, but also, characterisation of the geomechanical setting and changes resulting from any perturbation of the stress field. Indeed, it is recognised that there is unlikely to be an REV once a realistic system with impersistent fractures is modelled due to its intrinsic fractal nature.34 Future work will involve modelling of these more general networks which have fractures with arbitrary shapes, distributed sizes, and dispersed orientations. A 3D crack propagation model35 will be employed to capture the brittle deformation response including local concentrations of critically high tensile or differential stresses, together with realistic fracture opening and shearing behaviour on both pre-existing and newly propagated fractures. Such capability opens the way to modelling 3D flows in geomechanically realistic fractured rocks, including rock mass regions with locally higher density fractures or fracture corridors. Upscaling fracture network hydro-mechanical properties based on more physically realistic modelling of two-way coupling process31 is also a nontrivial issue to be resolved. 6. Conclusion To conclude, this paper presented a stress-induced variable aperture model to capture the effect of polyaxial stress conditions on the flow properties of 3D persistent fracture networks. Geomechanical behaviour of the rock mass was simulated by the FEMDEM solid model, where a fracture treated as the interface between discrete matrix bodies can open, shear and dilate in the heterogeneous stress field. Under the stress condition with relatively lower differential stress ratio, fracture apertures are mainly governed by the fracture-scale roughness effect. With the increase of in-situ stress ratio, fractures with favourable orientations are reactivated to shear and matrix blocks bounded by the shearing fractures are promoted to rotate, which generates significant fracture openings at the block boundaries. Such fracture openings tend to be the dominant contributor to the aperture field in the high stress ratio conditions. To prepare for the flow modelling required in this work, a new discrete fracture system indexing logic was developed based on a breadth-first search of ternary-tree structures to systematically identify the fracture network topology and its associated apertures. Based on a series of single-phase flow virtual experiments, equivalent permeability of the stressed fractured rock was computed, ranging over more than three orders of magnitude with respect to the variation of effective stress ratio. A near-isotropic permeability tensor was observed in the case with lower stress ratio, whereas the fractured rock under a critical stress state exhibits highly anisotropic features in its permeability. Fluid flow tends to localise in some critically stressed fractures that are associated with much higher hydraulic conductivity than other fractures that are not significantly reactivated for shearing. The large aperture

Q. Lei et al. / Geomechanics for Energy and the Environment 1 (2015) 34–47

channels that are optimally oriented with regard to the direction of pressure gradient provide a major pathway for fluid migration. The results of this study have important implications for upscaling permeability to grid block properties for reservoir flow simulation as well as other relevant engineering problems. Acknowledgements The authors would like to thank the sponsors of the itf-ISF project ‘‘Improved Simulation of Faulted and Fractured Reservoirs’’ and to acknowledge the Janet Watson scholarship, awarded to the first author by the Department of Earth Science and Engineering, Imperial College London. The authors would also thank Golder Associates Inc. for providing an academic licence of FracMan7 for this research. References [1] Zimmerman RW, Main I. Hydromechanical behavior of fractured rocks. Mechanics of Fluid-Saturated Rocks. In: Gueguen Y, Bouteca M, eds. London: Elsevier; 2004:363–421. [2] Warren JE, Root PJ. The behavior of naturally fractured reservoirs. SPE J. 1963;3:244–255 http://dx.doi.org/10.2118/426-PA. [3] Snow DT. Anisotropic permeability of fractured media. Water Resour. Res. 1969;5:1273–1289. http://dx.doi.org/10.1029/WR005i006p01273. [4] Oda M. Permeability tensor for discontinuous rock masses. Géotechnique 1985;35:483–495. http://dx.doi.org/10.1680/geot.1985.35.4.483. [5] Oda M. An equivalent continuum model for coupled stress and fluid flow analysis in jointed rock masses. Water Resour. Res. 1986;22: 1845–1856 http://dx.doi.org/10.1029/WR022i013p01845. [6] Dershowitz WS, Einstein HH. Characterizing rock joint geometry with joint system models. Rock Mech. Rock Eng. 1988;21:21–51 http://dx.doi.org/10.1007/BF01019674. [7] Dershowitz B, LaPointe P, Eiben T, Wei L. Integration of discrete feature network methods with conventional simulator approaches. SPE J. 2000;3:165–170 http://dx.doi.org/10.2118/62498-PA. [8] Pouya A, Fouché O. Permeability of 3D discontinuity networks: New tensors from boundary-conditioned homogenisation. Adv. Water Resour. 2009;32:303–314. http://dx.doi.org/10.1016/j.advwatres.2008.08.004. [9] Lang PS, Paluszny A, Zimmerman RW. Permeability tensor of threedimensional fractured porous rock and a comparison to trace map predictions. J. Geophys. Res. 2014;199:6288–6307. http://dx.doi.org/10.1002/2014JB011027. [10] Renard P, de Marsily G. Calculating equivalent permeability: a review. Adv. Water Resour. 1997;20:253–578. http://dx.doi.org/10.1016/S0309-1708(96)00050-4. [11] Zhang X, Sanderson DJ. Effects of stress on the two-dimensional permeability tensor of natural fracture networks. Geophys. J. Int. 1996;125:912–924. http://dx.doi.org/10.1111/j.1365-246X.1996.tb06034.x. [12] Min KB, Rutqvist J, Tsang C-F, Jing L. Stress dependent permeability of fractured rock masses: a numerical study. Int. J. Rock Mech. Min. Sci. 2004;41:1191–1210. http://dx.doi.org/10.1016/j.ijrmms.2004.05.005. [13] Latham J-P, Xiang J, Belayneh MW, Nick HM, Tsang C-F, Blunt MJ. Modelling stress-dependent permeability in fractured rock including effects of propagating and bending fractures. Int. J. Rock Mech. Min. Sci. 2013;57:100–112. http://dx.doi.org/10.1016/j.ijrmms.2012.08.002.

47

[14] Lei Q, Latham J-P, Xiang J, Tsang C-F, Lang P, Guo L. Effects of geomechanical changes on the validity of a discrete fracture network representation of a realistic two-dimensional fractured rock. Int. J. Rock Mech. Min. Sci. 2014;70:507–523. http://dx.doi.org/10.1016/j.ijrmms.2014.06.001. [15] Witherspoon PA, Wang JSY, Iwai K, Gale JE. Validity of cubic law for fluid flow in a deformable rock fracture. Water Resour. Res. 1980;16: 1016–1024 http://dx.doi.org/10.1029/WR016i006p01016. [16] Zimmerman RW, Chen D-W, Cook NGW. The effect of contact area on the permeability of fractures. J. Hydrol. 1992;139:79–96 http://dx.doi.org/10.1016/0022-1694(92)90196-3. [17] Brown SR. Fluid flow through rock joints: The effect of surface roughness. J. Geophys. Res. 1987;92:1337–1347. http://dx.doi.org/10.1029/JB092iB02p01337. [18] Tsang YW, Witherspoon PA. Hydromechanical behavior of a deformable rock fracture subject to normal stress. J. Geophys. Res. 1981;86:9287–9298 http://dx.doi.org/10.1029/JB086iB10p09287. [19] Renshaw CE. On the relationship between mechanical and hydraulic apertures in rough-walled fractures. J. Geophys. Res. 1995; 100:24629–24636 http://dx.doi.org/10.1029/95JB02159. [20] Barton N, Bandis SC, Bakhtar K. Strength, deformation and conductivity coupling of rock joints. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 1985;22:121–140 http://dx.doi.org/10.1016/01489062(85)93227-9. [21] Olsson R, Barton N. An improved model for hydromechanical coupling during shearing of rock joints. Int. J. Rock Mech. Min. Sci. 2001; 38:317–329 http://dx.doi.org/10.1016/S1365-1609(00)00079-4. [22] de Dreuzy J-R, Méheust Y, Pichot G. Influence of fracture scale heterogeneity on the flow properties of three-dimensional discrete fracture networks (DFN). J. Geophys. Res. 2012;117:B11207 http://dx.doi.org/10.1029/2012JB009461. [23] Munjiza A. The Combined Finite-Discrete Element Method. London: Wiley; 2004. [24] Munjiza A, Knight EE, Rougier E. Computational Mechanics of Discontinua. Chichester: John Wiley & Sons; 2011. [25] Xiang J, Munjiza A, Latham J-P. Finite strain, finite rotation quadratic tetrahedral element for the combined finite-discrete element method. Internat. J. Numer. Methods Engrg. 2009;79:946–978 http://dx.doi.org/10.1002/nme.2599. [26] Bandis SC, Lumsden AC, Barton NR. Fundamentals of rock joint deformation. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 1983;20: 249–268 http://dx.doi.org/10.1016/0148-9062(83)90595-8. [27] Asadollahi P, Tonon F. Constitutive model for rock fractures: Revisiting Barton’s empirical model. Eng. Geol. 2010;113:11–32 http://dx.doi.org/10.1016/j.enggeo.2010.01.007. [28] Bandis SC, Lumsden AC, Barton NR. Experimental studies of scale effects on the shear behaviour of rock joints. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 1981;18:1–21 http://dx.doi.org/10.1016/01489062(81)90262-X. [29] Lama RD, Vutukuri VS. Handbook on Mechanical Properties of Rocks: Testing Techniques and Results, Vol. 2. Clausthal: Trans. Tech. Publications; 1978. [30] Zoback M. Reservoir Geomechanics. New York: Cambridge University Press; 2007. [31] Rutqvist J, Stephansson O. The role of hydromechanical coupling in fractured rock engineering. Hydrogeol. J. 2003;11:7–40 http://dx.doi.org/10.1007/s10040-002-0241-5. [32] Sibson RH. Crustal stress, faulting and fluid flow. Geol. Soc. London Spec. Publ. 1994;78:69–84. http://dx.doi.org/10.1144/GSL.SP.1994.078.01.07. [33] Sanderson DJ, Zhang X. Critical stress localization of flow associated with deformation of well-fractured rock masses, with implications for mineral deposits. Geol. Soc. London Spec. Publ. 1999;155:69–81 http://dx.doi.org/10.1144/gsl.sp.1999.155.01.07. [34] Bonnet E, Bour O, Odling NE, Davy P, Main I, Cowie P, Berkowitz B. Scaling of fracture systems in geological media. Rev. Geophys. 2001; 39:347–383 http://dx.doi.org/10.1029/1999RG000074. [35] Guo L, Latham J-P, Xiang J. Numerical simulation of breakages of concrete armour units using a three-dimensional fracture model in the context of the combined finite-discrete element method. Comput & Structures 2014;146:117–142. http://dx.doi.org/10.1016/j.compstruc.2014.09.001.