Hydro-mechanical interaction effects and channelling in threedimensional fracture networks undergoing growth and nucleation

Adriana Paluszny, Robin N. Thomas, Maria C. Saceanu, Robert W. Zimmerman

Department of Earth Science and Engineering, Imperial College London, London, SW7 2AZ, UK

ABSTRACT The flow properties of geomechanically generated discrete fracture networks are examined in the context of channelling. Fracture networks are generated by growing fractures in tension, modelling the low permeability rock as a linear elastic material.Fractures are modelled as discrete surfaces which grow quasi-statically within a three-dimensional(3D)volume.Fractures may have their locations specified as a simulation input, or be generated as a function of damage, quantified using the local variation in equivalent strain. The properties of the grown networks are shown to be a product of in situ stress,relative orientation of initial flaws,and competitive process of fracture interaction and growth.Fractures grow preferentially in the direction perpendicular to the direction of maximum tension and may deviate from this path due to mechanical fracture interaction.Flow is significantly channelled through a subset of the fractures in the full domain,consistent with observations of other real and simulated fractures.As the fracture networks grow, small changes in the geometry of the fractures lead to large changes in the locations and scale of primary flow channels. The flow variability and formation of channels are examined for two growing networks,one with a fixed amount of fractures,and another with nucleating fractures. The interaction between fractures is shown to modify the local stress field, and in turn the aperture of the fractures. Pathways for single-phase flow are the results of hydro-mechanical effects in fracture networks during growth.These are the results of changes to the topology of the network as well as the result of mechanical self-organisation which occurs during interaction leading to growth and intersection.

Keywords:Fracture growth Fracture mechanics Stress intensity factor (SIF)Quasi-static growth Finite elements

1. Introduction

The study of fractured media in geological formations is challenging, mainly due to the large span of scales and the variety of physical processes involved. Fractures in porous media can act as connected flow paths for fluids, increasing the permeability of the medium and controlling the migration of fluids through the fractured subsurface.Modelling fluid flow in fractured porous rocks has practical applications in hydrogeology (Black et al., 2017), nuclear waste disposal (da Silveira et al., 2013; Kim et al., 2019), CO2sequestration (Bigi et al., 2013; Salimzadeh et al., 2018b),geothermal energy production(Fox et al.,2015;Wang et al.,2019),and fluid storage and retrieval from the subsurface (Bauer et al.,2017). The study of fluid flow in fractured porous media is, however, not limited to the subsurface but has a wider spectrum of applications, including applications in brittle and quasi-brittle engineered materials for civil engineering and industrial design(Carmeliet et al., 2004; Kim et al.,2016; Thomas et al.,2018).

In fractured rocks, channelling can be defined as the formation of a preferential pathway through a heterogeneous fractured rock(Tsang and Neretnieks,1998). A channel is regarded as the path of least resistance, influencing transport in fractured media by accelerating the migration of fluids in the subsurface, a phenomenon that has been studied extensively by many authors both experimentally and theoretically (e.g. Tsang and Neretnieks, 1998).Experimental studies of the permeability of single 0.05—0.3 m granitic fractures have shown the formation of channels as a function of contact area and aperture redistribution (Watanabe et al., 2015). A three-dimensional (3D) numerical single fracture study in which apertures are assigned in a stochastic manner also showed a significant channelling conducing to higher permeability for thermally stimulated fractures (Guo et al., 2016). Discrete fracture network (DFN) studies, primarily based on stochastically generated fracture networks, have investigated analytically and numerically the flow and transport properties of these networks,finding strong links between permeability, density and connectivity of these networks (Adler et al., 2012; Ishibashi et al., 2019).Fracture connectivity is generally regarded as one of the fundamental mechanisms of channelling, as the interconnected pathways are expected to form a series of competing routes for fluids to migrate through(Le Goc et al.,2010).Therefore,flow channelling in fracture networks in low permeability rocks is expected to depend not only on fracture geometry and network topology, but also on small-scale effects, such as surface roughness and material fill resulting from chemical reactions (Lang et al., 2015; Zambrano et al., 2019), and large-scale effects, such as in situ stresses that modify how apertures are distributed (Bisdom et al., 2016; Lang et al., 2018). In inter-bedded fractured media, the contrast between mechanically strong and weak layers may lead to aperture variations within fractures that induce channelling of flow along the weakest paths (Brenner and Gudmundsson, 2004). Channels frequently comprise a much smaller proportion of the full network,complicating the analysis of permeability and percolation in the subsurface, as not all fractures will provide a pathway for flow.

Channels are shown to evolve during growth as a function of both fracture interaction and aperture redistribution. The present work investigates channelling as a function of fracture nucleation and growth, in geomechanically generated DFNs. Fractures nucleate, interact, and grow in response to boundary stresses,and apertures and flow are resolved using a fully monolithically coupled poroelastic approach.The resulting flow patterns are influenced by fracture interaction, connectivity, and aperture redistribution effects simultaneously, without incorporating any further heterogeneity into the simulation.This yields strong channelling patterns at the fracture and network scale.

2. Modelling fracture growth

