Implementation of an Empirical Joint Constitutive Model into Finite-Discrete Element Analysis of the Geomechanical Behaviour of Fractured Rocks

An empirical joint constitutive model (JCM) that captures the rough wall interaction behaviour of individual fractures associated with roughness characteristics observed in laboratory experiments is combined with the solid mechanical model of the finite-discrete element method (FEMDEM). The combined JCM-FEMDEM formulation gives realistic fracture behaviour with respect to shear strength, normal closure, and shear dilatancy and includes the recognition of fracture length influence as seen in experiments. The validity of the numerical model is demonstrated by a comparison with the experimentally established empirical solutions. A 2D plane strain geomechanical simulation is conducted using an outcrop-based naturally fractured rock model with far-field stresses loaded in two consecutive phases, i.e. take-up of isotropic stresses and imposition of two deviatoric stress conditions. The modelled behaviour of natural fractures in response to various stress conditions illustrates a range of realistic behaviour including closure, opening, shearing, dilatancy, and new crack propagation. With the increase in stress ratio, significant deformation enhancement occurs in the vicinity of fracture tips, intersections, and bends, where large apertures can be generated. The JCM-FEMDEM model is also compared with conventional approaches that neglect the scale dependency of joint properties or the roughness-induced additional frictional resistance. The results of this paper have important implications for understanding the geomechanical behaviour of fractured rocks in various engineering activities.


