## 2 Materials and Methods
### Fundamentals of Brittle Failure Numerical Modelling
Although the introduction of a randomly sized polygonal block pattern, such as the Voronoi tessellation (concept of the grain-based model introduced by [21]; for an example, see [22]), together with a constitutive law accounting for shear and tensile strength along the block interface or contacts can help to overcome this drawback, very long computational times are needed when a dense tessellation is used. Another DEM approach is represented by particle flow codes, which can simulate materials at a smaller scale by means of spherical particles [23]. Hybrid, combined finite-discrete element methods (FDEM) can overcome the traditional boundary between continuum and discontinuity methods by capturing emergent discontinuities associated with brittle fracturing processes, progressive damage, and failure processes in rock. The idea is to start with a continuum representation of the domain by finite elements and to progressively move to a discontinuum by the explicit nucleation and growth of fractures as the rock is strained. As stated by [24], FDEM can capture the anisotropic and inhomogeneous effects of natural jointing, being more realistic than methods relying on just continuum or discontinuum. FDEM was first introduced in 1995 by Prof. [PERSON] [25]. Still, it was as of 2008 that the original open-source code Y-2D was further extended by including several essential features for geomechanical applications (Y-Geo code, [20]; Irazu, [26]). The latter will be used in this paper. Ref. [27] illustrates how the FDEM code with a hydro-mechanical coupled formulation in 2D and 3D was parallelised in Irazu using general-purpose Graphics Processing Unit (GPU) computing to reduce execution times and run faster, more extensive, higher-resolution, and multi-physics simulations. In the following, the theoretical fundamentals will be briefly recalled [28] concerning the two-dimensional case. Each intact body is discretised into three-noded, constant-strain triangular finite elements that interact with the neighbouring elements through a four-noded, strain-softening cohesive crack element (Figure 2). The equation of motion can be written as:  
\[\mathrm{M\ddot{x}+C\dot{x}+f_{int}(x)-f_{ext}(x)-f_{c}(x)-f_{h}(x)-f_{r}(x)=0} \tag{1}\]  
where M and C are mass and damping matrices, x is the vector of nodal displacements, \(f_{int}\) is the vector of internal resisting forces, \(f_{ext}\) is the vector of externally applied loads, and \(f_{c}\), \(f_{h}\) and \(f_{r}\) are the vectors of contact forces, fluid-pressure forces, and rock support forces, respectively. This equation is integrated step by step to determine nodal coordinates. Figure 2: (**a**) Example of discretisation of an intact body into finite elements; (**b**) zoom of two neighbouring finite elements together with the interface between them, the so-called cohesive crack element (modified by [28]). Finite elements behave according to linear elasticity, with an isotropic or transversely isotropic constitutive relationship used to derive the stress tensor (and then, equivalent nodal forces) from the strain tensor. The calculated forces are used to solve Equation (1) over the next integration time step. Contacts between all overlapping pairs of elements are detected, and interaction forces are applied: repulsive forces in the normal direction and friction forces in the tangential direction. Crack initiation and growth modelling is based on non-linear fracture mechanics (NLFM). When the material tensile strength starts being exceeded at the tip of a crack (Mode I failure, opening), non-linear behaviour comes into play by giving birth to a fracture process zone (FPZ), characterised by inter-locking and microcracking, albeit still able to transfer loads (Figure 3). A similar behaviour occurs in shear (Mode II failure, sliding). Four-noded cohesive crack elements (Figure 2) approximate the FPZ and its load-carrying capacity. They can yield and break under Mode I, II, or mixed I-II conditions. The fracture nucleation site does not need to be defined a priori.