Several methods have been developed to capture discrete fracture and fault discontinuities including extinction methods, debonding methods, and fracture-growth-based methods (Jing,2003). These can be implemented using a range of numerical methods, such as the finite element method (FEM), the boundary element method (BEM), the discrete/distinct element method(DEM), and extensions and specialisations thereof, such as the extended FEM(XFEM).In general,fracture growth methods tend to be implemented using FEM,BEM,and XFEM;de-bonding methods are associated with DEM implementations;and extinction methods can be implemented in both FEM and DEM approaches. FEM- and BEM-based methods solve the fundamental equations of linear elastic mechanics, by directly taking the elastic and inelastic material properties of the bodies as inputs(Cook et al.,2001;Kojic and Bathe, 2006). DEM requires extensive calibration of material and flow properties (Lisjak and Grasselli, 2014), with multiple triaxial tests performed to calibrate the properties of the sub-elements of the DEM mesh for a specific mesh size and sample size. FEM-DEM approaches seek to enrich the classic DEM formulation by modelling each sub-element as a finite element domain (e.g. Munjiza,2004), thus capturing inelastic, thermal, and flow behaviour within each sub-element.However,FEM-DEM does not remove the usual dependence of the DEM on the de-bonding growth approach,which is restricted to a pre-existing mesh.

In many applications, the growth path is not a variable of interest because it is known a priori, such as in delamination of an airplane wing component, or failure of a cement column. These applications utilise numerical methods in which fracture growth is guided by a pre-existing mesh or pre-defined path,such as the case for cohesive crack methods. However, when investigating the growth of multiple fractures in geological materials, crack paths directly influence the ensuing stress field, and strongly affect the resulting fracture pattern, in a manner that may or may not be predictable a priori. DEM and node and element extinction methods grow fractures along a selection of pre-existing element boundaries, often resulting in mesh refinement and mesh form bias. Fracture surfaces can be modelled independently from the mesh, as in BEM (Renshaw and Pollard,1994), FEM (Paluszny and Zimmerman, 2011; Davis et al., 2014), XFEM (Mo?s et al., 1999;Rabczuk et al.,2010),mesh-less methods(Bordas et al.,2008),and peri-dynamics(Ha and Bobaru,2010).FEM requires re-meshing to cope with the change in geometry,whilst XFEM uses a fixed ‘level set’ function with a restricted amount of possible fracture orientations to represent a fracture segment within each element.XFEM relies on a sub-element representation of discontinuity using enriched shape functions and is ideal for modelling the moving boundaries of two bodies, such as in fluid—solid interaction cases,and to represent material interfaces,such as when modelling solid inclusions, woven materials, or subsurface channels. When modelling fractures, XFEM elements containing tips have specific tip-enriched functions, and in two dimensions, elements that handle intersections have specific ‘blended’ functions.

Continuum mechanics methods allow defining a rock matrix that has a volumetric variation of properties and contains fractures as embedded discontinuities defined within the continuum. In all of these approaches, the volume is assumed to be discretised by a surface mesh in two dimensions, or a volume mesh in three dimensions.Two distinct numerical approaches can be used to model fractures within a continuum: (1) non-geometric methods, in which the fracture is represented as a material property of the mesh (including localised damage, smear, phase-field, and plastic models solved with FEM); and (2) geometric methods, in which each fracture is represented explicitly, and has a volumetric or surface domain that is distinct from the rock matrix (as in XFEM,DEM, and FEM with mesh discontinuities).

In geometric methods, specific properties can be directly assigned to the fracture, such as friction and cohesion. Processes such as fluid flow and frictional sliding can be directly examined.In geometric methods, discrete fractures are represented as lines,surfaces,or volumes,and have a distinct domain,which is assumed to be independent of the rock matrix definition. Geometric methods differ in the manner in which apertures are handled. In particular, when modelling coupled processes, apertures have a strong effect on flow, and should be tracked as a distribution of properties along fracture surfaces as opposed to single or averaged properties derived from length or remote stresses.

By modelling fractures using linear elastic fracture mechanics(LEFM), geological patterns have been reproduced, including fractures and faults in layered systems in both two dimensions(Renshaw and Pollard, 1994; Sch?pfer et al., 2007) and three dimensions(Paluszny and Zimmerman,2013),and single(Bremberg and Dhondt,2009)and multiple(Paluszny and Zimmerman,2011)fracture propagation and interaction in three dimensions. The present work models fractures geometrically. Fractures are represented discretely using a combination of smooth and tessellated curves and surfaces,and represented as a discontinuity in the flow field, which allows modelling the flow within fractures independently of, and in combination with, flow within the matrix. To handle growth along each fracture’s tip, three stress intensity factors (SIFs) are calculated, quantifying the energy that deforms the tip. The discretisation of the domain conforms to the evolving fracture geometry in each step and is automatically updated to follow the changing internal geometry of the domain,constituting an adaptively re-meshed domain. Therefore, fracture extension is not restricted to the existing discretisation.Instead,fractures grow along the paths determined by the local stress state at their tip,and growth is influenced by the interaction between fractures. This method has been applied in three dimensions to the analysis of fracture interaction(Thomas et al.,2017),single(Salimzadeh et al.,2017a)and multiple(Salimzadeh et al.,2017b)hydraulic fracturing within a poroelastic medium, and the analysis of thermo-hydromechanical effects of geothermal production on fracture aperture distribution (Salimzadeh et al., 2018a). In the simulator, fracture and fault geometries are stored independently from the mesh(Paluszny and Zimmerman, 2013), and thus the fracture growth is not restricted to conform to an existing mesh,nor is it constrained by any mesh. The FEM solves the deformation and flow in the system, taking into account the variable properties of the volumetric and surface bodies representing the rocks and surfaces in the model.Fig.1 shows an example of 100 growing fractures which interact to form a multi-scale pattern(further explained in Section 3.3).