Introduction
Fractures are ubiquitous in crustal rocks in the form of faults, joints, and veins over different length scales . Conceptually, there are two physical domains involved: fractures with relatively low stiffness and high porosity, and rock matrix with the opposite characteristics. The hydromechanical properties of fractured rocks are often governed by both the behaviour of individual fractures and the interactions among the fracture population (Zimmerman and Main 2004). The dominant role of fractures in crustal processes, such as in localising deformation and fluid flow (Sibson 1994), has important implications for various engineering applications including hydrocarbon production, geothermal energy, groundwater management, nuclear repository safety, and ground engineering.
To describe the behaviour of individual fractures associated with intrinsic surface roughness, many studies have been reported in the past few decades. Goodman (1976) proposed a hyperbolic relation to characterise the nonlinear closure of fractures under normal compression and studied the effect of mismatch between opposite rough joint walls. Barton and Choubey (1977) introduced an empirical system based on three main index parameters, i.e. joint roughness coefficient (JRC), joint wall compressive strength (JCS), and residual friction angle, to predict the shear strength of natural fractures. These parameters can be measured in the laboratory from tilt tests or shear box experiments. Bandis et al. (1983) summarised a series of empirical equations to interpret the deformation characteristics of rock joints in normal loading and direct shear experiments. Size effect on shear strength and deformation characteristics of individual fractures were further investigated based on the laboratory experiment results conducted on natural fracture replicas that were cast at different sizes (Bandis 1980;Bandis et al. 1981;Barton 1981). Fractures having the same roughness morphology but different sizes may exhibit distinctly different mechanical responses (Barton and Bandis 1980;Barton 2013). The empirical joint constitutive model was recently improved to capture the stress dependency of peak shear displacement (Asadollahi and Tonon 2010). In addition to the laboratory experiments, numerical studies have also been performed to model the behaviour of single fractures with the roughness profile represented explicitly. Karami and Stead (2008) used a finite-discrete element code to simulate the shearing process of rough joint specimens. Joint shear behaviour, e.g. peak shear displacement, shear strength, dilational displacement, and asperity degradation, was found to be dominated by the roughness geometry and the applied normal stress. Bahaaddini et al. (2014) studied the scale effects on joint shear strength and peak shear displacement using a particle-based discrete element model. Their numerical results reveal that with the increase in joint length, the peak shear strength and shear stiffness decrease, while the peak shear displacement increases, similar to the phenomena observed in the laboratory experiments. Grasselli (2012, 2015a) investigated the scale-dependent shear strength and the asperity damage mechanism of rough fractures in direct shear tests, using a numerical modelling approach that is based on a finite-discrete element model comprehensively calibrated for tensile and shear fracturing and employs explicit representation of fracture surfaces. The numerical solutions of shear strength and dilation were further compared with the experimental results from micro-CT imaging (Tatone and Grasselli 2015a).
The overall behaviour of fractured rocks involving multiple interconnected fractures has also been widely investigated based on various numerical methods. Sanderson and Zhang (1999) used the discrete element method (DEM) to analyse the deformation of naturally fractured rocks, and they found very large apertures are generated when the fracture system is under a critical stress state. Min et al. (2004) applied the DEM model to analyse the influence of far-field stresses on the distribution of fracture apertures with the effects of nonlinear normal deformation and shear dilation considered. Harthong et al. (2012) combined the particle-based DEM method with discrete fracture networks (DFNs) to study the strength characteristic of pre-fractured rocks under triaxial stress loading. Latham et al. (2013) employed the finite-discrete element method (FEMDEM) to simulate the geomechanical response of a natural fracture system with the consideration of bent fractures that accommodate high dilation, and crack propagation that can connect pre-existing fractures and increase network connectivity. Lei et al. (2014) examined the stress effect on the validity of synthetic networks for representing a two-dimensional (2D) naturally fractured rock in terms of geomechanical and hydraulic properties. The geomechanically induced fracture apertures captured by the 2D FEMDEM model were further scaled up for larger scale modelling of fluid flow in a self-referencing multi-scale system of fracture networks (Lei et al. 2015a). Some preliminary work on modelling the geomechanical behaviour of three-dimensional (3D) fracture networks has also been conducted using the 3D FEMDEM solver (Lei et al. 2015b. In summary, previous research has focused on two separate scales at which fracture properties affect rock mass strength and deformation: (1) the level of the individual fracture where surface roughness is represented in detail and (2) the level of fracture network with emphasis on the overall properties. For more realistic modelling of fractured rocks, it is of importance to integrate the detailed characteristics of individual fractures into the simulation of complex fracture systems (Jing and Stephansson 2007). In the past few decades, a few attempts have been made to bridge these two scales in numerical modelling (Saeb and Amadei 1992;Cai and Horii 1992;Jing et al. 1994). Some pioneer work has also been done by Mahabadi and Grasselli (2010) to implement an empirical shear strength criterion into the FEMDEM model. However, to better characterise the nonlinear deformation of rock fractures, some important aspects, such as the mobilisation of friction coefficient during shearing processes and the scale dependency of fracture roughness properties, may need to be more adequately considered.
The objective of this study is to develop a methodology to incorporate realistic joint constitutive characteristics in the numerical simulation of fractured rock masses involving pre-existing and propagating fractures. A novel feature of the proposed method is its capability of simulating the important size effect of fracture wall properties as observed in laboratory experiments through a systematic characterisation of fracture network topologies. The paper is organised as follows. In Sect. 2, approaches for mechanical modelling of multi-body systems are briefly discussed covering the issues of solid deformation, interaction, and fracture propagation. The constitutive models for rock fractures are reviewed with respect to normal closure, shear deformation, and dilatancy. A scheme to couple the empirical joint constitutive model (JCM) with the FEMDEM computation is then formulated. Section 3 presents a verification of the proposed numerical method. In Sect. 4, a natural fracture pattern involving intersections, bends, and roughness-induced initial apertures is presented. Numerical experiments are designed to illustrate various in situ stress conditions for which effects of far-field stresses on fracture system properties are investigated. A comparison with conventional joint modelling approaches is also presented. Finally, a brief discussion is given, and conclusions are drawn.

Governing Equation
The motions of elements are controlled by the forces acting on elemental nodes, and the governing equation is given by (Munjiza 2004;Xiang et al. 2009): where M is the lumped nodal mass matrix, x is the vector of nodal displacements, f int are the internal nodal forces induced by the deformation of triangular elements, f ext are the external nodal forces including external loads f l contributed by boundary conditions and body forces, cohesive bonding forces f b caused by the deformation of cohesive joint elements, and contact forces f c generated by the contact interaction via broken joint elements. The FEMDEM solid model is capable of modelling both the deformation and interaction of matrix bodies under prescribed boundary conditions (Munjiza et al. 2011(Munjiza et al. , 2015. The equations of motion of the FEMDEM system are solved by an explicit time integration scheme based on the forward Euler method.

Contact Force
Contact force between two discrete solids (one is named contactor and another target) is computed based on the penalty function method by integration over the boundary of penetration (Munjiza and Andrews 2000): Implementation of an Empirical Joint Constitutive Model into Finite-Discrete Element Analysis… 4801 where n is the outward unit normal to the penetration boundary C c , while u c and u t are potential functions for the contactor and target solids, respectively. In numerical implementation, the total contact force between two discrete solids is calculated as the summation of contact force between a set of couples of interacting finite elements. Interaction between two finite elements is further reduced into interactions between the contactor and the edges of target element. The normal contact force f n and tangential friction force f t exerted by a contactor onto a target edge are given by (Munjiza 2004;Munjiza et al. 2011): where L p is the penetration length, u is the potential function along the target edge, v r is the relative velocity (at the Gauss point) between the contactor and the target edge, p is the penalty term, and l mob is the mobilised friction coefficient which varies during the shearing process. No effort is required here to classify static and dynamic friction stages since they are inherently modelled by the mobilised friction coefficient (further details are given in Sect. 2.3.3).

Crack Propagation
Crack propagation induced by stress concentration is modelled by a smeared crack model (Munjiza et al. 1999) embedded in the FEMDEM formulation with both mode I and mode II brittle failure captured (Latham et al. 2013). Fracture initiation and propagation is characterised as occurring in three stages: (1) the continuum stage simulated by the finite element method through the solid constitutive law, (2) the transition stage described by the strain softening using the smeared crack model, and (3) the discontinuum stage in which elements along the new crack are physically separated with their interaction further modelled by the discrete element method (Munjiza 2004).

Joint Normal Deformation
Based on laboratory experiments, rock joints were found to exhibit nonlinear deformation response under compressive normal stress (Goodman 1976). An empirical hyperbolic model was proposed by Bandis et al. (1983) to represent this nonlinear relation: or where v n is the current closure (mm) under the normal stress r n (MPa), k n0 is the initial normal stiffness (MPa/ mm), and v m is the maximum allowable closure (mm).
Values of k n0 and v m are given by (Bandis et al. 1983): where a 0 is the initial aperture (mm), JRC is the joint roughness coefficient, JCS is the joint compressive strength (MPa). Coefficients derived from experimental measurements of numerous joint samples of five different rock types under a third loading cycle are adopted since in situ fractures are considered more likely to behave in a manner similar to the third or fourth cycle (Barton et al. 1985).
Both JRC and JCS are scale-dependent parameters (Bandis et al. 1981;Barton 1981), and their values in field scale, i.e. JRC n and JCS n , can be estimated using (Barton et al. 1985): where L n is the effective joint length (i.e. size of a block edge between fracture intersections) defined by the spacing of cross-joints, JRC 0 and JCS 0 are measured based on the laboratory sample with length L 0 . For the laboratory sample, the initial aperture a 0 may be estimated using an empirical relation (Bandis et al. 1983) as given by: where r c is the uniaxial compressive strength (MPa), and JCS 0 (MPa) can be set equal to r c , assuming the effect of weathering can be ignored. Under a varying normal stress condition, the joint normal stiffness k nn is given by (Saeb and Amadei 1992):

Joint Shear Deformation
Peak shear strength s p of fractures under different normal stress levels can be calculated by the following empirical law of friction (Barton and Choubey 1977): where r n is the normal compressive stress (MPa), and / r is the residual friction angle. The shear stress-displacement curve of rock joints in direct shear experiments shows two major phases, i.e. pre-peak and post-peak stages. Such relation can be empirically characterised by replacing JRC n in Eq. (12) with the mobilised value JRC mob : where s is the current shear stress, JRC mob can be calculated using the dimensionless model (Barton et al. 1985) as shown in Table 1, in which u is the current shear displacement, and u p is the peak shear displacement. The scale dependency of peak shear displacement u p can be characterised by (Barton et al. 1985): which was modified by Asadollahi and Tonon (2010) to further consider its stress dependency as given by For post-peak stage, JRC mob can also be estimated using a power-base empirical relation given by (Asadollahi and Tonon 2010): The joint shear stiffness k tt can be derived as the slope of the shear stress-shear displacement curve:

Joint Shear Dilatancy
During the shearing process under a normal stress, fractures contract first due to the compressibility of asperities and then dilate with roughness damaged and destroyed. Dilational displacement can be related to the shear displacement based on an incremental formulation given by (Olsson and Barton 2001): where dv s is the increment of normal displacement caused by shear dilation, du is the increment of shear displacement, and d mob is the mobilised tangential dilation angle. A quadratic equation was proposed to describe the prepeak dilational displacement with the tangential dilation angle given by (Asadollahi and Tonon 2010): where v p is the normal dilational displacement corresponding to the peak shear displacement u p and can be calculated from (Barton and Choubey 1977): where M is a damage coefficient that is determined by (Barton and Choubey 1977): For the post-peak phase, surface asperities of fracture walls begin to be damaged as shearing continues, and the variation in the tangential dilation angle can be captured by (Olsson and Barton 2001):

Coupled Joint Normal and Shear Behaviour
Fractures in crustal environment may experience complicated loading paths, e.g. shearing under a variable normal stress (Saeb and Amadei 1992). By combining Eqs. (11) and (18), the coupled behaviour of normal and shear  deformation can be modelled by an incremental formulation given as: or after rearrangement: It can also be written in a more compact form as: where k nn and k nt are the corresponding normal stiffness coefficients. A similar equation can be expressed for the relation between the increments of shear stress and displacement components: where the stiffness coefficient k tn is commonly assumed to be zero (Jing and Stephansson 2007) and k tt is derived using Eq. (17). A differential formulation for the rock joint deformability can be further expressed by a non-symmetric material tangent stiffness matrix as follows: 2.3 Combined JCM-FEMDEM Formulation

Joint Element
The combined fracture-matrix solid system (Fig. 1a) is represented by a discontinuous discretisation of the modelling domain using three-noded triangular elements and four-noded joint elements embedded between the edges of triangular elements (Fig. 1b). The deformation of the bulk material is captured by the linear elastic triangular finite elements with impenetrability enforced by a penalty function and continuity constrained by the constitutive relation of cohesive (unbroken) joint elements (Munjiza et al. 1999), while the interaction of matrix bodies through discontinuity interfaces is simulated by the penetration calculation along fracture (broken) joint elements (Munjiza and Andrews 2000). Construction of cohesive joint elements is achieved by a detachment algorithm based on the original continuous triangular mesh of the matrix domain, whereas formation of fracture joint elements is realised based on the initial overlapping configuration of the edges of opposite triangular elements along pre-existing fractures. Fracture initiation and propagation due to stress concentration is modelled by the transition from cohesive joint elements to fracture joint elements.
The mesoscopic (i.e. at the scale of grid discretisation) local normal and tangential displacements of a pre-existing Fig. 1 Representation of a a fracture-matrix system using a mesh consisting of b threenoded triangular elements and four-noded cohesive/fracture joint elements embedded between edges of triangular elements. c Fracture local opening and shearing displacements characterised by the geometrical deformation of the nodal system of joint elements fracture can be captured by the geometrical deformation of the nodal system of joint elements. As shown in Fig. 1c, deformation of the joint element AB-A 0 B 0 can be measured by a vector of coordinate difference between the midpoints (i.e. C and C') of the opposite edges, given bŷ based on which a local orthogonal coordinate system can be established with mutually unit base vectors defined bŷ Thus, the mesoscopic opening displacement w and shear displacement u can be calculated as Fracture aperture a m is further derived by combing the effects of mesoscopic opening and microscopic closure as given by where v is the accumulative closure derived from the incremental formulation, i.e. Eq. (26). The first part of the piecewise function corresponds to the scenario that the fracture joint element is mesoscopically opened, while the second part models the local two fracture walls are in contact at the scale of FEMDEM discretisation.

Characterisation of Fracture Systems Based on a Binary-Tree Search
Due to the scale dependency of fracture parameters such as JRC, JCS, and peak shear displacement (as discussed in Sect. 2.2), it is important to precisely characterise the distribution of effective fracture lengths (i.e. size of a block edge between fracture intersections) in the numerical modelling of a disordered, interconnected fracture system. One critical numerical difficulty related to effective fracture lengths is to distinguish the sophisticated topological relations of what is very often a complex system containing numerous joint elements, in which some pre-existing fracture joint elements may connect with each other to form a continuous fracture wall (i.e. block edge) and would act together as an equivalent individual fracture with two facing walls. A generic algorithm has been developed in this research for the topological diagnosis of general fracture networks involving bends, intersections, termination, and impersistence. Connectivity analysis is first implemented to recognise neighbours of each joint element based on the initial geometric coordinates, in which a joint element connecting the model boundary or a fracture intersection is considered having no neighbour on that side with a '-1' value assigned numerically, as shown by the schematic example in Fig. 2a. Binary-tree structures are constructed with the tree nodes representing joint elements (Fig. 2b). When scanning through the binary-tree system, previously visited tree nodes or unreal neighbour tree nodes are labelled to be dead (empty nodes in Fig. 2b) and will not grow in further loops. Block edges are identified as the connected chains of live tree nodes (solid nodes in Fig. 2b). Thus, scale effect of fractures can be modelled by relating the constitutive parameters of each local joint element to the length of its corresponding block edge.

Coupling Between JCM and FEMDEM
The JCM and FEMDEM modules are combined to achieve compatibility with respect to both stress and displacement fields. The displacement fields of JCM and FEMDEM are linked through Eq. (31), while the stress fields are coupled in both normal and tangential directions along the fracture interface. Normal stress of a joint element is extracted from adjacent finite elements of the FEMDEM solid model using where r is the Cauchy stress tensor of the finite element located on the opposite fracture walls, and n = [n x ,n y ] T is the outward unit normal vector of the finite element edge. By substituting the incremental value of normal stress and shear displacement into the JCM formulation, i.e. Eq. (26), the incremental normal displacement can be solved with the aperture further derived by Eq. (31). Friction angle between two rough fracture walls is often larger than the residual friction angle due to the effect of asperities (Barton and Choubey 1977). The friction coefficient also varies during the progression of shearing as a result of roughness degradation (Olsson and Barton 2001). Mobilised friction coefficient l mob of each fracture joint element can be calculated using its current parameters by: The updated friction coefficient is transferred to the FEMDEM solver in each time step for calculation of the Implementation of an Empirical Joint Constitutive Model into Finite-Discrete Element Analysis… 4805 tangential friction force between a contactor and a target edge as given by Eq. (4).

Numerical Verification
The empirical constitutive laws are implemented in the FEMDEM framework at the joint element scale, but the consistency between the simulated macroscopic fracture behaviour and the empirical formulations requires a detailed verification. The consistency between the empirical formulations and the laboratory experiments has been well demonstrated in the literature (Bandis et al. 1981(Bandis et al. , 1983Barton et al. 1985;Olsson and Barton 2001). Hence, the validity of the numerical model will be examined by comparing numerical results with the empirical solutions, i.e. Eq. (13) for the shear stress, the integral of Eq. (18) for the dilational displacement, and Eq. (5b) for the normal closure.
The model set-up is based on the physical experiment conducted by Bandis (1980) which used a series of cast replicas of natural joint surfaces prepared in different sizes, i.e. 6, 12, 18, and 36 cm. The material used for casting joints in the laboratory was made from the mixture of silver sand, alumina, barites, and water. The density of the analogue material q was 1850 kg/m 3 , and the Young's modulus E was 0.8 GPa. The Poisson's ratio was not provided in the reference (Bandis 1980), so a typical value of 0.3 is assumed for the numerical model. As shown in Fig. 3, the specimen consists of an upper portion and a longer lower portion, and is placed in a shear box made of steel having a density of 8030 kg/m 3 , a Young's modulus of 190 GPa, and a Poisson's ratio of 0.3. The bottom and right sides of the lower steel shell are constrained by the roller boundary conditions, while the upper one is free to move. The normal stressr n applied on the top of the shear box is designed to generate a constant normal stress r n : 24.5 kPa on the joint surface with the consideration of the gravitational forces of the upper block and the steel shell. The shearing of the two fracture walls is controlled by the velocity boundary condition applied on the upper half of the shear box. The input joint properties for the numerical models of different-sized joints were based on the smallest sample, i.e. L 0 = 6 cm, JRC 0 = 15.0, JCS 0 = r c = 2 MPa, and / r = 32°(the properties of the larger joints will be scaled up using Eqs. (8) and (9) based on their actual lengths Fig. 2 Characterisation of a fracture system, in which four block edges from two intersecting fractures are discretised into a number of joint elements, based on a connectivity analysis and b binary-tree search identified by the algorithm as described in Sect. 2.3.2). The penalty term p for the specimen is chosen to be 20 times that of the Young's modulus (Mahabadi 2012), i.e. p = 16 GPa. The damping coefficient g is assigned to be the theoretical critical value, i.e. g = 2 h (Eq) 1/2 , where h is the element size, to reduce dynamic oscillations. The numerical shear stress is derived as the quotient between the total tangential contact force integrated for all upper wall nodes and the length of the joint sample. In contrast to both the indirect measurement method which is used in laboratory testing of shear strength (i.e. by monitoring the horizontal forces loaded on the shear box in the laboratory) (Bandis 1980) and the method adopted for the numerical modelling of an explicit roughness profile (Karami and Stead 2008;Bahaaddini et al. 2014;Grasselli 2012, 2015a), for the verification of the proposed JCM-FEMDEM framework, the tangential force acting on the joint surface is directly extracted from the contact algorithm and emerges by virtue of the forces recorded by the joint element data structure. It also gives an unbiased measurement of the joint frictional forces. To mimic the quasi-static loading condition in the physical test, a staged loading scheme is applied, which is similar to the one used by Kazerani et al. (2012) for modelling uniaxial compression and Brazilian tests. In each stage, the upper block moves for 0.01 mm (say taking N time steps). Then, the block stops (zero shearing velocity) and the FEMDEM calculation is cycled to approach equilibrium by running for anther N time steps.
The sensitivity of the shear stress-shear displacement behaviour to the loading velocity is shown in Fig. 4a. The numerical models discretised by the same very fine mesh with an element size of 1 mm along the joint are loaded by different shearing velocities ranging from 1 to 5 mm/s. It can be seen that, as the velocity decreases, the oscillation of the modelling results is dramatically reduced and the numerical curve gradually approaches the empirical one. A velocity of 1 mm/s is considered to be an adequately small rate for simulating the quasi-static condition, which however requires a run-time of approximately 100 h on a desktop computer equipped with an Intel Core E5-1620@3.70 GHz. The effect of mesh size on the shear stress is assessed by comparing the modelling results with different element sizes along the joint, i.e. 1.5, 1.25, and 1 mm (Fig. 4b) under the same shearing velocity of 1 mm/ s. With the refinement of the mesh, the numerical result also gradually converges to the target empirical solution.
Based on the results of the sensitivity analysis, the models for larger joint sizes, i.e. 12, 18, and 36 cm, are discretised with an element size of 1 mm along the joint, and further sheared under the same loading velocity of 1 mm/s and the same constant normal stress of 24.5 kPa. The similarity between the empirical and numerical curves of shear stress-shear displacement (Fig. 5a) verifies the implementation of the JCM-FEMDEM model. The numerical predictions for the joint dilational behaviour also fit well to the empirical values (Fig. 5b). During the shearing process, the joint specimens exhibit a certain contraction in the pre-peak stage and a considerable dilation in the post-peak stage. It is reassuring that the scale effects on joint shearing behaviour observed in the laboratory test have been well captured by the JCM-FEMDEM model. With the increase in the joint sample size, the value of peak shear displacement increases, a transition from a 'brittle' to 'plastic' shear failure mode occurs, and a higher dilational displacement is generated.
In order to also examine the numerical model with respect to normal closure, the 6-cm joint sample is loaded with a normal stress gradually increased up to a value of 1 MPa, which is still smaller than the uniaxial compressive strength (i.e. r c = 2 MPa) of the analogue material and therefore will not cause breakage in the intact blocks. No shearing condition is imposed for this test of normal closure. As shown in Fig. 6, the numerical model gives consistent results with the empirically calculated values.
To sum up, the consistency of the numerical results with the empirical solutions demonstrates the performance of the combined JCM-FEMDEM formulation for capturing realistic shear strength and normal closure behaviour of single fractures, although it is recognised that it would be ideal to further test the model over a parameter space with different JRC, JCS, normal stresses, etc. In the following section, the numerical model will be applied to simulate the geomechanical behaviour of a complex fracture network under in situ stresses.  (Belayneh et al. 2009). The fracture pattern possesses a connectivity state of interest that is near the geometric percolation threshold under an intermediate fracture density (Lei et al. 2014). In this geological site, two systematic sets of vertical, layer-normal fractures were formed extensionally and filled with calcite minerals, striking approximately 100°(Set 1) and 140°(Set 2), respectively. The fractured limestone layer (*26 cm thick) is sandwiched by almost impervious shales, and the vein sets are layer bound. 2D plane strain analysis is used to capture the geomechanical response of the fractured rock under in situ stresses. Since the currently available processing power means it is very expensive in CPU time to compute very large domains, a 2 m 9 2 m subarea is extracted as a sample of the fracture system in the 2D study. Material properties for intact rock and fractures of the limestone layer are assumed as given in Table 2. The properties of intact rocks (e.g. density, Young's modulus, and strength properties) are generalised ones that are assumed to represent a typical type of limestone, according to the range of limestones found in the handbook by Lama and Vutukuri (1978). The mode I energy release rate G I is selected to be 100 J/m 2 and corresponds to a fracture toughness K IC of 1.8 MPa/m -1/2 for the plain strain problem, which is in the typical range for rocks, i.e. 1-5 MPa/ m -1/2 (Atkinson and Meredith 1987). The mode II energy release rate G II is chosen to be 500 J/m 2 , which also lies in the typical range for rocks (Cox and Scholz 1985). The parameters for the rough fractures are based on the reference of Olsson and Barton (2001). The potential correlation between intact rock strength properties and the energy release rates is not considered here for simplicity, but can be calibrated following the procedure developed by Tatone and Grasselli (2015b) when solving real engineering problems. Alternatively, one can find very useful empirical correlations between K IC and the physical properties and strength measures from the literature, such as Gunsallus and Kulhawy (1984).
In the geological setting, fractures often exhibit displacements perpendicular and/or parallel to the discontinuity surface, i.e. aperture and shear displacement. In this study, fractures are represented with no initial phase of shearing before the phases of far-field stress application. However, an initial aperture is considered significant and is assigned a priori to all fractures equally to enable the introduction of a potentially realistic joint aperture that is in turn controlled by the roughness characteristic. This aperture is assigned a value based on the empirical relation (Bandis et al. 1983) given by Eq. (10). Fractures are modelled as uncemented interfaces, i.e. they have no cohesion or tensile strength.

Model Set-Up and Simulation Results
The fractured rock is considered to be at a depth of *300 m with a pore fluid pressure ratio (i.e. the ratio of pore fluid pressure to lithostatic stress) equal to 0.35, producing an overburden effective stress of 5 MPa. Mechanical response of the rock mass is investigated based on a series of plane strain numerical experiments with biaxial effective stresses applied to the square-shaped model boundaries. Effect of pore fluid pressure is assumed here to be a second-order factor for aperture development and is not included in the simulation. Poroelastic effect of the Biot-type coupling of pore fluid pressure and solid elastic stress is only modelled for a particular scenario with the Biot coefficient for the solid skeleton compressibility equal to 1.0. The fractured rock is loaded in two consecutive phases with far-field stresses applied by a ramping stage to avoid artificial shock. First, an isotropic stress field (i.e. Phase I: r 0 x = r 0 y = 5 MPa) is imposed to consolidate the rock sample under the effective lithostatic stress (Fig. 8a). Second, deviatoric stress conditions are introduced with an increased r 0 x = 10 MPa (i.e. Phase II-A) or 15 MPa (i.e. Phase II-B), and a fixed r 0 y = 5 MPa to consider the evolution of corresponding tectonic regimes (Fig. 8b). Changes in fracture apertures in response to the applied stress conditions are computed by the combined JCM-FEMDEM formulation.
The geomechanical response of the fractured limestone under various in situ stress conditions is modelled in the numerical experiment. Heterogeneity of the fracture-dependent stress field, reactivation with shearing on pre-existing fractures, variation in displacement attributes as well as propagation of new cracks are captured. As shown in Fig. 9a-c, with the increase in stress ratio, pre-existing fractures experience more shearing, accompanied by initiation and propagation of kinked minor cracks if the stress at the tips of pre-existing structures exceeds the critical levels for tensile or shear mode failure. Fracture apertures under the hydrostatic stress condition (Fig. 9d) are quite uniformly distributed with low magnitude. However, in the deviatoric condition, e.g. Phase II-B (Fig. 9f), larger apertures are generated in fractures associated with intensive shearing due to two types of effects: network-scale mesoscopic opening and roughness-scale microscopic dilatancy. Mesoscopic opening occurs in dilational jogs and bends, and between boundaries of relatively rotated blocks as well as within propagated wing cracks, while microscopic dilatancy is accommodated between rough surfaces of shearing fractures.

Comparison with Conventional Joint Modelling Approaches
For simplistic purposes, conventional joint models often assume constant JRC and JCS values according to an average joint length (Kobayashi et al. 2001) and may also assume a constant friction angle based on the residual value (Min et al. 2004). To further demonstrate the potential importance of scale-dependent joint properties and roughness-induced additional frictional resistance, the  Figure 10 shows the simulation results of the three joint models for the fracture network under Phase II-B stress condition (r 0 x = 15 MPa, r 0 y = 5 MPa), loaded through the same two consecutive phases (Fig. 8). Compared to the JCM modelling results, Model-I exhibits a slight difference in the distribution of shear displacement, but overestimates the apertures of some long fractures (larger than the average length) under intensive shearing. This is caused by the exaggeration of dilational displacement of larger fractures by using the joint properties based on a smaller length (see Fig. 5).
Model-II that ignores the asperity effect significantly underestimates the strength of the rock mass and produces

Discussion
Stress-dependent heterogeneity of fracture opening and shear displacement in a naturally fractured rock has been captured by the 2D JCM-FEMDEM model that incorporates both the network-scale mesoscopic effect (e.g. orientations, spacing, junctions, dilational bends, and jogs) and the roughness-scale microscopic effect (e.g. roughness-controlled aperture closure and dilatancy). Integration of the realism of joint constitutive characteristics is considered to give more realistic results compared to conventional approaches that neglect the scale dependency of joint properties and/or the roughness-induced additional frictional resistance as well as its shearing-dependent degradation. The results of the model in Phase II-B (Fig. 9c, f) with a critical far-field stress ratio, i.e. r 0 x /r 0 y = 3, are of particular interest. The system finds equilibrium by activating sliding with local extremes of shear displacement on the favourably orientated joint set 2 as highlighted in Fig. 9c (see further discussion of joint orientation effects in Lei et al. 2014). Locally, the sliding on the two sets has created large apertures in some active fractures as well as their intersections (Figs. 9f, 10), which shows consistency with the field observation from boreholes that critically stressed faults with favourable orientations appear to have larger apertures and higher hydraulic conductivity (Zoback 2007).
The formation of large apertures along displacing and dilating fractures illustrated by the 2D model implies that localised flow might occur in the vertical direction and a higher permeability is expected in the third dimension of the strike-slip faulting system (Sibson 1994), as demonstrated in the work by Sanderson and Zhang (1999) using analytical solutions for vertical flow rate calculation based on the cubic law and the pipe formula. In the 3D geological setting of limestone-shale sequences, aperture variability and even impersistence may exist along the fracture walls normal to the layering, e.g. caused by inhomogeneous filling of calcite minerals.
Geological processes, such as episodes of delamination between layers and fracturing through shales, may make the flow in 3D even more complex, and furthermore, fluid flows are known to be channelized within the bedding planes and fractures, rather than flowing as if between parallel plates. Hence, further work is needed to integrate   has been also combined with the JCM-FEMDEM model to capture the brittle deformation response including local concentrations of critically high tensile or differential stresses, together with realistic fracture opening and shearing behaviour on both pre-existing and newly propagated fractures . Such capability opens the way to modelling 3D flows in geomechanically realistic multi-layer systems with both 'strata bound' and 'non-strata bound' fractures as well as plutonic rock masses.
One limitation of this research is the assumption that deformation of the solid skeleton was determined by the effective stress condition and the direct influence of local internal fluid pressure was not explicitly included. The immersed shell method (Viré et al. 2012(Viré et al. , 2015 and the multiphase flow modelling (Su et al. 2015) that have been recently developed in the research group at Imperial College will be coupled with the 2D and 3D JCM-FEMDEM geomechanical models to capture the non-trivial two-way coupling process involving the transient response of rock solid and pore fluid pressure as well as the dynamic fluidsolid interaction. Some preliminary results have been presented in Obeysekara et al. (2016).
Compared to other discrete element modelling approaches, e.g. the particle-based synthetic rock mass approach (Mas Ivars et al. 2011) and the grain-based Voronoi tessellation method (Damjanaca et al. 2007;Ghazvinian et al. 2014), the FEMDEM model is able to capture the realistic fracturing behaviour of brittle rocks governed by fundamental fracture mechanics principles associated with strength and fracture energy parameters. A detailed review about the FEMDEM method and various other discrete modelling techniques can be found in the paper by Lisjak and Grasselli (2014). The addition of the JCM module to the FEMDEM framework further permits the simulation of the sophisticated shearing behaviour of pre-existing rough fractures based on experimentally derived constitutive laws. Unlike the work conducted with an explicit representation of the fracture roughness profile (Karami and Stead 2008;Bahaaddini et al. 2014;Grasselli 2012, 2015a) that models the underlying process of asperity failure and roughness degradation, the proposed method integrates the well-established empirical joint constitutive laws directly as the criteria for implicit microscale modelling and can be advantageous in applications for large-scale engineering problems. However, these discrete modelling approaches based on an explicit time marching scheme may all suffer from potential dynamic effects in numerical experiments. Although a large damping coefficient can help significantly attenuate the dynamic oscillation and approximate a quasi-static

Conclusions
To conclude, an empirical joint constitutive model that captures the overall behaviour of sheared or compressed individual fractures as observed in laboratory experiments was implemented in the finite-discrete element analysis framework for 2D geomechanical modelling of fractured rocks. The combined JCM-FEMDEM model is able to achieve compatibility for both the fracture and matrix fields with respect to stress and displacement. The numerical model exhibits realistic shear strength and displacement characteristics with the recognition of fracture length influence, which was demonstrated by a comparison with the experimentally derived empirical solutions. 2D plane strain geomechanical modelling based on the combined JCM-FEMDEM formulation was conducted on an outcropbased fracture network. The fracture system response to different stress phases led to a wealth of different local fracture-dominated deformational behaviour. The numerical experiments include the specific local developments of fractures apertures due to the fracture closing, opening, shearing, dilatancy, and propagation. With the increase in stress ratio, significant deformation enhancement occurs in the vicinity of fracture tips, intersections, and bends, where large apertures can be generated. The JCM-FEMDEM model is considered to give more realistic results compared to conventional approaches that neglect the scale dependency of joint properties and/or the asperity effect. The results of this paper have important implications for many rock engineering applications where in situ stress and pore fluid pressure is disturbed including underground construction, geothermal energy, nuclear repository safety, and petroleum recovery.