The use of discrete fracture networks for modelling coupled geomechanical and hydrological behaviour of fractured rocks

We present a discussion of the state-of-the-art on the use of discrete fracture networks (DFNs) for modelling geometrical characteristics, geomechanical evolution and hydromechanical (HM) behaviour of natural fracture networks in rock. The DFN models considered include those based on geological mapping, stochastic generation and geomechanical simulation. Different types of continuum, discontinuum and hybrid geomechanical models that integrate DFN information are summarised. Numerical studies aiming at investigating geomechanical effects on fluid flow in DFNs are reviewed. The paper finally provides recommendations for advancing the modelling of coupled HM processes in fractured rocks through more physically-based DFN generation and geomechanical simulation. 2017 The Authors. Published by Elsevier Ltd. This is an openaccess article under the CCBY license (http:// creativecommons.org/licenses/by/4.0/).


Introduction
Fractures such as joints, faults, veins and bedding planes are ubiquitous in crustal rocks. These naturally occurring discontinuities often comprise complex networks and dominate the geomechanical and hydrological behaviour of subsurface rocks [1]. Understanding of the nontrivial effect of fractures is a challenging issue which is relevant to many engineering applications such as underground construction, enhanced geothermal systems, unconventional shale gas production (fracking), groundwater management and radioactive waste disposal [2][3][4]. The importance of the presence of natural fractures, which can result in heterogeneous stress fields [5] and channelized fluid flow pathways [6] in highly disordered geological formations, has promoted the development of robust fracture network models for numerical simulation of fractured rocks [7].
The first difficulty in modelling fractured rocks is the geometrical representation of complex three-dimensional (3D) discontinuity systems. Natural fractures form under certain mechanically self-organised dynamics, where breakage and fragmentation can occur at all scales [8]. They are subject to in-situ stress fields at depth and can form intricate topologies, such as cross-cutting, abutting, branching, termination, bends, spacing and clustering [9,10]. However, direct observation of the detailed 3D structure of fracture networks deep in the crust is impossible. Field data are usually collected from lower dimensional limited exposures, e.g. one-dimensional (1D) borehole logging and two-dimensional (2D) outcrop mapping [11]. Seismological surveys may be able to locate 3D large-scale structures but the current technology can hardly detect widely-spreading medium and small fractures due to the resolution limit. The description of natural fracture geometries, therefore, has to largely rely on extrapolations, from 1D/2D to 3D and from small samples to the whole study domain. Hence, the question of how to create realistic fracture networks remains an unresolved issue. Often simplifications have to be made by ignoring some details of less importance for the particular questions at hand.
The second fundamental issue in modelling fractured rocks is to solve the fracture and solid mechanics problems of the discontinuous geological media under complex boundary conditions, i.e. the fractured rock mass mechanical response to boundary stresses or displacements. When the spacing of natural fractures is comparable to the scale of interest of the problem to be modelled, the conventional continuum approach may not be adequate to capture some important mechanical behaviour of fractured rocks [12], such as fracturing of intact rocks [13], interaction of multiple fractures [5], and opening, shearing and dilatancy of rough fracture walls [14]. Thus, many computational schemes based on extended continuum, discontinuum, or hybrid continuum-discontinuum methods have been developed to solve a numerical system with fracture geometries explicitly represented.
Another important question is about the hydrological behaviour of fractured rocks under geomechanical conditions. Fractured rocks may deform in response to geological and/or man-made perturbations (e.g. tectonic events, underground excavations), resulting in changes of bulk permeability and fluid migration [2]. The intricacy of such coupled hydromechanical (HM) processes is increased if the presence of natural fractures associated with topological complexities is to be considered [1,15]. The quest for a means of quantifying the influence of in-situ stresses on the permeability of fractured reservoirs has been driven largely by the motivation from petroleum engineering [16]. The understanding of contaminant migration through tectonically strained fractured formations is also crucial for the groundwater community [17] and nuclear waste management [3,4].
The objective of this paper is to review the current state-of-theart of fracture network models and to provide some discussions and recommendations for HM modelling of fractured rocks. In this context, the issues arising in the aforementioned three key subject areas (i.e. geometry, geomechanics and hydromechanics) will be examined progressively. The paper will present a summary of various approaches used in developing fracture network models that represent natural fracture geometries often with different degrees of simplification, and different numerical frameworks that integrate discrete fracture representations for modelling geomechanical and HM behaviour of fractured rocks. One important focus of this work is to discuss the features or characteristics of a fracture network model needed to properly simulate coupled HM processes in fractured rocks. The rest of the paper is organised as follows. Section 2 reviews the methods of representing natural fracture geometries by geological mapping, stochastic generation or geomechanical simulation. Section 3 provides a brief overview of continuum and discontinuum models that integrate explicit fracture geometries for geomechanical modelling of fractured rocks. A short summary is presented in each of Sections 2 and 3 with a view towards HM modelling. Section 4 summarises the numerical studies of geomechanical effects on fluid flow in fractured rocks and further provides some suggestions on choosing appropriate DFNs and geomechanical models for HM modelling. Finally, outstanding issues are discussed and some concluding remarks are made. It has to be mentioned that this review is focused on HM behaviour of fractured rocks and is not exhaustive. Only limited references are cited for brevity. For more extensive descriptions of fracture network models, the reader is referred to the reviews by Dershowitz and Einstein [18] and Liu et al. [19], and the textbooks by Adler and Thovert [20] and Adler et al. [21]. Extensive reviews of different numerical methods developed in the field of rock mechanics can be found in Jing and Hudson [22], Jing [12], Yuan and Harrison [23], Bobet et al. [24], and Lisjak and Grasselli [25]. More in-depth discussions about stress effects on fluid flow in fractured rocks can be found in Zhang and Sanderson [26], and Rutqvist and Stephansson [2].

Representation of natural fracture networks
A ''discrete fracture network" (DFN) refers to a computational model that explicitly represents the geometrical properties of each individual fracture (e.g. orientation, size, position, shape and aperture), and the topological relationships between individual fractures and fracture sets. Unlike the conventional definition of DFNs that corresponds to stochastic fracture networks, the term DFN here represents a much broader concept of any explicit fracture network model. A DFN can be generated from geological mapping, stochastic realisation or geomechanical simulation to represent different types of rock fractures including joints, faults, veins and bedding planes.

Geologically-mapped fracture networks
Fracture patterns can be mapped from the exposure of rock outcrops or man-made excavations (e.g. borehole, quarry, tunnel and roadcut). These geologically-mapped fracture networks were widely used to understand the process of fracture formation [5,27], interpret the history of tectonic stresses [28][29][30], and derive the statistics and scaling of fracture populations [8,31,32]. However, digitised outcrop analogues (Fig. 1) can also be used to build DFNs for numerical simulations. For example, a series of discrete fracture patterns were mapped from limestone outcrops at the south margin of the Bristol Channel Basin, UK [33]. The traced DFNs were used to study the connectivity [34], multiphase flow [35][36][37], solute transport [38] and HM behaviour [15,[39][40][41] of natural fracture systems. Similar outcrop-based DFNs have also been constructed by many other researchers for modelling natural fracture systems [42][43][44][45][46][47][48]. It is noteworthy that the terrestrial remote sensing techniques (e.g. LIDAR and photogrammetry) have recently been introduced to better capture large-scale trace features [49][50][51][52]. Fracture apertures may be determined from a detailed field mapping [53,54] and further calibrated by comparing flow simulation results with in-situ measurements [55]. However, apertures were more commonly assumed to be constant or to follow an a priori statistical distribution (sometimes correlated with trace lengths).
Advantages of such an outcrop-based DFN model include preservation of natural fracture features (e.g. curvature and segmentation) and unbiased characterisation of complex topologies (e.g. intersection, truncation, arrest, spacing, clustering and hierarchy). However, it is often constrained to 2D analysis and may not be applicable for general 3D problems involving obliquely dipping fractures. The LIDAR survey combined with assumptions of fracture shape/persistence [56,57] may provide an estimation of 3D structural variations in near surface, but it is not applicable to buried geological complexities. Extrapolation from borehole imaging can provide an estimation of 3D fracture distributions but confidence can only be guaranteed for the areas close to boreholes [58]. Seismic data can be used to build 3D maps of large-scale geological structures, for which, however, the limited resolution often obscures detailed features such as the segmentation of faults and impedes the detection of small cracks widely spreading in subsurface rocks [59].

Stochastically-generated fracture networks
Due to the difficulty of performing a complete measurement of 3D natural fracture systems, stochastic approaches using statistics from limited sampling have been developed and widely used [18]. The stochastic DFN method emerged in the 1980s with the aims to study the percolation of finite-sized fracture populations [60][61][62] and fluid flow in complex fracture networks [63][64][65][66][67].
The general stochastic DFN approach assumes fractures to be straight lines (in 2D) or planar discs/polygons (in 3D), and treats the other geometrical properties (e.g. position, frequency, size, orientation, aperture) as independent random variables obeying certain probability distributions [68] derived from field measurements such as scanline [69] or window sampling [70] of outcrop traces as well as borehole imaging [71]. The orientation data can be processed using a rosette or stereogram so that fractures can be grouped into different sets with their orientations characterised by e.g. a uniform, normal or Fisher distribution [11]. Fracture sizes may exhibit a negative exponential, lognormal, gamma or power law distribution [72,73]. Fracture frequency can be described in terms of fracture density or fracture intensity and denoted using the P ij system, where i gives the dimension of the sample and j the dimension of the measurement [74]. Fracture density measures the number of fractures per unit volume (P 30 ), area (P 20 ) or length (P 10 ), while fracture intensity represents the total fracture persistence per unit volume (P 32 ), area (P 21 ) or length (P 10 ). Fracture spacing is related to fracture frequency (P 10 ) [75], which is both a density and intensity measure. Fracture spacing may follow a negative exponential, lognormal or normal distribution depending on the degree of fracture saturation in the network [76]. Fracture apertures usually obey a lognormal [77] or power law distribution [78,79], and may be related to fracture sizes by a power law [73] with a linear [5] or sublinear [80] scaling relationship. The measurement data (e.g. those based on outcrop sampling) may be biased under the truncation and censoring effects and requires to be amended to determine the underlying statistical distributions [81][82][83]. The statistics of 3D fracture systems (e.g. size, spacing, density distributions) may be extrapolated from 1D/2D sampling data based on stereological analysis [84,85]. In the stochastic simulation, fractures can be assumed randomly located (represented by their barycentres), while the geometrical attributes can be sampled from the corresponding probability density functions [18]. Such a random fracture network modelling approach, termed the conventional Poisson DFN model (or Baecher model) (Fig. 2), has been implemented within the commercial software FracMan [86] and also adopted by many research codes to study the connectivity, fragmentation, deformability, permeability and transport properties of fractured rocks in the past three decades  (only a few among many others).
However, the conventional Poisson DFN model tends to have large uncertainties due to its assumption of a homogeneous spatial distribution, simplification of fracture shape using linear/planar geometries, and negligence of the correlations between different geometrical properties as well as disregard of the diverse topological relations (e.g. ''T" type intersections). Several researchers have examined the conventional Poisson DFN model by comparing it with an original natural fracture network with respect to geometrical, geomechanical and hydrological properties, and significant discrepancies were observed [36,39,[108][109][110]. Several improvements on the Poisson DFN model have been developed and include considerations of: (i) the inhomogeneity of fracture spatial distribution based on a geostatistically-derived density field [67,111] or a spatial clustering process [112][113][114], (ii) the correlation between fracture attributes (i.e. length, orientation and position) based on an elastic energy criterion [115,116], (iii) the unbroken areas inside individual fracture planes based on a Poisson line tessellation and zone marking process [18,117], (iv) the topological complexity based on a characterisation of fracture intersection types [118] and the termination of fractures at intersections with pre-existing fractures [86], and (v) the mechanical interaction between neighbouring fractures based on a stress shadow zone model [76,119,120]. It has to be mentioned that an alternative to the Poisson DFN model may use the spacing distribution to locate fractures in the stochastic generation [121], which is however considered more suitable for highly persistent fracture systems. For more detailed summary of the Poisson model and other improved models, the reader is referred to Dershowitz and Einstein [18], Staub et al. [122] and Fadakar-Alghalandis [114].
A more systematic characterisation of the hierarchy, clustering and scaling of natural fracture systems may involve the methods of fractal geometry and power law models [73]. Extensive field observations suggest that fracturing occurs at all scales in the crust and creates a hierarchical structure that can exhibit long-range correlations from macroscale frameworks to microscale fabrics [8,124,125]. The spatial organisation of natural fracture networks can be characterised by the fractal dimension D, which quantifies the manner whereby fractals cluster and spread in the Euclidean space and can be measured using the box-counting method [31,126,127] or the two-point correlation function [128,129]. The density and length distribution of a fracture population can be then described by a statistical model [32,130] given as: n(l, L) = aL D l Àa , for l 2 ½l min ; l max , where n(l, L)dl gives the number of fractures with sizes belonging to the interval [l, l + dl] (dl ( l) in an elementary volume of characteristic size L, a is the power law length exponent, a is the density term, and l min and l max are the smallest and largest fracture sizes. The exponent a, which defines the respective proportion of large and small fractures, can be derived from the cumulative distribution or density distribution of fracture lengths [72,82]. In theory, D is restricted to the range of [1,2] for 2D and [2,3] for 3D, and a is restricted to [1,1] for 2D and [2, 1] for 3D. However, extensive measurements based on 2D trace maps reveal that generally D varies between 1.5 and 2 and a falls between 1.3 and 3.5 [73] in natural fracture systems. The D and a measured from 1D/2D samples can be extrapolated to derive 3D parameters based on stereological relationships [123]. The density term a is related to the total number of fractures in the system and varies as a function of fracture orientations [130]. The extent of the power law relation is bounded by an upper limit l max that is probably related to the thickness of the crust, and a lower limit l min that is constrained by a physical length scale (e.g. grain size) or the resolution of measurement [73]. For numerical simulations, the model size L usually meets l min ( L ( l max [131]. A fractal spatial distribution of fracture barycentres can be modelled through a multiplicative cascade process governed by a prescribed D value, while fracture lengths can be sampled from a power law distribution having an exponent a [131]. Fracture orientations can be assigned isotropically or based on a statistical distribution. Fractal fracture networks can then be generated through a synthesis of the different geometrical attributes modelled by independent random variables (Fig. 3). A D value of 2 (in 2D) and 3 (in 3D) represents a homogeneous spatial distribution, i.e. ''space filling". As D decreases, the fracture pattern becomes more clustered associated with more empty areas. A small a value corresponds to a system dominated by large fractures, while a ? 1 relates to a pattern with all fractures having an equal size (i.e. l min ). The D and a values as well as their relationship may control the connectivity, permeability and strength of fractured rocks [131][132][133][134]. More interestingly, when a = D + 1, the fracture network is self-similar and the connectivity properties are scale invariant [131]. A self-similar fracture pattern statistically exhibits a hierarchical characteristic whereby a large fracture inhibits the propagation of smaller ones in its vicinity, but not the converse [32,130]. Implementation of such a hierarchical rule together with subcritical fracture growth laws leads to a new DFN model that can simulate the sequential stages of fracture network formation associated with nucleation, propagation and arrest processes (Fig. 4a) [135]. The conventional Poisson DFN models: (a) a 2D random fracture pattern conditioned by field data from the Sellafield site, Cumbria, UK [95] and (b) a conceptual model (no scale specified) of 3D random fracture network with three orthogonal sets of disc-shaped fractures [64].
The assumption of a linear/planar fracture shape in the Poisson model and fractal DFN model may be simplistic. Many field observations have shown that natural fracture geometries can be curved and irregular [9]. To simulate the non-planarity of natural fractures, Lei et al. [40] developed a novel DFN model that accommodates discrete-time random walks in recursive self-referencing lattices to extrapolate fracture geometries into larger domains based on a small sample (Fig. 4b). This model was developed mainly for generating cross-cutting fracture sets. For more complex networks involving abutting relations, the random walk tech-nique proposed by Horgan and Young [136] for simulating polygonal crack patterns in soil may be employed. Furthermore, the invasion percolation method that has been used to model channel networks [137] may provide a way to simulate some highly branched and tortuous fracture systems.
The stochastic DFN method, in essence, treats problems in a probabilistic framework and regards the real physical system as one possibility among simulated realisations sharing the same statistics. Hence, a sufficient number of realisations based on a Monte Carlo process are necessary to predict a bounded range. In Fig. 4. Some new stochastic DFN models: (a) trace map views of a 3D sequential DFN model that simulates the nucleation, growth and arrest processes of natural fractures (the ratio of the domain size L to the characteristic nuclei length l N is 100) [135] and (b) a 2D self-referencing DFN model that extrapolates a larger scale pattern recursively from a smaller one through random walks with the preservation of natural fracture curvatures [40]. practice, a balance exists between the benefits of collecting detailed information to create more representative DFNs and the increased cost of field measurements. The uncertainty can be reduced by constraining the random process with deterministic data, e.g. forcing the 3D DFN generator to reproduce a 2D trace pattern such as one exposed on tunnel walls [138]. Calibration and validation of stochastic DFN models are important for solving real problems and can be conducted based on the in-situ data from field mapping and/or hydraulic tests [138][139][140][141][142][143][144]. The random nature of the stochastic DFN method may be regarded as an advantageous aspect, because uncertainty is unavoidable when analysing complex geological systems, for which single-valued predictions from deterministic methods may be more risky [7]. However, it is still very important to continue improving the realism and accuracy of stochastic DFN models, because the predicted results (e.g. permeability, strength, fragmentation, etc.) featured by upper and lower bounds based on unrealistic DFNs may be systematically biased from the truth. Understanding of the variability and uncertainty in data collection and DFN generation are thus critical for engineering applications [145]. Developments are needed towards (i) a more thorough characterisation of the underlying statistics, such as multifractals for which a single scaling exponent is not sufficient [110], and (ii) a more precise and efficient generator to create DFNs respecting more details of real fracture systems, such as the diversity of individual fracture shapes and morphology [9], the topological complexity in fracture populations [114,146] and the correlation between geometrical properties [129,147,148]. The important difference between 2D and 3D fracture networks with respect to connectivity and permeability [64,90,105] renders another advantage of the stochastic method, which has an intrinsic capability of generating 3D networks.

Geomechanically-grown fracture networks
Extensive studies have been conducted to interpret the geological history and the formation mechanism behind field observations of natural fracture systems (e.g. patterns, statistics and minerals) [8,9,[28][29][30]149]. The increased knowledge of fracture mechanics [5] promoted the development of geomechanicallybased DFN models that incorporate the physics of fracture growth and simulate fracture network evolution as a geometrical response to stress and deformation. By applying a geologically-inferred palaeo-stress/strain condition, natural fracture patterns may be reproduced by such a DFN simulator that progressively solves the perturbation of stress fields and captures the nucleation, propagation and coalescence of discrete fractures. Different numerical methods have been proposed and the one based on linear elastic fracture mechanics (LEFM) is frequently adopted.
In the LEFM model, fracture patterns can be simulated by four main steps in an iterative fashion [150,151]: (1) generation of initial flaws to mimic the process that natural fractures initiate from microcracks, (2) calculation of the perturbed stress field in the rock caused by the presence and evolution of fractures under an imposed boundary condition, (3) derivation of the stress intensity factor (K I ) at the tip of each fracture, and (4) propagation of fractures which satisfy a growth criterion, e.g. a subcritical law K O 6 K I 6 K IC , where K O is the stress corrosion limit and K IC is the material toughness [152]. The stress field and stress intensity factor can be calculated based on analytical solutions [150] or (most commonly) numerical methods such as the boundary element method (BEM) [153] and finite element method (FEM) [151,154]. The propagation length in each growth iteration can be derived according to a power law relation with the energy release rate G (related to K I ) associated with a velocity exponent j (or subcritical index n) [152], while the propagation angle may be computed if the curvature and coalescence effects are important, especially when the tectonic stress field is quite isotropic [29].
The development of fracture networks is a sophisticated feedback-loop process, in which the complexity of growth dynamics is directly related to the complexity of the developing structures. Specifically, the propagation of a fracture is influenced by the mechanical interaction with others, and the propagated fracture geometries can conversely generate stress perturbations into the system. The mechanical interaction of fractures was found strongly dependent on the velocity exponent j: an increased j tends to promote a localised fracture pattern [150,153,[155][156][157]. The fracture pattern evolution is also affected by the density [150,157] and orientation [158] of the initial flaws, and the 3D layer confinement effect [155,156]. Such a fracture mechanics model has been developed to mimic the evolution of a 2D single set of straight fractures (Fig. 5a) [150,153,155], 2D orthogonal sets of straight fractures [157], 2D curved fracture patterns (Fig. 5b) [29,151,156,158,159], and 3D curved fracture geometries (Fig. 5c) [160]. The generated DFN pattern can be further used to evaluate the connectivity [91,157], permeability [161] and solute transport properties [162] of natural fracture systems.
Apart from the LEFM approach, fracture patterns can also be simulated using other numerical methods. Cowie et al. [163,164] developed a lattice-based rupture model to simulate the antiplane shear deformation of a tectonic plate and the spatiotemporal evolution of a multifractal fault system. Tang et al. [165] used a damage mechanics FEM model to simulate the evolution of parallel, laddering or polygonal fracture patterns formed under different boundary conditions, i.e. uniaxial, anisotropic and isotropic tectonic stretch, respectively. Spence and Finch [166] employed the discrete element method (DEM) to simulate the fracture pattern development in a sedimentary sequence embedded with stratified nodular chert rhythmites. Asahina et al. [167] coupled the finite volume multiphase flow simulator (i.e. TOUGH2) and a latticebased elasticity and fracture model (i.e. Rigid-Body-Spring Network) to simulate the desiccation cracking in a mining waste material under a hydromechanically coupled process. The finite-discrete element method (FEMDEM) has also been used to simulate the formation and evolution of geological fracture patterns [168,169].
The geomechanically-grown DFN model, as a process-oriented approach, has the advantage of linking the geometry and topology of fracture networks with the conditions and physics of their formation. Another merit is the automatic correlation between the geometrical attributes of individual fractures (e.g. length, orientation, aperture and shear displacement) linked by the governing physics. To solve practical problems, such a DFN generator can be constrained by the measurement of rock properties (e.g. the subcritical index measured from core samples) and the information of geological conditions (e.g. stress, strain, pore pressure and diagenesis) to achieve rational predictions [170]. However, difficulty and uncertainty still exist in creating fracture patterns consistent with the real systems for which coupled tectonic, hydrological, thermal and chemical processes may be involved.

Summary with a view towards HM modelling
The three types of DFN models have each its own strength and also suffer from some limitations, as listed in Table 1. The geologically-mapped DFNs can preserve many realistic features of fractures but it is hard to apply this method to characterise deep rocks and 3D structures. The stochastic DFN approach has the merits of simplicity and efficiency as well as applicability for 3D problems. However, the strong geological hypotheses of fracture geometries and topologies used in its construction may ignore some important underlying mechanical and tectonic constraints and thus possibly lead to errors or large uncertainties. The geome-chanical DFN models, which may capture some mechanical characteristics of natural fractures, are sensitive to the assumed/measured rock properties and inferred palaeostress fields. If they are to be improved to become the preferred useable model, coupling with hydrological, thermal and chemical mechanisms may be needed to reproduce actual geological systems. To generate better fracture networks, a future research direction that is attracting much effort is the development of hybrid DFN models that assimilate the advantages of different approaches. Some of the current DFN models have already exhibited some of such features. For example, Kattenhorn and Pollard [59] used mechanical simulation to correct the 3D fault structures interpreted from seismic survey data. In the sequential stochastic DFN model proposed by Davy et al. [135], fractures develop following the subcritical growth law and their interactions are governed by an arrest mechanism. In the self-referencing DFN model by Lei et al. [40], a geologically-mapped fracture pattern is populated to larger domains based on a random walk algorithm that preserves key geometric and geomechanical attributes.
The created DFN geometries can serve as important inputs for modelling the geomechanical and/or hydrological behaviour of fractured rocks. To reduce the numerical difficulties induced by   DFN complexities whilst capturing the key characteristics, simplifications may be applied. In conventional hydrological studies, inactive fractures (i.e. dead-end or isolated fractures) and very small fractures are often removed from original DFNs based on an assumption of their negligible contribution to the overall connectivity and permeability [67,87,95,99]. Such treatments may be applicable for low permeability rocks with well-connected fracture topologies, inside which the matrix flow can be neglected. However, when the fractures are poorly connected and/or the rock matrix is highly permeable, the effects of fracture-matrix fluid transfer and flow in the matrix can significantly affect (or even govern) the hydrological processes in fractured rocks [171] and thus the features of inactive fractures may need to be preserved. Attentions may also be drawn when setting a truncation threshold for deleting small-sized fractures which were commonly believed to have little contribution to the flow. However, several studies have suggested that when a > D + 1, the connectivity of fracture networks is ruled by fractures much smaller than the system size [131,172] and thus the removal of very small fractures may lead to important biases. Apart from these potential hydrological effects, even more importantly, the presence of inactive fractures may also influence the mechanical behaviour of fractured rocks (e.g. elastic or brittle deformation). New fractures may propagate from the dead-end tips of pre-existing fractures under stress concentrations and link with other fractures to form critical pathways for fluid flow. More discussions about choosing appropriate DFNs for geomechanical and HM modelling will be presented in the following sections.

Modelling of the geomechanical behaviour of fracture networks
The numerical methods for geomechanical modelling of fractured rocks can be categorised as continuum and discontinuum approaches with the classification based on their treatment of displacement compatibility [12]. The preference for a continuum or discontinuum modelling scheme depends on the scale of the problem and the complexity of the fracture system [22]. The continuum approach has the advantage of its greater efficiency to handle large-scale problems, whereas the discontinuum method can more straightforwardly incorporate complex fracture networks and model multi-fracturing and fragmentation processes. In this section, commonly used models for simulating geomechanical behaviour of fractured rocks will be reviewed: continuum/extendedcontinuum, block-type discontinuum, particle-based discontinuum and hybrid finite-discrete element approaches. It is worth mentioning that the classification here is not intended to be absolute, since the boundary between the continuum and discontinuum methods has become very vague. Some advanced continuum techniques have included contact algorithms and fracture mechanics to consider discontinuities, while many discontinuum models are able to deal with continuous deformations.

Continuum/extended-continuum models
The conventional continuum approach treats a rock domain as a continuous body that can be solved by the finite element method (FEM) or finite difference method (FDM). It may be applicable to a fractured rock with only a few or a large number of fractures [12]. If the system consists of only a few discontinuities associated with only a small amount of displacement/rotation, the discrete fractures can be modelled by special ''interface elements" (or ''joint elements") that are forced to have fixed connectivity with the solid elements [173]. However, such a treatment renders difficult to handle the dynamics and large displacement problems of natural fracture systems. When the density of DFN fractures is very high, the modelling domain may be divided into a finite number of grid blocks assigned with equivalent properties derived from homogenisation techniques (Fig. 6a). The equivalent properties, such as bulk modulus and strength parameters, are usually calculated using empirical formulations that consider the degradation effect caused by pre-existing fractures [174,175] or analytical solutions based on the crack tensor theory [176,177]. The crack tensor theory calculates the volume averaged parameters of all fractures with respect to their geometrical properties (e.g. length, orientation and aperture) and was extended to consider the coupling between stress and fluid flow [46,178]. Such a crack tensor method has been integrated into the FEM [179,180] and FDM [101,181] solvers to model the geomechanical and HM behaviour of fractured rocks. The simulation results may be sensitive to the grid block discretisation especially when the block size is significantly smaller than the representative elementary volume (REV) [101]. The homogenisation-based continuum model may not adequately consider the connectivity effect of very long fractures that penetrate numerous grid blocks and the results may be even worse if apertures are very heterogeneous and positively correlated with fracture sizes. Thus, it may not be applicable to a fractal fracture system with high variability in its density distribution (i.e. a small fractal dimension D) and/or a large proportion of long fractures having a size comparable to the problem domain (i.e. a small power law length exponent a). The two conventional continuum schemes (i.e. the one employing interface elements and the one using homogenised properties) may be combined to explicitly model large discontinuities (e.g. dominant faults) using interface elements and then to characterise each isolated block as continuum bodies with their bulk properties parameterised according to the distribution of small fractures. Furthermore, the crack tensor method cannot consider the interaction between fractures and blocks as well as the resulting localised deformation and damage in the rock.
To more explicitly capture the effects of discrete fractures, an extended-continuum model has been developed by assuming fractures to have certain numerical width (for connectivity preservation) and representing them as arrays of grid elements with softening and weakening properties based on a very fine finite difference mesh (Fig. 6b) [41,182]. It treats the fractured rock as a composite elasto-plastic solid system, in which the failure of intact rocks or stress-displacement behaviour of fractures can be modelled by a Mohr-Coulomb criterion with tension cut-off. Similar ''weak material" representation of fractures has also been implemented in the rock failure process analysis (RFPA) code (a damage mechanics FEM model) [183] and the cellular automation model [184]. Such a composite continuum model with explicit DFN representations may be more suitable for simulating weakly cemented fractures (i.e. mineral filled veins), whereas the physical rationale is not intuitive if being applied to unfilled discontinuities having clean wall surfaces. Furthermore, another important extendedcontinuum model is the extended finite element method (XFEM), which can represent fracture propagations without re-meshing and has been recently developed for simulating fluid-driven fracturing [185,186].

Block-type discontinuum models
The block-type discontinuum models include the distinct element method (DEM) with an explicit solution scheme and the discontinuous deformation analysis (DDA) method with an implicit solution form. In this discontinuum modelling framework, the fractured rock is represented as an assemblage of blocks (i.e. discrete elements) bounded by a number of intersecting discontinuities. The geometry of the interlocking block structures can be identified first by e.g. employing the techniques of combinatorial topology [187]. In the subsequent mechanical computations, these blocks can be treated as rigid bodies or deformable subdomains (further discretised by finite difference/volume grids) with their interactions continually tracked by spatial detections during their deformation and motion processes.

Distinct element method (DEM)
The DEM method was originated by Cundall [188,189] and gradually evolved to the commercial codes UDEC and 3DEC for solving 2D and 3D problems, respectively [190,191]. Its basic computational procedure can be summarised as four steps [192]: (1) the contacts between blocks are identified/updated through a space detection, (2) the contact forces between discrete bodies are computed based on their relative positions, (3) the acceleration induced by force imbalance for each discrete element is calculated using Newton's second law, and (4) the velocity and displacement are further derived by time integration with new positions determined. An explicit time marching scheme is applied to solve the problem iteratively until the block interaction process to be simulated has been completed. The mechanical interaction between blocks is captured by a compliant contact model that accommodates virtual ''interpenetrations" governed by assumed finite stiffnesses to derive normal and tangential contact forces. The empirical joint constitutive laws derived from laboratory experiments [14,193] can be implemented into the interaction calculation in an incremental form to simulate joint normal and shearing behaviour [194][195][196][197]. A viscous damping parameter may be introduced to reduce dynamic effects for modelling quasi-static conditions [198].
The DEM approach is able to capture the stress-strain characteristics of intact rocks, the opening/shearing of pre-existing fractures and the interaction between multiple blocks and fractures. Combined with DFN models, it has been widely applied to study the mechanical behaviour of fractured rocks. Zhang and Sanderson [44] studied the critical deformation behaviour of fractured rocks and observed an abrupt increase in the deformability associated with a power law scaling when the fracture density exceeds the mechanical percolation threshold (slightly higher than the geometrical threshold) (Fig. 7). Min and Jing [199] examined the scale dependency of the equivalent elastic properties of a fractured rock based on multiple DFN realisations conditioned by the same fracture statistics (Fig. 8). In their study, a technique to derive the fourth-order elastic compliance tensor has also been developed for equivalent continuum representations. Min and Jing [200] further found that the equivalent mechanical properties (i.e. elastic modulus and Poisson's ratio) of fractured rocks may also be stress dependent. With the increase of stress magnitudes, the equivalent elastic modulus significantly increases, while the Poisson's ratio generally decreases but can be well above 0.5 (i.e. the upper limit for isotropic materials). Noorian-Bidgoli et al. [201] extended this DEM-DFN model to a more systematic framework to derive the strength and deformability of fractured rocks under different loading conditions, which is further applied to study the anisotropy [202] and randomness [203] of the strength/deformability of stochastic DFNs. Recently, Le Goc et al. [204] integrated 3D DFNs into the 3DEC simulator and investigated the effects of fracture density, sizes and orientations on the magnitude and scaling of the equivalent elastic modulus of fractured rocks. Havaej et al. [205] conducted 3D DEM-DFN modelling to analyse the effects of geological structures on the quarry slope failure mechanism. In addition, a considerable number of similar DEM-DFN models have been developed where the main motive is to study the effect of stresses on fluid flow. Such models were applied to explicitly capture the fracture opening, closing, shearing and dilational characteristics in complex fracture networks under in-situ stresses, after which the fluid flow implications were investigated [43-45, 96,98,206-209] (more details are presented in Section 4). Furthermore, some models that incorporate the pore fluid pressure effect show that the pore pressure level can also significantly influence the mechanical and hydrological properties of fractured rocks [47,210].
The classic DEM formulation cannot simulate the propagation of new fractures in intact rocks driven by stress concentrations, although plastic yielding can capture some aspects of the rock mass failure process [211]. Such a shortcoming was recently addressed by introducing a Voronoi or Trigon discretisation in matrix blocks that allows fracturing along the internal ''grain" boundaries governed by tensile and shear failure criteria [212][213][214][215][216]. The DFN representation can also be integrated into such a DEM model using a dual-scale tessellation, in which the primary grid represents natural fractures and the secondary discretisation mimics microstructures [217]. The uncertainty of the Voronoi/Trigon DEM model due to mesh dependency indicates the necessity of realistic representation of natural grain shapes and textures [218].

Discontinuous deformation analysis (DDA)
The DDA method was proposed by Shi and Goodman [219,220] to compute the deformation and motion of a multi-block system. The discretisation of a DDA system is quite similar to the one for the DEM, i.e. the medium is dissected into blocks by intersecting discontinuities. However, a fundamental difference between the two methods lies in their computational frameworks. The DEM treats the kinematics of each block separately based on an explicit time-marching scheme, while the DDA calculates the displacement field based on a minimisation of the total potential energy of the whole blocky system and an implicit solution to the established system of equations through a matrix inversion. Thus, the DDA method has an important advantage of fast convergence with unconditional numerical stability compared to the DEM method that requires a time step smaller than a critical threshold [12]. Important extensions of the original DDA method include the finite element discretisation of rock matrix [221], the sub-block technique (similar to the Voronoi DEM) for simulating fracturing processes [222], the formulation for modelling coupled solid deformation and fluid flow [223,224], and the development of 3D models [225]. The mechanical behaviour of fractured rocks has also been investigated based on DDA-DFN simulations, with emphasis on studying the stability of underground excavations and slope engineering ( Fig. 9) [226][227][228]. Recently, Tang et al. [229] combined the RFPA (a damage mechanics FEM model) with the DDA method to capture crack propagations and block kinematics of a rock slope system.

Particle-based discontinuum models
The particle-based discontinuum model was originally introduced by Cundall and Strack [230] to simulate granular materials such as soils/sands and gradually evolved to a commercial code, i.e. particle flow code (PFC) [231]. Similar to the block-type DEM method, PFC calculates the inertial forces, velocities and displacements of interacting particles by solving Newton's second law through an explicit time-marching scheme. The discrete particles are assumed to be rigid with a circular (in 2D) or spherical (in 3D) shape, and can have variable sizes which are usually much larger than the physical grain scale. To extend PFC to model rock materials, Potyondy and Cundall [232] developed a bondedparticle model (BPM), in which intact rocks are represented as assemblages of cemented rigid particles and the macroscopic fracturing is simulated as a result of the breakage and coalescence of the inter-particle cohesive bonds. The interaction between particles can be characterised by two types of bond models in PFC, i.e. the contact bond model and the parallel bond model. A contact bond serves as a linear elastic spring with normal and shear strengths, and transmits forces via the contact point between two particles. A parallel bond with certain normal and shear stiffnesses joins two particles to resist against separation under tension, shear and rotation. The parallel bond model is more suitable for simulating rock materials as it can capture the cement-like behaviour with resistance to bending moments. To overcome the original deficiency of PFC in reproducing a realistic rock strength ratio (i.e. the ratio of uniaxial compressive strength to tensile strength) and macroscopic friction angle, a hierarchical bonding structure can be built based on a cluster logic [232] or a clump logic [233]. The cluster/clump approach mimics the interlocking effect of irregular grains by defining a higher value of intra-cluster bond strength (i.e. the strength between particles in the same cluster) than that of the strength between cluster boundaries.
The BPM representation using particles with idealised circular/spherical shapes can introduce unphysical asperities on discontinuity surfaces, resulting in an additional resistance to frictional sliding. To suppress such an artificial roughness effect, a smoothjoint contact model (SJM) was proposed to simulate fracture wall behaviour based on the actual geometry and morphology of discontinuities and independent of the arrangement of local particles in contact [234]. The smooth contact is assigned to all particle pairs lying on the fracture interface but belonging to opposite matrix blocks, so that they can overlap and pass through each other (Fig. 10a). The contact forces are calculated based on the relative displacements and the smooth-joint stiffness in the normal and tangential directions of the local surfaces. The SJM was found to be able to capture the shear strength and dilational behaviour of natural fractures associated with significant scale effects [235,236]. By simulating intact rock matrix using the BPM and discrete fractures using the SJM, a synthetic rock mass (SRM) modelling approach (Fig. 10b) has been developed to characterise the mechanical properties of fractured rocks including peak strength, damage, fragmentation, brittleness, anisotropy and scale effects [234]. Compared to the Hoek-Brown empirical approach for presumed isotropic rock masses, the SRM method that integrates explicit DFN representations has an advantage to derive the orientation-specific strength of naturally fractured rocks and consider the influence of fracture length distribution and connectivity [237]. The SRM model has been applied to reproduce the failure behaviour of veined core samples under uniaxial compression tests [238], to estimate the mechanical REV of a jointed rock mass near an underground facility [239], to simulate fracture propagation in jointed rock pillars [240,241], and to evaluate the stability of wedges around a vertical excavation in a hard rock [242].
An open source particle-based modelling platform named YADE [243,244] has recently been developed as an alternative to the commercial code PFC. YADE represents intact rocks using glued discs/spheres and models the fracturing process based on the rupture of inter-particle bonds with the contact bond algorithms following a similar logic to PFC. To reproduce the high ratio of the compressive strength to the tensile strength and the non-linear failure envelope of brittle rocks, the concept of ''interaction range" was introduced by Scholtès and Donzé [245]. They mimic the microstructural complexity by assembling constitutive particles in neighbouring zones (not only the particles in direct contact). A joint contact logic equivalent to the SJM has also been imple-mented in YADE to avoid the particle interlocking effect between sliding fracture surfaces [246]. By integrating 3D fractal DFNs into the YADE BPM model associated with the smooth joint contact treatment (Fig. 11), Harthong et al. [134] studied the influence of fracture network properties (i.e. fractal dimension D, power law length exponent a and fracture intensity P 32 ) on the mechanical behaviour of fractured rocks. The strength and elastic modulus of rock masses decrease if P 32 increases (i.e. more fractures) or a decreases (i.e. higher proportion of larger fractures), while the spatial heterogeneity and scaling of the mechanical properties are   11. Integration of (a) a fractal DFN into (b) the YADE bonded-particle model (BPM) for mechanical modelling of fractured rocks [134]. affected by D. Such a combined BPM-DFN model has also been applied to analyse the stability of fractured rock slopes, which is controlled by the strengths of both pre-existing fractures and intact rocks [247,248].

Hybrid finite-discrete element models
The hybrid finite-discrete element method (FEMDEM or FDEM) combines the finite element analysis of stress/deformation evolution with the discrete element solutions of transient dynamics, contact detection and interaction. In such a discontinuum modelling scheme, the internal stress field of each discrete matrix block is calculated by the FEM solver, while the translation, rotation and interaction of multiple rock blocks are traced by the DEM algorithms. Preexisting fractures in rocks are treated as the internal boundaries of rock volumes. The FEMDEM approach also provides a natural solution route to modelling the transitional behaviour of brittle/quasibrittle materials from continuum to discontinuum (i.e. fracturing processes) by integrating fracture mechanics principles into the formulation. This section will review the two most commonly used FEMDEM models, i.e. the commercial software ELFEN [249] and an open source platform Y-code [250], which have been broadly used to simulate the mechanical processes in geological media containing pre-existing discontinuities. Although the two FEMDEM models originated from the same concept [251], a fundamental difference exists in their discontinuum implementations: ELFEN is based on a partially-discontinuous discretisation with joint elements inserted only in pre-existing or propagated fractures, while Y-code is based on a fully-discontinuous mesh where joint elements are embedded between the interfaces of all finite elements (i.e. along fractures and inside intact domains). Some other differences also lie in the computation of intact material deformation and determination of local failure. For more details, the reader is referred to the work by Lisjak and Grasselli [25].

ELFEN
The ELFEN code models the degradation of an initial continuous domain into discrete bodies by inserting cracks into a finite ele-ment mesh. A nodal fracture scheme was introduced by constructing a non-local failure map for the whole system [252]. The feasibility of local failure is determined based on the evolution of nodal damage indicators. The fracturing direction (if failure occurs) is calculated based on the weighted average of the maximum failure strain directions of all surrounding elements. A new discrete fracture is then inserted along the failure plane with the local mesh topology updated through either the ''intra-element" or ''interelement" insertion algorithm with adaptive mesh refinement applied if necessary [253]. ELFEN provides various material constitutive models including the elastic, elasto-plastic and visco-plastic laws, and many brittle/quasi-brittle failure models including the rotating crack model, the Rankine material model, and the compressive fracture model (i.e. Mohr-Coulomb failure criterion coupled with a tensile crack model) [253,254]. The capability of ELFEN for modelling 3D fracturing processes has also been demonstrated through the simulation of laboratory strength tests [255]. The explicit DFN fracture geometries generated from e.g. FracMan can be imported into the ELFEN platform by embedding fracture entities into rock solids and representing each fracture as opposed free surfaces [256,257]. To mesh such complex systems, special geometrical treatments may be involved to avoid ill-posed elements caused by subparallel fractures intersecting at a very small acute angle or a fracture tip terminating at the vicinity of another fracture [249]. Both pre-existing and newly propagated fractures are assigned with contact properties, e.g. fracture stiffness and friction coefficient, to simulate solid interactions through discontinuity surfaces [257]. The degradation of natural fractures during shearing can also be modelled by introducing roughness profiles [258].
The combined FEMDEM-DFN model has been applied to tackle the geomechanical problems for various engineering applications [259]. The presence of natural fractures may dominate the strength of slender pillars but have a reduced influence for wider pillars (Fig. 12a) [256]. The orientation and length distribution of DFN fractures also affect the failure mode of the pillar structures, which can exhibit splitting with lateral kinematic releases or shearing of critically inclined pre-existing fractures linked by new cracks Fig. 12. Integration of DFN geometries into the FEMDEM model of ELFEN for modelling strength of (a) a prefractured pillar [256] (note: r ci is the uniaxial compressive strength of intact rocks, P 21 denotes the fracture intensity, i.e. total length of fractures per unit area) and (b) an open pit slope [265]. through intact rock bridges [260]. This synthetic numerical model has also been used to investigate the progressive failure of rock slopes (Fig. 12b) [261,262], which in reality is usually triggered by both the reactivation of natural fractures and the propagation of new cracks [263]. The FEMDEM model is well suited to mimic the staged failure processes of rock slopes including initiation, transportation/comminution and deposition, which involve yielding and fracturing of intact materials, shearing of fracture surfaces and translational/rotational instabilities [264]. The rock mass fabrics and rock bridge properties can have important influences on the stability of large-scale open pit slopes [261]. The cavinginduced rock mass deformations and associated surface subsidence may be controlled by the orientation of joint sets and the location/ inclination of major faults [265]. The FEMDEM-DFN model has also been applied to capture the progressive failure of the roof strata of a coal mine tunnel [266]. All these engineering applications highlighted the advantage of the FEMDEM-DFN technique with the capability of simulating the reactivation/interaction of preexisting fractures and initiation/propagation of new cracks.

Y-code
During the 1990s, many algorithmic solutions for 2D and 3D FEMDEM simulation were developed by Munjiza et al. [251,267] and Munjiza and Andrews [268,269]. Extensive developments and applications of the FEMDEM model have been conducted after the release of the open source Y-code [250], and different versions have emerged including the code developed collaboratively by Queen Mary University and Los Alamos National Laboratory [270][271][272], the Y-Geo program by Toronto University [25,273,274], and the VGeST (recently renamed ''Solidity") platform by Imperial College London [15,39,[275][276][277][278][279][280][281]. The Y-code FEMDEM model accommodates the finite strain elasticity coupled with a smeared crack model and is able to capture the complex behaviour of discontinuous rock masses involving deformation, rotation, interaction, fracturing and fragmentation.
In the Y-code, the fractured rock is represented by a discontinuous discretisation of the model domain using three-noded triangular (in 2D) or four-noded tetrahedral (in 3D) finite elements and four-noded (in 2D) or six-noded (in 3D) joint elements embedded at the interface between finite elements. An important difference with the ELFEN code is that the joint elements are inserted for all edges (in 2D) or surfaces (in 3D) of finite elements. The deformation of the bulk material is captured by the linear-elastic constant-strain finite elements with the impenetrability enforced by a penalty function and the continuity constrained by a constitutive relation [267], while the interaction of matrix bodies through discontinuity interfaces is simulated by the penetration calculation [269]. The joint elements are created and embedded between triangular/tetrahedral element pairs before the numerical simulation, and no further remeshing process is performed during later computations. Pre-existing fractures can be represented by a series of joint elements which are initially overlapped (but opposite sides are separately defined) free surfaces [39]. The brittle failure of intact materials is governed by both fracture energy parameters (for mode I and mode II failure) and strength properties (e.g. tensile strength, internal friction angle and cohesion) [25]. A numerical calibration can be conducted to achieve consistency between the input material strength parameters and simulated macroscopic response [282]. Code development for modelling 3D crack propagation has also been achieved by different research groups [272,[283][284][285][286].
The advantage of the FEMDEM model for simulating the degradation of continuum into discrete pieces promoted the application to tackle various engineering problems, such as rock blasting [287], fracture development around excavations in isotropic/anisotropic intact rocks [288][289][290][291] and mountain slope failure [292]. It has also been used to model the mechanical behaviour of fractured rocks embedded with pre-existing fractures. Lisjak et al. [293] integrated DFN crack arrays into the FEMDEM model to imitate the anisotropy of an argillaceous rock. The FEMDEM technique has been applied to model geologically-mapped DFNs under various stress conditions [15,39,40,280,294]. It can capture realistic geomechanical phenomena such as deformation and rotation of matrix blocks, opening, shearing, and dilation of pre-existing fractures as well as the propagation of new cracks [281]. Due to the presence of natural fractures, strongly heterogeneous stresses are distributed in the rock mass (Fig. 13a) [39]. The in-situ stress can introduce new cracks into the natural fracture system which potentially increases the network connectivity (Fig. 13b) [15,280]. Highly variable apertures can also be generated (Fig. 13c), which is related to the complexity of DFN geometries and the condition of far-field stresses [39,278,280]. The stress-dependent mechanical response of fractured rocks (e.g. deformation and rotation of matrix blocks, opening, shearing and dilation of pre-existing fractures and new crack propagation) can significantly affect the fluid flow in fracture networks, which will be discussed in more details in Section 4.

Summary with a view towards HM modelling
As can be seen from above, geomechanical modelling of fractured rocks can be achieved by continuum or discontinuum approaches, which have important differences in the conceptualisation of geological media and the treatment of displacement compatibility [12]. Table 2 presents a comparison of the continuum and discontinuum models. The continuum modelling scheme mainly reflects the material deformation of a geological system from a more overarching view and attempts to bypass the geometrical complexity by using specific constitutive laws and equivalent material properties derived from homogenisation techniques. However, it cannot adequately consider the effects of stress variations, fracture interactions, block displacements and rotations. More importantly, the applicability of a homogenisation process is based on the assumption of an REV, which may not exist for natural fracture systems [73]. On the other hand, the discontinuum scheme treats the system as an assemblage of interacting individual components and permits the integration of complex constitutive laws for rock materials and fracture interfaces. A discontinuum model can be established at a specific scale of investigation without presuming the existence of an REV. However, some of the input parameters (e.g. bonding strength, joint stiffness) may need to be determined by indirect numerical calibrations rather than from physical measurements. Furthermore, the computational time for solving discontinuous problems can be considerably larger than that for continuum models. To take advantage of the two modelling technologies for tackling practical issues, a discontinuum model may be used to derive the REV size (if it exists) as the onset to treat a geological system as a continuum.
In the conventional block-type DEM and DDA modelling, the rock mass is dissected into blocks through topological identifications, during which the dead-ends and isolated fractures are usually removed for simplicity [44,199,206,221,227]. Such treatments may be acceptable for modelling extremely hard rocks or highly fragmented rocks, for which the effect of new crack propagations may be negligible. However, for most rock types, fractures are likely to propagate from the tips of pre-existing fractures driven by concentrated local stresses. This fracture propagation and coalescence process can create new pathways for fluid migration and thus affect the hydrological properties (e.g. permeability) of fractured rocks, as has been revealed in numerical simulations [15,39]. Furthermore, the curved natural fractures, which are often assumed to be straight line (in 2D) or planar disc (in 3D) in the Poisson DFN model, may accommodate localised dilations under high differential stresses with localised increases in aperture and this may have important HM effects [15,47,280]. Hence, towards a better HM modelling of fractured rocks, the capability of precisely modelling the geomechanical response is prerequisite, for which the important characteristics of fracture dead-ends and curvatures need to be appropriately treated.

Modelling of the hydrological behaviour of fracture networks under geomechanical effects
The presence of fractures can generate stress perturbations in the rock, such as rotation of stress fields, stress shadows around discontinuities and stress concentration at fracture tips [5]. The resulting heterogeneous stress distribution may lead to variable local normal/shear stresses loaded on different fractures having distinct sizes and orientations, and produce various fracture responses such as opening, closing, sliding, dilatancy and propagation. Since the conductivity of fractures is critically dependent on the third power of fracture apertures [295], the geomechanical conditions can considerably affect the hydrological properties of fractured rocks including fluid pathways, bulk permeability and mass transport [207]. Numerical models that integrate explicit DFNs and non-linear rock/fracture behaviours provide powerful (and so far irreplaceable) tools to investigate the geomechanical effects on fluid flow in complex fracture networks [296].

Fluid pathways
Fracture networks usually serve as the major pathways for fluid transport in subsurface rocks, especially if the matrix is almost impermeable compared to the fractures [297]. The partitioning of fluid flow within a fracture population relies on the spatial connectivity of fracture geometries and the transmissivity of individual fractures, both of which can be affected by the geomechanical conditions.
Zhang et al. [45] used the UDEC DEM code to study the deformation of a fractured rock based on a geologically-mapped DFN pattern and found the closure of fractures under applied in-situ stresses can re-organise the fluid pathways. They also found the Variation of fracture apertures in a geologically-mapped DFN network that accommodates further new crack propagations in response to a biaxial stress condition [15]. (c) Highly inhomogeneous aperture distribution within a single fracture of a 3D persistent DFN system under a polyaxial stress condition [278]. closed joints of one set can hinder fluids from passing through not only themselves but open fractures of another set due to the net effect. Zhang and Sanderson [43] applied the same technique to different types of outcrop DFN patterns with/without systematic fracture sets. The fluid flow in systematic fracture networks tend to be dominated by the primary joint set if it is relatively open, while for non-systematic networks containing fairly randomly oriented small fractures, the flow channels tend to align the direction of the maximum principal stress. Min et al. [206] incorporated the fracture shear dilation behaviour in the DEM modelling of a stochastic DFN having a power law distribution of fracture lengths and a uniform distribution of initial (i.e. zero stress) apertures. They observed that, under high differential stresses, a small portion of fractures which have critical/near-critical orientations, good connectivity and long lengths would dilate and form large flow channels (Fig. 14a). The ''critical orientations" here correspond to a range of fracture orientations that would allow discontinuities with no cohesive strength to slide under the given in-situ differential stress condition. The localised features would be augmented if the initial apertures are broadly distributed (e.g. following a lognormal distribution) and correlated with fracture lengths [98]. Latham et al. [15] employed the FEMDEM method combined with a smeared crack model to study the geomechanical response and fluid flow in an outcrop-based DFN. They found that bent natural fractures under high differential stresses may exhibit evident dilational jogs and can be linked by newly propagated cracks to form major fluid pathways. Similar localised flow channels created by the connection of pre-existing fractures have also been observed by Figueiredo et al. [41] using a FDM simulator. Sanderson and Zhang [47,96] calculated the fluid flow in the third dimension of sedimentary rocks using an analytical pipe formula based on the deformed 2D fracture networks. They found vertical flow becomes extremely localised when the pore fluid pressure exceeds a critical level and very large aperture channels emerge (especially in some preferentially-oriented and well-connected curved fractures) (Fig. 14b). The flow distribution also exhibits significant multifractality when the loading condition approaches the critical state. Lei et al. [278] studied the fluid flow in a 3D persistent fracture network under various polyaxial stress conditions. They observed a transition from homogeneous flow under an isotropic stress field to highly localised patterns under a deviatoric stress condition ( Fig. 14c). Recently, Lei et al. [280] also studied the impacts of stress on reorganising vertical fluid flow through a 3D sedimentary layer.

Permeability
There are two different notions of rock mass permeability, i.e. equivalent permeability and effective permeability. The equivalent permeability is defined as a constant tensor in Darcy's law to represent flow in a heterogeneous medium, while the effective permeability is an intrinsic material property based on the existence of an REV at a large homogenisation scale [298]. Permeability here mainly refers to the equivalent permeability of a fractured rock at a specific study scale.
The permeability tensor was observed to be highly dependent on both the geometrical attributes of fracture networks (e.g. density, lengths and orientations) and the in-situ stress conditions (e.g. direction, magnitude and ratio of the principal stresses) [45]. When the differential stress ratio is relatively low, the permeability decreases with the increase of burial depth (or mean stress) of the fractured rock due to the closure of most fractures [43,206]. The non-linear relationship between normal stress and fracture closure results in a phenomenon that the permeability is more sensitive at shallower depths (i.e. smaller mean stresses) and approaches a minimum value when most fractures are closed to their residual apertures under high mean stresses (Fig. 15a) [206]. The permeability anisotropy of a fracture network with non-systematic fractures is more dependent on the ratio and direction of the applied principal stresses than that of a network with systematic fracture sets which is more controlled by the fracture set orientations [43]. With the increase of differential stresses, the permeability exhibits a decrease and then an abrupt increase separated by a critical stress ratio that begins to cause continued shear dilations along some preferentially oriented fractures (Fig. 15b) [206]. A similar variation of permeability occurs when the pore fluid pressure is elevated [41]. Simultaneously, the permeability anisotropy is also enlarged by the increased stress ratio [206]. If initial apertures are correlated with fracture lengths, the permeability of fractured rocks is dominated by larger fractures with wider apertures. This model tends to exhibit a permeability value much higher than the model assigned with a constant initial aperture (Fig. 15c) [98]. The permeability tensor may be destroyed but then reestablished with the increase of differential stresses [98]. In addition to the shear dilation of fractures, the increase of network connectivity caused by brittle failure and crack propagation under geomechanical loading can also significantly raise the permeability of fractured rocks [15,41,151,157]. More interestingly, the emergence of dilational jogs/bends associated with curved fractures in response to high differential stresses may lead to even more increase of permeability in the third dimension [47,96,280] and result in a highly anisotropic permeability tensor [278]. The permeability of 3D fracture networks also exhibits a transition stage having steep growth that occurs when the stress ratio is approaching the critical state (Fig. 15d) [278]. It has to be mentioned that an increased fracture density may not always lead to an increased permeability under some tectonic conditions (e.g. an extensional regime), because fractures may be closed due to mechanical interactions when they are too densely spaced [299].

Transport
Mass transport in fractured rocks is governed by various mechanisms including advection, dispersion, matrix diffusion, interface sorption and chemical reaction [17,297]. The heterogeneous fluid velocity fields in geological media consisting of distributed fractures and porous rocks can result in complex transport phenomena in the system. Computational models employing solute components or tracked particles have been developed to simulate migration processes in fractured rocks [6,300]. Recently, numerical studies have also been conducted to investigate the effects of the stress/deformation on the transport properties of fracture networks, as summarised below.
Zhao et al. [208] coupled the UDEC DEM code and a random walk particle tracking code PTFR and investigated the stress effects on the hydrodynamic dispersion of contaminant solutes in a stochastic DFN system. They found that compressive stresses can close fracture apertures and attenuate the dispersivity, but an increased differential stress ratio could greatly intensify the spreading phenomenon if it exceeds certain threshold for triggering shear dilations. Zhao et al. [209] extended this modelling technique to incorporate the effects of matrix diffusion and sorption, and conducted a systematic study of the solute transport under various stress conditions. Their results showed that the stress can significantly affect the solute residence time in the fracture network (Fig. 16a). When the stress ratio is increased but not very high (<3), the breakthrough curve shifts to the right side (i.e. the average residence time increases) compared to that of the initial zero-stress condition, due to the closure of fractures under relatively isotropic stresses. However, as the stress ratio exceeds 3, fluid velocity is raised drastically in some dominant channels formed by dilated fractures due to shearing, and thus the residence time decreases with the breakthrough curve shifting backward. They also observed that the breakthrough curve for interacting tracers (i.e. with matrix diffusion) exhibits longer tails than that for non-interacting tracers (i.e. without matrix diffusion) due to the meandering of a small amount of particles passing through tortuous fluid pathways. Such long tail phenomena were more significant when the pressure gradient is small, for which the matrix diffusion tends to play a dominant role in solute transport (Fig. 16b). Rutqvist et al. [101] used an extended multiple interacting continua model combined with the crack tensor approach to simulate the advection-dominated transport (under a high hydraulic gradient) and diffusion-retarded transport (under a low hydraulic gradient). In addition to the stress-dependent transport  [206]. (b) The vertical fluid flow through a jointed layer exhibits a highly localised pattern when the fractured rock is deformed under a critical stress state [47]. (c) The vertical flow in a 3D persistent fracture network shifts from a homogeneous to localised pattern with the increase of the horizontal stress ratio [278].  [206]. (b) Permeability change with the increase of stress ratio for a DFN with a constant initial aperture [206]. (c) Permeability change with the increase of stress ratio for a DFN with a lognormally distributed and length correlated apertures under rotated stress fields [98]. (d) Increase of vertical permeability with the increase of horizontal stress ratio for a 3D persistent fracture network [278]. behaviour, they also observed a delayed breakthrough in low pressure gradient scenarios due to the residence of solutes in the porous rock matrix. Wang et al. [301] further applied this multicontinuum method to demonstrate the contribution of inactive fractures (i.e. isolated cracks or dead-ends of fractures) for stagnating solutes by providing additional surface areas for diffusive transfer into/out of matrix pores. Zhao et al. [102] compared the stress-flow models of five different research groups, which showed consistency in predicting the stress dependency of mass transport in a fractured rock. Recently, Kang et al. [302] studied the impact of in-situ stresses and aperture changes on the anomalous transport in a fractured rock based on an outcrop-based DFN model. Apart from the stress-induced aperture change that can affect the transport behaviour, Nick et al. [162] found that the propagation of fractures under tectonic loading can also vary the breakthrough properties due to the increased fracture density and network connectivity.

Comments on the use of DFNs for HM modelling
The geomechanical effects on fluid flow in fractured rocks as observed in numerical simulations demonstrate the importance of using appropriate fracture network representations and conducting systematic geomechanical computations for characterising the coupled HM behaviour of fractured geological formations. Simulation results rely strongly on the suitability and accuracy of the constructed DFNs and the geomechanical models of rocks and fractures.
Natural fracture networks have been observed to exhibit fractal characteristics and long-range correlations [8,31,32,73,124-127,1 30]. The fractal (or clustered) characteristics have been found to have significant impacts on the connectivity, permeability and flow pattern of fracture networks [131][132][133]. The scattered distribution of fracture sizes (e.g. power law distribution) also considerably affects the hydrological [92][93][94], geomechanical [134] and HM [40] behaviour of fractured rocks. These important features may need to be prudently incorporated when creating DFNs. The hydraulically inactive fractures (i.e. dead-end and isolated fractures) were found playing important roles in propagating new fractures that can link with pre-existing fractures [15,39,41,151,161,280], and in transferring fluids and solutes into/out of matrix pores [101,171,301]. Thus, they may need to be preserved in DFNs for HM modelling of fractured rocks. The bent fractures that can accommodate large shear dilations and induce localised fluid flows [15,47] may also need to be adequately characterised in DFN representations. Such curvature characteristics may provide additional surface areas for diffusive transfer and are therefore important for solute transport modelling. Special attention is needed when defining fracture apertures in DFNs for HM modelling. Fracture apertures have been suggested to have intrinsic correlations with fracture sizes [5,80], which tend to exhibit a significant power law scaling over a broad spectrum of length scales [8,73]. This long-range correlation has been found to have critical influences on the permeability and flow distribution of fracture networks [93,94,97,99] as well as their stress dependency [40,98]. In addition, the geomechanical simulation tool also needs to be carefully adopted or developed for HM modelling. The following aspects may be of concern: (1) whether the elastic/elastoplastic deformation of rock matrix can be modelled; (2) whether the non-linear deformation of natural fractures (opening/closing, shearing and dilation) can be captured; (3) whether the fracture propagation can be simulated; and (4) whether the geomechanical model can be coupled with fluid flow solvers.
It can be noted that this section is focused on the research of modelling geomechanical impacts on the hydrological behaviour of fractured rocks. The reverse effect of fluid flow on rock mass deformations is also a very important issue which is relevant to many engineering problems such as hydraulic fracturing, CO2 sequestration and reservoir compaction. Recent developments in exploitation of shale gas resources have highlighted the demand for a better understanding of the interaction between hydraulic fractures and natural fractures. Several studies have been done to develop coupled HM models to simulate fluid-driven fracture propagations in pre-existing DFN networks [185,186,[303][304][305][306][307][308].

Outstanding issues and concluding remarks
Modelling the geomechanical evolution and the resulting hydrological characteristics of fractured rocks is a challenging issue. The previous HM simulation results have shown consistency with field measurements: (i) only a small portion of fractures are conductive [6,144], (ii) permeability is less sensitive to depth in deep rocks [2], and (iii) critically stressed faults tend to have much higher hydraulic conductivity [16,309]. However, a number of important issues of difficulty still exist in applying DFNs for HM modelling of fractured rocks. The first task will be the development of realistic DFNs that can represent the important geometrical and geostatistical characteristics of natural fracture systems (e.g. clustering, scaling, connectivity, intersection, termination, curvature and correlation), and that can further be used to capture the key geomechanical and hydrological characteristics of fractured rocks (e.g. strength, deformation, permeability and mass transport). The geometry of the created DFNs needs to be consistent with the geomechanics such as the in-situ stress field so that the geometrically-and geomechanically-dependent fluid flow can be properly modelled. Another critical issue is to develop appropriate upscaling approaches to DFN models to evaluate the large-scale behaviour, and this would require the preservation of geostatistical and geomechanical characteristics. Some techniques have been proposed to construct large-scale heterogeneous models based on upscaled stress-dependent permeability tensors of local grid blocks [310][311][312] or using self-referencing random walks associated with stress-dependent aperture information [40]. More efforts are also needed in the future for many other aspects, such as developing more advanced ''two-way" coupling schemes, modelling geomechanical effects on multiphase flow and, importantly, extension to complex 3D systems.
The development of calibration and validation procedures for DFN models based on field measurements is a crucial issue if the numerical models are to be used for practical applications. The DFN model can be examined by comparing the geometrical properties of a simulated 2D trace network with those of the trace pattern exposed on an outcrop or excavation wall [117,138,141]. Mine-by experiments in underground research laboratories [313,314] provide also valuable patterns of fracture traces that intersect tunnel walls and measurement data of rock mass mechanical responses (e.g. breakout characteristics, spalling depth, tunnel radial deformation, etc.). Such data might be very useful for the DFN calibration study. Hydraulic test data such as groundwater inflow recorded at boreholes/shafts, pressure variations and derivatives obtained from well/interference tests, and particle concentration monitored during tracer tests can also be used to tune the hydrological performance of DFN models [138][139][140][142][143][144]. Furthermore, the oil recovery case history information in reservoir engineering [315] offers useful references for the calibration of DFN models used for multiphase flow simulation. Recently, an inverse model constrained by in-situ hydraulic data from sequential pumping tests has been developed to derive the spatial distribution of flow channels on the bedding plane of a karstic formation [316]. This new technique can be employed to calculate best-fit aperture distributions for DFN models based on an optimisation process using multiple objective functions. Other possible calibration approaches may include those based on heat transfer experiments and electrical resistivity measurements of fractured rocks.
In summary, the present paper began by presenting an overview of various discrete fracture network (DFN) models for simulating the geometry and topology of natural discontinuity systems. Different continuum and discontinuum approaches that integrate DFN geometries to model the geomechanical behaviour of fractured rocks were then surveyed. Numerical results of the fracture-dependent mechanical response and stress-dependent hydrological characteristics of fractured rocks suggest that it is important to use appropriate DFN representations and conduct geomechanical computations to better characterise the bulk behaviour (e.g. strength, deformation, permeability and mass transport) of highly disordered geological media embedded with naturally occurring discontinuities. Some of the outstanding issues that need to be addressed in the near future will be those related to generation of more realistic DFNs, improvement of computational efficiency, study of coupled multiphysics processes, calibration and validation of numerical models as well as upscaling for largescale practical applications.