2.1. Governing equations

Patterns are generated by a finite element-based discrete fracture propagation simulator, the Imperial College Geomechanics Toolkit(ICGT).The library implements a geometric fracture growth method, in which deformation is numerically calculated based on measurable material properties,such as Young’s modulus,Poisson’s ratio, and friction coefficient. The library also performs fully coupled thermo-poroelastic simulations to model hydraulic fracturing in a poroelastic matrix,with full 3D treatment of fluid leakoff, in layered systems, and in the context of geomechanical implications of sub-surface fluid injection and extraction. The displacement field is calculated by solving 3D linear elasticity governing equations with the FEM, considering surface tractions and friction between fracture walls:

Fig.1. Fractures in a 40 m×40 m×20 m granitic rock,for(a,b)for growth step 5(c,d)growth step 10(e,f)growth step 15,and(g,h)growth step 20(a—d)Top view and(e—h)side view of a fracture network with 100 fractures, which nucleated and grew as a function of remote anisotropic stresses.

whereDis the drained stiffness matrix; ε is the strain tensor;α is the Biot coefficient;pmis the matrix fluid pressure;Iis the second order identity matrix;Fis the body force per unit volume;pfis the fluid pressure in the fracture;andncis the unit normal vector to the fracture surface.

The volumetric mesh is split at the fractures, creating a conforming surface mesh that discretises the fracture domain. The contact between these fracture walls is resolved using a gap-based augmented Lagrangian approach(Nejati et al.,2016),that captures both stick and slip regimes. To facilitate computation of relative displacement and contact between the two surfaces of a fracture,each fracture surface element is paired with the equivalent surface on the other side of the fracture, with one side referred to as the‘master’ and the other as the ‘slave’. Discretisation of fracture surfaces always generates two coincident meshes, thus master—slave pairs always have the same undeformed shape, and only share nodes if they are along a fracture’s tip. Fracture element penalties are defined as a function of local relative fracture node displacement, with a minimum at the tips. Once the displacement field is computed, stresses and strains are calculated assuming that the matrix is homogeneous, linear elastic and isotropic.

The flow properties of the fractured rock are calculated numerically for geomechanical 3D discrete fracture patterns.FEM is used to compute flow in the system, assuming Darcy flow in the matrix and laminar flow in the fractures,as described by Lang et al.(2014). Fluid flow through the matrix is given by

wherekmis the permeability tensor of the rock matrix; μwis the fluid viscosity;ρwis the fluid density;gis the gravitational acceleration vector;uis the displacement vector;cwis the fluid compressibility; α is the Biot coefficient; φ is the porosity of the rock;Ksis the bulk modulus of the matrix;is the flow rate through the matrix;knis the matrix permeability normal to the fracture surface; andis the pressure gradient perpendicular to the fracture.Fracture apertures are determined with the separation between the two sides of each fracture, and are influenced by the distribution of stress, fracture orientation, and mutual fracture interaction:

whereafis the fracture aperture,u+andu-are the displacements of the two faces of the fracture. Apertures are tracked for each fracture node and stored in the form of a scalar variable. Local fracture permeability is calculated using the cubic law, for each fracture triangle element. Darcy flow simulations are used to examine the distribution of flow and domain equivalent permeability.The governing equation for fluid flow through the fracture is given by

Fracture apertures are responsive to the history of deformation of the fracture surface and result in an uneven distribution of flow along the fracture surface in the form of channelling (Lang et al.,2016). This leads to a lower permeability of the fracture, as compared to the length-dependent permeability assumption(Lang et al.,2015).There are multiple mechanisms,in addition to fracture network topology and connectivity, which may strongly affect the enhancement or reduction of permeability during growth and intersection (Paluszny and Matthai, 2010). In order to investigate these effects, detailed multi-fracture growth simulations are performed yielding mechanically consistent aperture distributions,and subject them to representative in situ stresses(Paluszny et al.,2015).

2.2. Flow and permeability

The equivalent permeability tensor,K, is calculated using the method of Lang et al. (2014), which relies on element-wise averaging of pressure and fluid velocity. The pressure is obtained from the finite element solution to the Laplace equation, assuming incompressible, single-phase flow:

where Δpis the pressure gradient andkeis the element permeability. Using Darcy’s law, pressure and fluid velocity at each element are related to the permeability:

where νeis the velocity vector at the centre of an elemente,μ is the fluid viscosity, Δpeis the element pressure gradient. Eq. (6) can then be extended to three dimensions to determine the full equivalent permeability tensor,K,as explained by Lang et al.(2014).Three flow experiments must be conducted to obtain the permeability tensor in three dimensions. The boundary conditions for these flow experiments are a fluid pressure gradient of 200 Pa—100 Pa between opposite boundaries, and no flow on the other boundaries. In the case of a box domain, velocity and pressure gradient averages are calculated in the directions of the three global coordinate axes. In the case of a radial domain, the effective radial permeability is determined by computing flux and pressure gradient averages for the radial direction (Thomas et al.,2018).

SIFs for the three principal deformation modes (KI,KIIandKIII)are calculated using the interaction integral and the virtual integration disk method(Nejati et al.,2015a).This method solves theI-integral numerically over a two-dimensional (2D) virtual disk domain, embedded in the 3D matrix domain, providing the strain energy release rate,G,from which the separate values ofKI,KIIandKIIIare calculated(Nakamura and Parks,1989).In the case of closely spaced fracture fronts,the displacement correlation method is used to compute the SIFs (Nejati et al.,2015b).

TheI-integral, solved by the disk-shaped domain integral method, is formulated from the concept of two states of equilibrium for a cracked body. The FEM solution (Eq. (1)) is the actual state, with displacement (u), strain (e) and stress (σ) values throughout the domain. The actual state will be accurate in the domain and close to the fracture,but very close to the tip,it will not reproduce the near tip stress singularity due to the limits of the discretisation.A second state,the auxiliary state,is given by known asymptotic fields, which are functions of the SIFs, and has its own forms of the displacement, strain and stress fields (uaux,eaux, and σaux).These values are found from plane strain fields,which prevail close to the tip and are relatively independent of the far-field conditions (Nakamura and Parks,1989; Nejati et al., 2015a). Each field will have its own energy release rate,and thereforeJ-integral value, known as the actual (Jact) and auxiliary (Jaux) fields. The superimposedJ-integral is the combination ofJact,JauxandI:

wheresis a point near the fracture tip, andIis the interaction integral,relating the actual and auxiliary states.Eq.(7)can be written in terms of SIFs as

whereKi(i=I，II，III)is the SIF for modei= I， II， III andKauxis the SIF for modeidue to the auxiliary state,E′=E/(1- ν2)(whereEis the Young’s Modulus and ν is the Poisson’s ratio), and μ is the shear modulus of the rock.

To compute theI-integral on an unstructured mesh, a 2D diskshaped surface is resampled with a virtual mesh. Fig. 2a shows a schematic of this disk at a fracture tip, whereAis the integration domain andRdis the domain radius of the integration.The terms Γ,Γ0, Γ+, and Γ-are the boundaries enclosingA, with Γ being the point at the centreA,Γ0being the outer boundary ofA,and Γ+and Γ-as the top and bottom edges whereAmeets the fracture surfaces, respectively. The virtual 2D elements which discretise this disk are shown in Fig.2b,along with the local tip coordinate system withx1- (along the fracture plane) andx2-axis (normal to the fracture plane, see Fig. 2b). Line elements are positioned whereAtouches the fracture plane(on Γ+and Γ-).This mesh is virtual as it is discarded after solving theI-integral and does not need to conform to the volumetric mesh.

TheI-integral is evaluated over domainAusing the following formulation (Nejati et al., 2015a):

where σij,eij, anduiare the Cartesian components of the stress tensor, strain tensor and displacement vector in thex1,x2coordinate system,respectivelyare the components of the stress tensor,and strain tensor,and displacement vector which are given by the auxiliary field,respectively;WIis the mutual strain energy densityδijis the Kronecker delta;qis a smooth scalar function in the region enclosed by Γc= Γ0+Γ++ Γ-- Γ,equalling unity on Γ and vanishing on Γ0;andmis a vector normal to the plane ofA,such thatm= -1 on Γ+.The two terms in Eq.(9)are solved over virtual triangular and line elements,respectively, and then are combined. Use of the smooth scalar functionqavoids the particular form of the contour integrals in the originalJ-integral;therefore,the last term in Eq.(9)is transformed into an area integral overA, together with contour integrals over the fracture surface. This is both more convenient and more accurate in a numerical implementation (Nikishkov and Atluri, 1987).The full numerical implementation of Eq.(9)is given in Algorithm 1 in Nejati et al.(2015a).Three auxiliary states are used,one for each deformation mode(pure modes I,II and III),producing three forms of Eq. (9) (II(s),III(s) andIIII(s). The SIFs can then be related toIby combining Eqs. (7) and (8).

In each pure mode auxiliary state, the other SIF terms vanish,thus one expression can be formed per SIF:

2.3. Fracture nucleation

Microscopic defects within the rock lead to the formation of a‘damaged zone’,where a local reduction in the rock strength occurs(Boukharouba et al., 2009). A damage model is appropriate to describe fracture nucleation,since the process occurs at microscale,where the rock is loaded beyond its elastic limit.Fractures are set to nucleate in the elements having the largest damage, where the material stiffness has been degraded the most due to the applied stresses.The model of Mazars(1984)is used to capture the process of stiffness degradation, as it has been shown to predict quasibrittle behaviour accurately:

Fig.2. Diagram showing the disk-shaped domain integral used to calculate stress intensity factors and strain energy release rates along fracture tips.(a)Disk-shaped domain shown in the context of a fracture surface.The terms Γ,Γ0,Γ+,and Γ- show the edges of the integration domain A.Rd is the domain radius.(b)Discretisation of one quadrant of the diskshaped domain, used to resample the stress and displacement fields. The dotted line shows the location of nodes and mid-side nodes along the edges of the elements in the quadrant, highlighting the quarter point element adjacent to the tip (modified after Thomas (2019)).

The damage variable,d, is found by calculating an equivalent strain, ~ε, which is the L2norm of the positive part of the strain tensor:

For a general stress state,the damage parameter is defined as a linear combination of damage under tension and under compression, if the equivalent strain is larger than the strain threshold:

where εtiand εciare the principal strains due to positive and negative stresses,respectively;the exponentb,a damping factor,is equal to 1.06 (Jirasek, 2011); ε0is the strain threshold, defined as the ratio of undamaged Young’s modulus and rock tensile strength(ft);andAt,Bt,Ac,andBcare material parameters that modulate the shape of the uniaxial stress—strain curves in tension and compression, respectively. In elements where equivalent strains that have not yet exceeded the threshold,the damage variable is set to zero and no fractures nucleate.

2.4. Fracture growth

Fracture growth is determined at each fracture tip using the calculated SIFs andGas inputs to three growth criteria. First, the failure criterion(Section 2.4.1)determines whether a particular tip extends. If extension occurs, a vector is generated using the angle(Section 2.4.2)and length criteria(Section 2.4.3).New surfaces are added between growing fracture tip node pairs, based on the generated vectors.At the start of the next growth step,the domain is re-meshed,with new fracture tip nodes placed at the end of the growth vectors (Thomas et al., 2020).

2.4.1. Mixed mode propagation criterion

Fracture propagation can occur due to pure modes I, II, or III deformation, or a combination of all three modes. A mixed mode failure criterion is therefore used, accounting for tension,compression,and shear simultaneously(Richard et al.,2014).Using the criterion, fracture tips grow based on a comparative SIF (Kν),derived from the values ofKI,KIIandKIII. IfKν exceeds the local fracture toughness of the material, fracture extension occurs(Richard et al., 2005). Richard et al. (2014) found that the most accurate condition for mixed-mode growth is

whereR1andR2are the ratios of the mode I fracture toughness to the modes II and III fracture toughness(R1=KIC/KIIC,R2=KIC/KIIIC,where subscript C refers to the critical SIF). For brittle rocks,KIICtends to be higher thanKIC,andKIIICtends to be lower.In this work,KICvaries between different numerical experiments,KIC/KIIC= 1/1.3268 andKIC/KIIIC= 1/0.6297 (Chang et al., 2002; Aliha et al.,2015; Wei et al., 2016).

2.4.2. Propagation angle criterion

The mixed mode criterion of Richard et al.(2005)finds the angle of a fracture tip’s propagation vector as a function of the SIFs. It calculates two angles perpendicular and tangent to the crack front,and has been experimentally validated (Richard et al., 2013). The perpendicular angle(φ),or fracture deflection angle along its plane is calculated from the three SIFs as follows:

whereX=140°andY=-70°.Angle φ is negative whenKII>0 and positive whenKII＜ 0 andKI≥ 0.For simplicity,the tangent angle is not incorporated into the growth of the presented fracture networks.

2.4.3. Propagation length criterion

Extension lengths of quasi-statically propagating fractures can be predicted using Paris-type propagation criteria (Paris and Erdogan, 1963), which scale extension based on the value ofGalong fracture tips (Lazarus, 2003; Renshaw and Pollard, 1994;Pugno et al., 2006). Thomas et al. (2020) introduced a modified Paris-type approach for handling propagation in domains containing 3D multiple fractures, motivated by the need to scale the relative propagation rate across many fractures in order to reproduce growth patterns.This two-step approach assigns each fracturefwith a maximum extension length in the first step(Δlf),and each tip nodenwith an extension length in the second step (Δln). The distribution of extension lengths between 0 and Δlνis controlled by two growth exponents,βfand βn,used in the first and second step,respectively. The first step is performed once per fracture, as

whereGfis the largest value ofGon fracturef,andGνis the largest value ofGin the domain. The second step is calculated once per fracture tipn:

whereGnis the value ofGat tipn.

2.4.4. Validation

The ICGT simulator, which can model a range of physical processes, has been widely validated during its development. Where possible, this validation has been performed using LEFM-based analytical solutions of fracture properties. A representative chronological list of published validations of the methodology is as follows:

(1) SIF computation:validated by Nejati et al.(2015b)for quarter points and displacement correlation,Nejati et al.(2015a)for the disk-shaped domain integral method, and for all combinations of implemented computation methods and quarter points, on meshes optimised for growth (Thomas et al.,2020).

(2) Fracture growth: validated by Paluszny and Matth?i (2009)in two dimensions on interacting fractures in PMMA, and in three dimensions for isolated fractures under tension and compression in rock (Thomas et al., 2020). The growth and aperture evolution of hydraulic fractures is also validated by Salimzadeh et al. (2017b).

(3) Friction and contact:validated in Nejati et al.(2016)in three dimensions and through the validation of SIFs.

(4) Fracture interaction: validated in Thomas et al. (2017) for static coplanar fractures, as a function of the change in SIFs along a fracture’s tip. The mutual growth of interacting fracture arrays is also analysed in Thomas et al.(2020)and in Section 3.1.

(5) Flow and permeability: validated by Paluszny and Matthai(2010) in two dimensions and Lang et al. (2014) in three dimensions, comparing equivalent permeability with analytical solutions for sugar cube models and parallel fractures.

3. Results and discussion

In this paper, several numerical experiments are performed using the methods described in Section 2 and implemented in the ICGT. These numerical experiments demonstrate how fracture networks nucleate, grow, interact, and develop distinct flow channels. Exploring these phenomena numerically, with confidence that they produce realistic results, is possible through the separate validation of each of the physical processes being modelled, and their implementation together in the ICGT.

The first set of results presents the growth and characterisation of a five-fracture array (Fig. 3). This is followed by the analysis of apertures and channelling (Figs. 4 and 5) and within a larger geomechanically grown network, consisting of 90 growing fractures.Two datasets demonstrate how a network can form when fractures initiate using the nucleation approach, as described in Section 2.3.The first of these is limited to 100 fractures, which nucleate and grow forming a pattern dominated by few large fractures (Fig.1).The second has a lower maximum number of fractures,resulting in a smaller range of fracture sizes (Fig. 7). The methodology for solving for flow in fractured domains, as described Section 2.2, is also demonstrated for this network, showing fluids migrating through the system when apertures are determined by the separation between fracture surfaces. Table 1 summarises the material properties and boundary conditions used in each numerical experiment in the present work.

3.1. Geomechanically generated fracture array

Fig.3. Grown fracture array comprised of five fractures at the centre of a 20 m×20 m×20 m domain.(a)and(b)show perspective and plan views,respectively,visualising fracture surfaces as opaque planes.Four cut-planes(numbered 1 to 4)are shown as transparent planes,and their intersection with the fractures is shown in(c)—(f)as ‘virtual outcrop trace maps’.The simulation uses E=20 GPa,ν=0.213,KIC =106 Pa m1/2,βn =2,βf =0.35.The boundary condition is a tensile stress on the+z boundary of 8 × 107 Pa in the+z-direction,with the -z boundary being fixed (no displacement).

Fig.4. Aperture distribution on fractures in a geomechanically grown network containing 90 initial fractures.The P32 fracture intensities are 0.044 m-1(step 0),0.352 m-1(step 9),0.396 m-1 (step 12), and 0.426 m-1 (step 17). Fractures are visualised as opaque planes coloured according to local variations in aperture.

Fig.5. Flow channelling through a fracture network.Aperture distributions and additional details are given in Fig.4.Fractures are visualised as transparent planes.Arrows indicate flow direction in each element, with colour and size indicating magnitude. Flow occurs in the rock matrix at a much lower velocity; therefore vectors are only visible on the fractures. Flow is from the -x boundary to the +x boundary, with a fluid pressure gradient of 100 Pa and 10 Pa, respectively.

Fig.6. Stress intensity factors of all fracture tips during nucleation and growth in the domain shown in Fig.1.(a)KI,(b)KII values greater than zero,and(c)KII values less than zero.These data are presented as box plots, where horizontal lines represent the median, and rectangles represent the range between the upper (q3) and lower (q1) quartiles. Vertical lines end at the last SIF which is not an outlier, defined as being either above q3 + 1.5 (q3 - q1) or below q1 -1.5 (q3 - q1), respectively, and crosses represent outliers. The distributions for KIII closely resemble KII.

Fig. 3a and b shows the 3D geometry of five fractures, grown geomechanically in tension in an array over 36 steps, demonstrating mechanical fracture interaction independently from hydraulic effects. The five fractures are initially circular with 1-m radius. The middle fracture is centred at (0, 0, 0) and the other fracture centres have a 1-m separation in both thex-andz-axis.They are oriented perpendicular to the direction of tensile stress,such that if they are unaffected by mechanical interaction, they would maintain a circular shape. However, their proximity results in mutual interaction which changes the SIFs along the fracture tips,with the largest change occurring where the fractures overlap.This interaction reduces the growth of overlapped tips and makes the growth mixed-mode, resulting in hooking towards the interacting surface. No additional fractures are nucleated during the simulation.Therefore,the only influences on growth are the tensile stress boundary condition and interaction between existing fractures.

Fig. 3c—f shows the intersection between four cut-planes and the array, simulating how the fractures may appear in an outcrop trace map. The four resulting sets of fracture curves are broadly similar, but demonstrate the variety of curvatures which can be observed depending on the inclination of the visualising cut-plane.Fig.3c and e is the most alike,demonstrating fracture reorientation occurring in the overlapped zones, particularly on the fractures in the middle of the array. The apparent curvature of the fractures is different from that in Fig.3d,whose plane is offset from the central line where the initial flaws are placed. Additionally, the top and bottom fractures now curve away from, rather than towards, the middle fractures.The plane which creates Fig.3f is inclined,causing the apparent separation between the fractures to be more significant.These small differences could lead to different interpretations of fracture behaviour in outcrops or experiments, despite all the cut-planes originating from the same 3D array.

3.2. Apertures and channelling

Fracture apertures, which control flow and macro-scale permeability, respond to the deformation of fracture surface,resulting in variation of flow within each individual fracture as well as across the network.In the coupled processes,the actual growth and extension of fracture surfaces during interaction may result in an additional controlling mechanism for flow and transport of solute through a network.Figs.4 and 5 demonstrate the influence of fracture growth and interaction on aperture and flow distributions for a geomechanically growing network. The domain contains 90 fractures which grow in response to a tensile stress at the top(+z)boundary with magnitude 2 × 106Pa, and the bottom (-z)boundary is fixed. The elastic rock matrix has material properties matching crystalline rock(E=20×109Pa,ν=0.2,KIC=106Pa m1/2,Km= 1×10-15m2).In the initial domain, 90 fractures,discretised with 20 tip nodes,are randomly inserted with centres in the central 14 m×14 m×14 m portion of the domain.Seventy-five fractures have a 1-m radius and are oriented with normal vectors to their surfaces of(0,0,1),representing flaws which grow immediately in response to the tensile stress regime.Fifteen fractures have a 2.5-m radius and normal (4, 5, 1), representing a previous phase of deformation which is initially inactive, i.e. they do not meet the failure criterion until fracture interaction causes them to become critically stressed. Fracture growth parameters are Δlν= 2m,βn= 2,and βf= 0.35, with no additional fractures nucleating.

In the first few growth steps, the fracture network grows quickly, with horizontal fractures being activated first. In step 4,percolation occurs between the horizontal axes due to fracture and boundary intersections. Once intersected, fracture tips deactivate,resulting in gradually less extension per step. The steps visualised in Fig. 5 show the network in a mature state, where significant anisotropy has been developed due to growth perpendicular to the orientation of tensile stress. The orientation of stress also means that apertures are highest at the centre of horizontally oriented fractures and very low on the initially inclined fractures. Mechanical fracture interaction also influences apertures, as fractures which are less overlapped towards the edges of the domain (i.e.outside of the stress shielding region of more fractures)also acquire higher apertures.

The three growth steps shown in Fig. 5 demonstrate multiple percolating flow paths, whose locations are determined by evolving fracture connectivity and aperture distributions throughout the network. Fractures with large apertures and geometrically small intersections have the highest flow velocities.The flow path with the highest magnitude of flow changes between steps 9,12 and 17. In step 9, most of the flow occurs near the+yboundary.This flow path is presented in step 12,along with a second path near the -yboundary. In step 17, the second flow path becomes the primary flow channel. The fixed bottom boundary causes apertures to gradually increase from the -zto the +zboundaries, therefore, most of the flow paths are found near the stressed boundary.

Fig.7. Aperture(m)and flow velocities(m/s)on a fracture network.In each growth step,isometric and plan views of the domain are shown.Surfaces are visualised as transparent planes with colour based on the local fracture aperture. Flow vectors are visualised as arrows, with magnitude being represented by size and colour.

Table 1Material properties and boundary conditions of the numerical experiments shown in this work in Figs.1 and 3—5 and 7.

For flow in a geomechanical network with uniform fracture apertures, the simulation results yield a different qualitative distribution of flow velocities. When assuming uniform fracture permeability, velocities are uniformly distributed, with higher values on the fractures, and without the formation of channels as depicted in Fig. 5. In particular, when a uniform aperture distribution is assumed, the inclined fractures acquire relatively higher flow velocities, resulting in much less channelling along horizontal fractures, as observed here. Flow channels are sensitive to slight changes in aperture distributions (Tsang and Neretnieks,1998). In addition to channel changes due to the redistribution of loosened filling material,and changes to the compressive stresses,this work provides some evidence that channel redistribution may occur through fracture growth, inducing small changes to the interactions between fractures and connectivity of the network.

3.3. Hydro-mechanical interaction effects during growth

Two datasets are generated to investigate the flow properties of a fracture network,formed from nucleated fractures as described in Section 2.3, rather than randomly placed flaws. As opposed to the previously simulated cases, fractures nucleate during the simulation in this dataset. The rock matrix is assumed to be cuboid,measuring 40 m × 40 m × 20 m. The first fracture network, as shown in Fig.1,contains 100 fractures at its final stage(grown with Δlν= 2 m, with initially disc-shaped fractures of radius 0.5 m),while the second dataset contains 20 fractures at its final stage(grown with Δlν= 1 m, with initially disc-shaped fractures of radius 1 m). In both cases, βn=1 and βf= 0.35. In the first case,fewer initial fractures dominate the final pattern, while in the second case,more initial fractures grow throughout the simulation.

In Fig. 1, the model is initially intact and remote principal extensional horizontal stresses of σx= σy= 1 MPa, and vertical compressive stress σz=-0.3 MPa,are applied(Case 3.3 A),where a negative stress indicates compression. The rock is assumed to be linear elastic, homogeneous, and isotropic granite. Thus, the domain undergoes bilateral extension, driving nucleation and growth of fractures,with growth competing for alignment with thex-andy-axis, and strong interaction as a result. Further material properties are as follows:E= 1 GPa,ft= 0.25 MPa, friction coefficient of 0.6, andKICis set intentionally low at 0.1 MPa m1/2, to encourage nucleated fractures to grow.Each fracture tip is initially discretised by 20 tip nodes. This growth simulation aims to benchmark the growth behaviour for a case in which fracture growth displays strong interaction.

During simulation, the refinement of tip fractures is kept constant, resulting in refined fronts as fracture growth progresses. As shown in Table 2, the number of fractures and nodes increases during the simulation, ranging from 3031 nodes to up to 107,179 nodes. Fracture radius is approximated based on the total surface area, assuming that the shape of the fracture is approximately a disc. As fractures intersect the boundaries, the regions of the fracture requiring refinement reduce, and the total amount of nodes further reduces.Once fractures saturate the medium,the number of nodes stabilises. Furthermore, once fractures intersect the boundaries and fewer larger fractures are dominant, the number of growing nodes decreases. The dynamic evolution of the mesh and its adaptation to the evolving geometry allow the nodes to be redistributed to emphasise detail in the tips of the fractures.

Table 2Progression of geometric data during fracture growth in the domain shown in Fig.1.For each growth step, the number of fractures, number of nodes, number of elements(tetrahedra and triangles),and largest fracture radius at each step are shown.The total degree of freedom equals the number of nodes times three for mechanicsonly simulations,and times five for hydro-mechanical simulations.Similar numbers of nodes and elements are used in the other numerical experiments presented in this work.

Growth is concentrated on fracture tip nodes which have the highest strain energy release rates(G),as described in Section 2.4.3,and in broadly tensile regimes,KIcorrelates with fracture size.Thus,most of the growth in the domain is on a few large fractures,particularly those nucleated in the domain at earlier growth steps.This is consistent with observations of rock and material failure,where initial stress concentration points continue to concentrate stress during the progressive failure of a sample.Fractures generally maintain a straight growth path,with exceptions for fractures that propagate close to others.In these cases,slight curvature develops where fractures hook towards one another,as observed in the array in Fig. 3. Hooking does not always occur as a result of interaction.Boundary stresses also influence fracture curvature,as fracture SIFs are modified by the presence of a free surface (e.g. Pollard and Holzhausen,1979; Al-Ostaz et al.,1998). Stress concentration also influences fracture nucleation, as new fractures tend to be added near boundaries and near the tips of fractures, in ‘process zones’,where stress is concentrated (Thomas et al., 2017).

Fig. 6 shows the change of SIFs in modes I and II during the nucleation and growth of fractures in the domain. The number of data points increases gradually as more fractures nucleate. The trend in mode I (KI) captures how opening mode deformation gradually increases over time,whereas modes II and III increase at a slower rate.These graphs show that the fractures are subjected to a complex multi-modal stress regime, resulting in the local reorientation of the fracture to grow in extension. SIF magnitudes are spread over several orders of magnitude due to the large variations in fracture sizes and orientations emerging during the process of nucleation and growth.

Fig.7 shows the evolution of the second pattern under the same conditions, where the maximum number of nucleated fractures is set to 20 rather than 100.Top and side views show the geometries of the final stages of growth, before the end of the simulation.Relative to Fig.1, fractures become larger with more even sizes, as there is less stress shadowing.As expected,at early stage of growth,the matrix evenly distributes flow amongst the fractures,and flow velocity is highest at the centre of fractures. As growth and interaction progresses, this relationship becomes more complex, as depicted in Fig. 7. Flow does not merely follow the aperture distribution,and highest flow velocities are accrued in a sub-region of a subset of fractures. This is consistent with results in Section 3.2.The steps shown depict a fracture network at growth steps 14,18,19 and 21. As fractures interconnect, preferential pathways form, and channels become apparent. In this example, channelling is significantly controlled by the geometric intersections between the fractures in the domain and is enhanced as growth advances. In step 14, the most flow occurs at a single intersection of two fractures near the-zboundary.Due to growth,additional intersections have occurred in step 18,creating multiple regions with high flow rates. In steps 19 and 21, many intersections have occurred, connecting the fracture network and providing long channels for high velocity flow.

At each step of growth shown in Fig. 7, the network’s hydromechanical properties can be quantified using the effective permeability tensor. In particular, the variation of the equivalent permeability is examined as a function of growth. Matrix permeability is defined as 1×10-15m2and porosity is 15%.Thus,flow is modelled through both matrix and fractures simultaneously. Apertures are mechanically determined for each node on the surface of each fracture, as a result of applying a load of 0.1 MPa, with extension in they-direction and compression in thex-andz-direction. The largest apertures range from 1.0 × 10-4m to 1.1×10-4m.Table 3 shows the variation of permeability along thex-,y-, andz-direction at late stages of growth, highlighting the strong variations occurring during the formation of the channels.The emergence of flow channels disrupts fluctuations which do not necessarily follow the monotonic increase of density.

Table 3Permeability of the fractured network visualised in Fig.1 in the x-,y-and z-direction.

All simulations run in a 44-core logical processor shared memory machine with 32 GB of RAM. Simulations take between minutes and hours to generate the patterns presented in this paper,depending on the number of nodes in the domain and the degrees of freedom associated with each simulation (see Table 2).

4. Conclusions

This work examines the development of channelised flow in numerically modelled, geomechanically-generated discrete fracture networks. A finite element based LEFM model is described in detail, including the methods used to quantify fracture tip stress states, propagate fractures, and compute flow and network permeability.The methods are implemented in the ICGT,which has been validated against a range of important analytical solutions.The numerically generated networks contain many non-planar fractures which form as a result of remote stresses and local interactions in the network. Mechanical fracture interaction affects the evolution of the networks, controlling the relative positions where fractures are nucleated, and determining which fractures become large and intersect with one another. Furthermore, interaction affects aperture distributions and the effective permeability of networks.For complex networks,preferential flow patterns form as a function of both topology and stresses.These fracture networks exhibit strong preferential flow pathways that are more pronounced at specific sub-regions of the interconnected fracture geometry, resulting from stronger interaction behaviour during growth. For growing fracture patterns, channels become more pronounced during growth as fractures become intersected and wider apertures form in favourably stressed regions. Results show that channels can form without explicitly defined fracture surface heterogeneity, which if modelled, would further focus flow into distinct channels. Patterns show that self-organisation during fracture growth has an effect on aperture distribution and subsequently on flow, an effect that results from intersection between fractures before and during intersection.

Declaration of Competing Interest

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

Acknowledgments

The authors thank Yvonne Tsang and Chin-Fu Tsang for valuable conversations related to this paper. The authors thank the Natural Environment Research Council for the funding received for project NE/R018065/1.The authors thank the Royal Society for funding this research through fellowship UF160443.The authors also thank the Swedish Nuclear Fuel and Waste Management Company (SKB) for partially funding this research.