Coupled Thermo-hydro-mechanical Simulation of Hydraulic Fracturing in Deep Reservoirs Using Finite-Discrete Element Method

Hydraulic fracturing (HF) is one of the most effective stimulation techniques to enhance reservoir permeability. The efficiency of an HF fluid injection depends on the pre-existing discontinuities or sources of heterogeneities and these features need to be considered in a HF operation treatment. Moreover, deep reservoirs are usually located in hot dry rocks (HDR). Hence, thermal conduction through the rock and fluid and advection and convective heat transfer in the fluid can affect the fluid–rock interaction. This study focuses on HF development in deep reservoirs under a high-temperature field. Two separate scenarios are considered: a reservoir containing discrete fracture networks (DFN) and another considering blocks in a matrix as conglomerate reservoirs (there is no relation between the scenarios considered). The study discusses each reservoir separately and simulates their thermo-hydro-mechanical (THM) behaviour using the combined finite-discrete element method (FDEM). First, the capabilities of the FDEM are verified against the existing analytical solutions, and then the FDEM is employed to model HF development. The effects of controlling factors, including flow rate, fluid kinematic viscosity and DFN aperture for jointed reservoirs and flow rate, fluid kinematic viscosity and block strength in conglomerate ones, are studied. The results show that the high fracture density DFNs strongly affect the HF propagation pattern and fluid pressure rise. Moreover, the DFN’s aperture significantly alters the HF treatment behaviour. The controlling factors are observed to influence the HF pattern strongly, and a successful HF treatment requires careful consideration of all the factors. In the conglomerate reservoirs, the strength of the blocks strongly dominates the HF mechanism, in which soft blocks break and allow for uniform fluid pressure distribution and longer HFs, while hard blocks stop fluid from flowing over longer distances accumulating high fluid pressure around the injection. This mechanism excessively breaks the matrix and reduces HF efficiency. Crack branching frequently occurs in conglomerate reservoirs due to thermal exchange between the blocks, matrix, and fluid. Coupled THM FDEM is utilized to study hydraulic fracturing in HDR. Reservoirs containing blocks in matrix or DFN are considered. The effects of controlling factors are evaluated. Recommendations to improve the efficiency of a fluid injection are suggested. Coupled THM FDEM is utilized to study hydraulic fracturing in HDR. Reservoirs containing blocks in matrix or DFN are considered. The effects of controlling factors are evaluated. Recommendations to improve the efficiency of a fluid injection are suggested.


Introduction
With high energy consumption and rapid climate change, the utilisation of alternative energy resources, including shale gas, tight gas and oil, and geothermal energy in deep reservoirs, has attracted more attention and is becoming more and more vital. Since these energy resources are in deep reservoirs, they are under high in situ stresses and high temperatures, making their behaviour complex for energy extraction purposes. On the other hand, these types of reservoirs usually have extremely low porosity and permeability (usually below 0.1 mD) (Lu et al. 2015;Hou et al. 2015). Therefore, it is a challenging task to deal with such complex reservoirs.
It is necessary to enhance the permeability of the abovementioned complex reservoirs to allow for efficient energy extraction. Hydraulic fracturing, commonly known as fracking, injects water, sand, and chemicals into a well to break up underground bedrock and free up oil or gas reserves. Hydraulic fracturing is a well-stimulation technique commonly used in low-permeability rocks like tight sandstone, shale, and some coal beds to increase oil or gas flow to a well from petroleum-bearing rock formations. Hydraulic fracturing began as an experiment in 1949 (Clark 1949), and the first commercially successful application followed in 1950 (Suchy and Newell 2011). Since 1950, it has been used widely to extract gas, hot water, and storage and extraction of thermal energy (Sherratt et al. 2021;Zhao et al. 2020;Ren et al. 2020;Yu et al. 2020;Hu and Ghasemi 2020;Umar et al. 2021;Valley and Evans 2019). This method has also attracted attention in mining and civil engineering with different applications for relieving high stress, extending excavation size, enhancing cavability of hard coals, avoiding rock burst, in situ rock stress measurement and pre-conditioning of rock in block caving mining (Amadei and Stephansson 1997;Bai et al. 2017;Chen et al. 2020a, b;He et al. 2016;Li et al. 2021a, b;Liu et al. 2020a, b;Zhai et al. 2012).
Hydraulic fracturing is a complex operation in which its successful application depends upon many factors, more importantly on the reservoir's geological conditions. In homogenous and isotropic rock without any discontinuity, hydraulic fractures (HFs) tend to propagate in the direction with the least resistance, which is normal to the minimum principal stress direction (Cheng et al. 2018;Zhong et al. 2021;Sharafisafa et al. 2023). However, natural rock masses usually contain a particular source of anisotropy which causes complex deformation and failure behaviour. Amongst the anisotropic sources, discontinuities such as flaws (Sharafisafa et al. 2018), cracks, joint sets, bedding, blocks and fault are very common and are always confronted in civil and mining projects (Aliabadian et al. 2021;Sharafisafa et al. 2020a, b;Sharafisafa and Shen 2020). These geological features not only affect the overall deformation behaviour of the rock masses but also degrade their fracture characteristics and impose strong influences on the crack propagation mechanism (Li et al. 2021a, b;Park and Bobet 2010). These planes of weaknesses also significantly affect the HF propagation mechanism and can lead to a complex hydraulic fracture geometry (Kamali and Ghasemi 2018;Fu et al. 2022). Given the interaction of a HF with the existing discontinuities, fracture branching, offset crossing and non-planar fracture growth are widely observed in the field and laboratory experiments (Dehghan et al. 2015;Gracia et al. 2013;Hu and Ghasemi 2021;Lamont and Jessen 1963). The interaction between the natural discontinuities and the generated HFs have already been well studied using various methods, including analytical solutions (Sarmadivaleh and Rasouli 2014;Zhao et al. 2019), experimental (Hu and Ghasemi 2021;Dehghan 2020;Ham and Kwon 2020;Ma et al. 2019), artificial neural network (Silveria et al. 2022) and numerical modelling Yan et al. 2022a, b;Sun et al. 2022;Dong et al. 2021;Liu et al. 2020a, b;Xu et al. 2019;Zhao et al. 2021;Zhou et al. 2020) and many more.
Besides the geological conditions of reservoirs, temperature of the reservoirs also plays a significant role in the development and efficiency of hydraulic fracturing. Most of the existing studies on the HF development have been conducted on low-temperature rocks at room temperature, while in EGS usually the rock temperature exceeds 200 °C. Temperature increase in rocks causes changes in physical and mechanical characteristics, reducing their tensile and compressive strengths and elastic modulus and degrading their properties Quansheng 2000;Zhang et al. 2011). Moreover, the large temperature difference between the injecting fluid and the host rock leads to a complicated thermo-hydro-mechanical coupling between the rock grain and the fluid in a hydraulic fracturing operation. This further causes intense thermal shock to the rock and its components, which is then followed by a thermal failure (Chaki et al. 2008;Fu et al.2007). The presence of anisotropy sources such as joints, flaws and blocks make the interaction between the injecting fluid and host rock more complex. Thermal conduction between the injecting fluid within fractures and the host rock has been the subject of some studies (Wang and Papamichos 1999;Tran et al. 2013;Enayatpour and Patzek 2013). It is usually assumed that the rock temperature at the surface of a HF is constant and equal to the temperature of the fluid (Tran et al. 2013). However, such an assumption is not valid and does not meet the conservation of energy, since it does not consider the fact that the thermal exchange between the injecting fluid and the surrounding rock leads gradually to the thermal equilibrium with the rock matrix. As a result, unrealistic modelling is carried out, which cannot accurately predict the real behaviour of rocks under hightemperature conditions (Salimzadeh et al. 2016).
On the other hand, several methods, including laboratory experiments, field experiments, analytical solutions, and numerical simulations have been used to study the hydraulic fracturing mechanism and efficiency. Amongst the utilised methodologies to study the HF propagation mechanism, numerical simulations have been extensively used due to their capabilities in modelling complex behaviours. Several numerical techniques can be utilized to study the behaviour of HF under various environments and reservoir conditions. For instance, the finite element method (FEM) is the most common and developed methodology in the field. However, this method encounters a fundamental problem when dealing with a fracture, and complex remeshing is usually required, reducing the simulation efficiency (Belytschko and Black 1999). To comply the FEM with a fracture, a cohesive zone model (CZM) was implemented in the FEM to avoid remeshing; however, the fracture development path needs to be predefined which restricts the applicability of the method (Chen et al. 2009). The CZM was then extended to an unknown fracture trajectory to overcome this problem by modelling cohesive interface elements between all the finite elements (Wang 2019). The extended finite element method (XFEM) has shown strong capabilities in modelling fractures under the assumptions of a FEM analysis without the need for remeshing, which increases the efficiency of this method compared to the conventional FEM. This method has also been used to model the HF development in rocks (Dogn et al. 2021;Jafar et al. 2021;Zhou et al. 2021;Sharafisafa and Nazem 2014). However, modelling complex fracture networks and three-dimensional non-planar fractures by XFEM is still challenging and of low accuracy (Lacampion et al. 2018). Another method is the displacement discontinuity method (DDM) which discretises the region by fracture surface to enable fracture propagation ). However, the DDM faces fundamental issues when heterogeneities and inelastic deformation dominate the behaviour of a rock mass. The discrete element method (DEM) is a robust and powerful simulation tool developed by Potyondy and Cundall (2004), which is capable of modelling the evolution of fractures by using bonded particles. Moreover, DEM can model fluid flow in fractures and the complex interaction between HFs and fluid flow in a discontinuum rock mass (Ghader et al. 2018;Kong et al. 2021). However, the computation efficiency of the DEM decreases with the increase in the number of particles. Thus, the application of the DEM in large-scale models is restricted. There are some other numerical tools; each has its advantages and disadvantages when dealing with complex fracture network propagation in rock masses. Therefore, a suitable numerical method should be carefully selected in accordance with the rock mass characteristic and features.
The abovementioned numerical methods all suffer from some disadvantages when dealing with fracture initiation and propagation in geologically complex rocks. Hence, the combined finite-discrete element method (FDEM) has been proposed to overcome the shortages of the existing numerical methods. The FDEM proposed by Munjiza (2004) combines the advantages of the FEM for solving problems of solid deformation and the DEM for simulating contacts between deforming blocks. Furthermore, fluid flow can also be considered by generating flow channels at the interfaces between adjoining finite elements and these flow channels connect virtual cavities corresponding to the nodes of the original finite elements. At the same time, FDEM models the fracture development via the breakage of embedded joint elements. The FDEM has recently been widely used to model fracture propagation in rock masses (Yang et al. 2022;Chen et al. 2020a, b;Deng et al. 2021;Liu et al. 2021;Wang et al. 2020). The FDEM has also shown strong capabilities in modelling fluid flow in rock matrix by considering hydromechanical or thermo-hydro-mechanical assumptions (Sun et al. 2019;Yan et al. 2022a, b).
The deep jointed or conglomerate reservoirs are increasingly getting noticed for the purpose of energy extraction. Despite their widespread existence, their failure behaviour and fracture propagation mechanism in a hydraulic fracturing treatment under the fluid injection process are largely ignored. Moreover, the complexity of such reservoirs becomes even more pronounced when a high-temperature field also exists. The complicated interactions between the DFNs or blocks in a matrix, in situ stresses, high-temperature field, and cold injected water make the issue very complex. Furthermore, a few studies have made some efforts to explore the HF development behaviour considering a change in controlling factors, including flow rate, fluid kinematic viscosity, and in situ stresses. Therefore, a thorough analysis is required to investigate the HF development under such conditions. Since the FDEM has proven its robust capabilities in modelling complex grounds considering complex conditions, it is used in this study to systematically investigate the HF propagation pattern in deep hot jointed or conglomerate reservoirs when cold water is injected. The main objective is to evaluate HF propagation and its interaction with random DFNs or blocks under high temperatures and in situ stresses. In this study, the combined-finite-discrete element method (FDEM) is utilised to study the complex interaction between the pre-existing natural sources of heterogeneity and discontinuity and the HF under high temperatures. First, the mathematical background and theory of the thermal logic in FDEM and thermo-hydro-mechanical coupling are briefly explained. Then the validity of the FDEM is verified by comparing the FDEM results and the analytical solutions. Next, the interaction of the HF and some sources of heterogeneities such as low and high distribution density DFNs and blocks of different strength in a matrix are investigated. Finally, the results are analysed to make several recommendations for a successful HF treatment in deep hot rocks.

Coupled Thermo-Hydro-Mechanical (THM) Model in FDEM
The combined Finite-Discrete Element Method (FDEM), a numerical technique primarily proposed by Munjiza et al. (1995), allows the modelling of bodies having various interactions. A modelling stage can start with a single intact region or a combination of discrete intact domains. By progressing the modelling, these domains (bodies) can undergo elastic deformation, rotate, interact and translate, and experience fracturing once some fracture criterion are satisfied, thereby generating new discrete bodies or regions. The newly formed bodies can then undergo more displacement, deformation, interaction, and fracture. The method utilises the Finite Element Method (FEM) to investigate the solid deformation and assess the failure criterion for fracturing and the Discrete Element Method (DEM) to detect newly formed contacts and elaborate the translation, rotation, and interaction of discrete bodies. This method takes advantage of the continuum-based FEM and the discontinuum-based DEM. The FDEM is not only capable of modelling the damage of brittle materials prior to fracture but also is able to simulate the complex interaction of the fractures. During every hydraulic time step, the following main calculations are sequentially calculated; determination of volume and saturation of the cavities, as well as a hydraulic aperture of the channels; application of Darcy's law and parallel-plate flow model at each channel; application of the continuity equation at each cavity; and computation of the fluid pressure field. The fundamental of the FDEM can be found in Munjiza et al. (1995) and is not repeated here. However, the fundamentals of the coupled THM in FDEM is given here briefly. THM modelling uses the same modelling framework adopted for the mechanical and hydraulic calculations while adding new thermal-related variables (i.e., thermal energy and temperature) to the nodes and cavities and new functions capturing several thermal transfer processes. This paper employs the 2D FDEM Irazu V. 5.1.0 for the HF development (Geomechanica 2022). Irazu software has strong capabilities for a full thermo-hydro-mechanical coupling, which enables us to precisely model the HF development under complex subsurface conditions (Lisjak et al. 2017). The thermal solver of Irazu accounts for the heat transfer processes as outlined below (reproduced from Geomechanica 2022).
For thermal conduction in solids, thermal energy is transferred between the nodes of each triangular finite element. In terms of thermal conduction in fluids, thermal energy is transferred between the fluid-filled virtual cavities connected by flow channels. For thermal advection, thermal energy is transferred between the fluid-filled virtual cavities via the bulk movement of the fluid inside the flow channels. For convective thermal transfer, thermal energy is transferred between the fluid-filled cavities and the nodes of the triangular finite elements. In terms of contact thermal transfer, thermal energy is transferred across broken crack elements between the nodes of the contacting pairs of finite element edges. The heat transfer processes listed above are illustrated in the context of the computational framework of Irazu using the schematics of Fig. 1a.
The formulation for the thermal transfer at the contact between finite elements is based on the topology of (broken) crack elements, which are pre-inserted and can be activated when the fracture criterion is met. Therefore, it is strictly valid only for the cases of small slip along fractures. In other words, only edge-to-edge contacts between initially jointed finite elements are considered (Fig. 2a). The thermohydro-mechanical (THM) coupling is based on a two-way, explicit coupling approach obtained by alternating between the three solvers while marching forward in time (Fig. 2b). In a typical simulation, the hydraulic step is followed by the thermal step, which is followed by enough mechanical steps to guarantee quasi-static conditions. Hydraulic and thermal calculations are typically synchronised (i.e., the smallest of the two critical time step sizes is adopted for both solvers), while the mechanical solver uses its own (much smaller) time step size.
The thermo-mechanical (TM) coupling is based on a two-way explicit coupling approach obtained by alternating between the mechanical and the thermal solvers while marching forward in time for solid (Fig. 2c). In the thermalto-mechanical (T → M) coupling, the solid temperature field causes the expansion or contraction of the finite element mesh and thermally induced stresses may naturally emerge in the model (depending on the boundary conditions and/or inhomogeneous thermal deformation).
The initiation and propagation of explicit fractures are simulated in Irazu using concepts of non-linear fracture mechanics (NLFM) (Fig. 2d). According to this theory, as the tensile strength of a material is exceeded at the tip of an opening (Mode I) crack, a zone characterized by non-linear behaviour, called the fracture process zone (FPZ), begins to form. In brittle geomaterials, the FPZ is characterized by interlocking and micro-cracking and, albeit damaged, is capable of transferring load across the fracture walls. Within Irazu, an analogous FPZ behaviour is also assumed to exist ahead of crack tips subjected to Mode II loading.
Depending on the local stress and relative crack wall displacements, the crack elements may yield and break under Mode I (i.e., opening mode), Mode II (i.e., sliding mode) or mixed-Mode I-II conditions. According to a cohesive model, Mode I fracture initiation and propagation occurs when the opening of crack element, o, reaches a critical value, o p , corresponding to the intrinsic tensile strength of the element, f t (Figure a below). As a crack element is opened beyond o p , the normal stress, σ, is assumed to gradually decrease until o exceeds the residual opening value, or, at which point a traction-free surface is formed.
According to a slip-weakening model in the software, Mode II fracture initiation and propagation occurs when the tangential slip of a crack element, s, reaches a critical value s p , corresponding to the intrinsic shear strength of the element. As a crack element sustains further slip beyond s p , the tangential stress, τ, is assumed to gradually decrease until s exceeds the residual slippage, s r . At this point, a physical discontinuity is formed and τ becomes equal zero. Further frictional resistance along broken crack elements is calculated by the contact interaction model. Instead of pure Mode I or pure Mode II displacements, crack elements often experience a combination of the two (mixed-Mode I-II displacement). With mixed-Mode I-II displacement, the o and s of a crack element may be less than o r and s r , respectively, but the resultant crack displacement can be quite large.
In the mechanical-to-thermal (M → T) coupling, the failure of one or more crack elements (cohesive element) around a node/cavity location triggers the following effects ( Fig. 3a): (1). the number of flow channels contributing to the cavity increases, (2). the original group of nodes is split (based on the new topology of the broken crack elements) and the nodes are re-grouped to form new groups, (3). the number of finite element-flow channel pairs across which convective thermal transfer takes place increases, and (4). the number of finite element pairs across which contact thermal transfer takes place increases. The thermohydraulic (TH) coupling is based on a two-way explicit coupling approach, obtained by alternating between the hydraulic and the thermal solvers while marching forward in time (Fig. 3b).
During each thermal time step, the following main operations are sequentially performed: 1: computation of heat transport at each finite element and flow channel; 2: application of continuity equation at each node and cavity; and 3: update of solid and fluid temperatures. The mathematical formulation at the basis of the three steps is described in the following subsections.

Thermal Conduction in Solids
The transport of thermal energy is accounted for by applying Fourier's law of heat conduction, a constitutive law relating the thermal flux vector, q, to the solid temperature gradient, ∇T s ; for a homogeneous and isotropic solid: where s is the thermal conductivity of the solid. This equation is numerically integrated over each triangular element surface ( Fig. 4a) to compute the thermal flow rate, Q, in each node assuming a linear variation of temperature within the triangle and applying Gauss' divergence theorem. For instance, using a local coordinate system (x,y), the expression for the thermal flow rate at node 0 is given by: where T A , T B , T C are the temperatures at the three midedge points defined by coordinates (x A , y A ), (x B , y B ), and (x C , y C ), respectively, and A is the surface area of the triangle.

Thermal Conduction in Fluids
Thermal conduction inside fluid-filled fractures is governed by an equation analogous to the one used for the solid conduction but with the thermal conductivity of the solid, s , replaced by the thermal conductivity of the fluid, f . The transport equation is applied to each flow channel (Fig. 4b) to compute the variation of thermal energy in the virtual cavities: (2) where a and L are the hydraulic aperture and length of the flow channel, respectively; ad T f,0 and T f,1 are the fluid temperatures in the two cavities.

Thermal Advection
Movement of thermal energy via the bulk flow of the carrying fluid inside fractures is described by the advection equation: where v is the fluid velocity vector. Using again the divergence theorem, this equation is integrated for each flow channel ( Fig. 4c) to compute the thermal flow rate into the downstream cavity: where T u and T d are the fluid temperatures in the upstream and downstream cavities, respectively; f and C f are the fluid density and specific heat capacity; Q f is the fluid flow rate vector; and n is the flow channel axis vector.

Convective Thermal Transfer
The heat transfer at the solid-fluid interface is assumed to be governed by Newton's law of cooling: where h conv is the convective heat transfer coefficient and ΔT is the differential temperature between solid and fluid. Convective heat transfer is assumed to occur between fluid-filled flow channels and solid finite elements (Fig. 4d). Equation 6 is integrated over the length of the flow channel walls, L, in order to update the thermal flow rate of the connected nodes and cavities. For example, the thermal flow rate at node 0 is given by: where T s,0 and T f,0 are the solid and fluid temperatures at node 0 and cavity 0, respectively.

Contact Thermal Transfer
Upon breakage of a crack element, the thermal transfer between contacting finite elements ( Fig. 4e) is assumed to be governed by an equation similar to Eq. 6 with the coefficient h conv replaced by the contact conductance, h cont . Such an equation is integrated over the contact length of the two finite element edges, L c , to determine the thermal flow rate, Q conv,0 , at each node.

Energy Balance and Temperature Update
The variation of thermal energy is computed by integrating each flow rate with an explicit temporal integration scheme: where ΔT t is the thermal time step size. The conservation of energy is then applied to update the thermal energy stored at each node: where ΔQ t cond,s , ΔQ t conv and ΔQ t cont are the node net thermal energy variations due to solid conduction, convective transfer, and contact transfer, respectively. Similarly, the thermal energy stored in each cavity is updated as: where ΔQ t cond,f , ΔQ t conv , and ΔQ t cont are the cavity net thermal energy variations due to fluid conduction, advection and convective transfer, respectively. Using the updated thermal energy values at the end of the time step, the node (solid) and cavity (fluid) temperatures are updated as: where m s and m f are the (solid) nodal and (fluid) cavity masses and C s and C f are the solid and fluid specific heat capacities.

Validation Examples
It is essential to ensure that the employed method can accurately simulate all the THM processes. Therefore, this section validates the FDEM method by benchmarking the results of the HM, TM and TH couplings against the corresponding analytical solutions.

HM Coupling Validation
This numerical example considers the elastic response of a 21.4 m-long linear discontinuity of zero initial width embedded in a linear elastic, homogeneous, isotropic, impermeable rock mass subjected to a constant, uniform internal fluid pressure, P = 5 MPa (Fig. 5a). This model aims to verify the correctness of the coupling function used to apply an inter-element fluid pressure in the mechanical calculations as well as the overall elastic response of the solid material. The domain is meshed with 59,810 triangular finite elements with a nominal size of 0.15 m along the discontinuity, graded to 3 m along the exterior domain boundary. A fully coupled mechanical-hydraulic simulation is carried out. The following values are assumed for the rock's elastic constants: Young's modulus, E = 40 GPa, and Poisson's ratio, υ = 0.22. Artificially high strength parameter values are used to enforce a linear elastic behaviour. To minimise the effect of the crack element compliance on the overall model stiffness, the penalty values, p n , p t and p f , are set equal to 10 × the Young's modulus value, respectively. The mechanical aperture along the discontinuity can be analytically calculated as (Parker 1981): where L is the length of the discontinuity (21.4 m) and x is the horizontal axis with origin at the centre of the discontinuity. As shown in Fig. 5a, b good agreement was achieved between the FDEM numerical and analytical results.

TM Coupling Validation
This section is devoted to the verification of the FDEM for modelling the heat transfer problems in rocks. Two examples are considered here. The first verification example considers a hollow disc of initial temperature T t0 (= 20 °C) subjected to internal heating by a prescribed inner temperature T i (= 100 °C) at the inner circle and a constant temperature of T o (= 20 °C) at the outer circle. Both the inner and outer circles are free boundaries. A two-way-coupled T ↔ M simulation is carried out with the geometry and boundary conditions shown in Fig. 6a and thermo-mechanical properties listed in Table 1. The steady-state temperature distribution as a function of the radial distance r is given by (Eslami et al. 2013): where r i and r o are the inner and outer radii of the hollow cylinder, and T i and T o are the temperatures at the inner and outer circles. The radial and hoop stresses are given by:  Figure 6b shows the temperature distribution in the disc. The comparison between the analytical and numerical results are reported in Fig. 6c, which shows a good agreement between the two solutions, thus verifying the thermal conduction in solid and the TM coupling process.

Fluid-Solid Heat Exchange
This verification problem deals with the transient heat exchange between the fluid-filled fracture and surrounding matrix. It considers a single horizontal fracture of applying P A = 24 kPa (inlet) and P B = 0 (outlet) to the left and right ends of the fracture, respectively, resulting in a flow velocity 0.01 m/s. Then, a constant fluid temperature equal to T A = 20 °C is applied to the fluid at point A, and the advection in the fracture, conduction in the matrix, and fluid-solid heat exchange across the fracture walls were simulated. Figure 7b illustrates the model meshes. The rock and fluid properties are listed in Table 2.
The analytical solution to this problem is given by Lauwerier (1955), in which the fluid temperature distribution as a function of time and distance from the injection point is given by: while the solid temperature distribution as a function of time and position is given by: where T i s is the initial temperature of the system, s is the thermal conductivity of the matrix rocks, s, f , C s andC f are the density and heat capacities of water and matrix rock respectively, v f is fracture aperture, and erfc is the complementary error function. Figure 7c shows the simulated temperature distribution in the solid after 46 days of fluid injection. As highlighted in the plots of Fig. 8, a good match is obtained between the numerical and analytical results of the fluid temperature along the fracture as well as the matrix temperature along a horizontal monitoring line. It is worth noting that the (17)  analytical solution assumes the fluid and matrix temperature to be equal along the fracture walls. In the numerical model, however, this condition is not explicitly enforced and a convective heat transfer model is adopted instead.
A sensitivity analysis to the convective heat transfer coefficient shows that the numerical solution converges with the analytical solution with increasing values of h conv and that for h conv greater than 100 the difference becomes negligible.

Thermal Conduction and Advection Through a Fracture
This example verifies the thermal transport through a saturated fracture of uniform aperture due to fluid advection and conduction. The geometry and boundary conditions of the one-way TH coupled simulation problem are depicted in Fig. 9a. The length and width of the model are 1 m and 0.1 m, respectively. The steady-state flow conditions are first obtained by applying a pressure gradient to the fracture. Then, constant fluid temperatures equal to 50 °C and 20 °C are applied to the two fracture ends, respectively, while the rock block is kept in adiabatic conditions. Fluid density, dynamic viscosity, specific heat capacity and thermal conductivity are set to 1000 kg/m 3 , 0.001 Pa.s, 4186 J/kg/ o C, and 0.5918 W/m/ o C, respectively. The cases of advection and conduction acting in the same (Case A) and opposite (Case B) directions are considered by changing the direction of the hydraulic pressure gradient. The simulated steady-state flow velocity through the 1 m-long, 0.001 mm-wide fracture is equal to 4.2e-7 m/s. The steady-state temperature distribution for the verification example described is given by (Ogata and Banks 1961): where T A and T B are the temperatures at the fracture end points, L is the length of the fracture, and Pe is the Peclet number defined as: where ρ f , C f and q f are the density, specific heat capacity and thermal conductivity of the fluid, respectively; λ f is the specific discharge of fluid through the fracture. As depicted in Fig. 9a, b, perfect match was obtained between the analytical and numerical fluid temperature profiles along the fracture. Results also show that thermal advection acts by shifting the temperature distribution in the direction of the fluid flow.

Simulation Methodology
In this section, the validated FDEM is employed to simulate the process of HF initiation and propagation and its interaction with pre-existing sources of anisotropy and heterogeneity in deep rocks. All the simulations will include the coupled thermo-hydro-mechanical (THM) analysis, enabling accurate modelling of the injection fluid into rocks. All the 2D models have the dimensions of 5 m × 5 m with two meshing regions being used. A mesh refinement zone with mesh size of 0.1 m is introduced within the primary model with the size of 4 m × 4 m and the outer surface is meshed gradually to 1 m with a transition speed of 2. By this mesh arrangement, a fine mesh size inside the mesh refinement zone is introduced around the injection point and the pre-existing features to increase the analysis precision. To simplify the analysis, the fluid injection borehole is modelled by an injection point in a node of a mesh at the center  Table 3 (adopted partially from Lisjak et al. 2014) and the parameters assigned for the microcracks in DFN are given in Table 4. The models for both DFN and blocky ones are shown in Fig. 1a. Note that the rock is assumed impermeable with no leak-off.
In the models with DFNs, the fractures are considered closed, and the total volume of the DFNs is much smaller than the rock mass volume. As a result, HFs will propagate to a short distance; thus, the interaction between the HFs and the DFNs cannot be investigated precisely. Moreover, usually in deep reservoirs containing DFNs, the permeability of the rock itself is very low, and in most studies (Li et al. 2022;Hu et al. 2022;Liu et al. 2022;Pu et al. 2021), these reservoirs are considered impenetrable or with extremely low permeability. Furthermore, water-rock interactions at elevated temperatures may significantly decrease its permeability due to the healing of existing fractures. Thus, in these models, the rock mass is considered impermeable, allowing the fluid to flow only through the DFNs and newly formed HFs for a better investigation of the interaction of the fluid with such discontinuities under the high-temperature field. However, usually in conglomerate reservoirs, the strength of the matrix is low compared to the blocks, while the permeability of the matrix is higher than the blocks. In studies on hydraulic fracturing in conglomerate reservoirs, the matrix is considered permeable, while the blocks are usually considered impermeable or have extremely low permeability. This assumption makes realistic modelling of conglomerate reservoirs so that a fluid can flow through the matrix and newly formed HFs and interact with the blocks. Therefore, the interaction of fluid and blocks can be simulated accurately.
The initial conditions are given as follows: the initial temperature of the rock is 200 °C, the flow rate at the injection point is 1 kg/s, the initial injection fluid temperature is 20 °C, and the in situ stresses are initially set as isotropic 6 MPa in both directions. Two models are simulated and analysed: a model with sets of DFNs and a model with blocks in a matrix as shown in Fig. 10a. For each model, its specific conditions and modelling procedure will be explained in its related section.

Model with DFNs
It is revealed that the HFs interact with pre-existing natural discontinuities, such as in the stimulated fracture pattern in dilation, crossing or activation. Pre-existing natural discontinuities are an important source for multistranded and non-planar HFs since some natural discontinuities might be reactivated or crossed by the propagation of the HFs, and thus form a complex fracture network instead of a straight fracture. The interactions between the natural discontinuities and the propagating HFs have been well studied using various methods, including analytical solutions (Cheng et al. 2014), neural networks (Silveria et al. 2022), experimental   (Fatahi et al. 2017), and numerical (Cheng et al. 2019).
All the analysis methods have confirmed that five types of interactions occur between the natural discontinuities and HFs: crossing, arresting without crossing, dilation, activation and offset crossing (Fig. 10b). Moreover, the mixed types of these interactions have been observed. The type of interaction depends on many factors, including shear strength of natural discontinuities, in situ stresses, injection flow rate, intersection angle, the spatial distribution of the natural fractures and others. Interested readers are referred to Cheng et al. (2015) for a detailed explanation of the types of interaction.
To numerically reproduce and verify the mentioned interaction types, a model with a set of DFNs is built, and the assigned values of the DFNs are given in Table 3. The distribution of the cracks follows a Gaussian distribution model with a standard deviation of 10°. DFNs spacing is not considered and only angle varies to make a random DFN. The mean length of the cracks is 5 m and the mean fracture number density values is 0.2 1/m 2 with a standard deviation of 0.05 1/m 2 . An injection point is placed at the centre of the model. A fluid with 20 °C temperature is injected into the model at the injection point. Figure 10(c1-c7) illustrates the identified interaction types in the model. A new interaction type is also observed in Fig. 10(c8) which is here called dilation plus activation. In this interaction a HF arrives at one side of a DFN with fluid infiltrating within the joint and the dilation causes the initiation of a new HF from the other side of the joint. Figure 10(c1-c8) shows that all the interaction types occurred during the injection and were captured successfully by the FDEM. Note that the dilation mode is identified by the opening of DFNs and the high fluid pressure in a joint, as illustrated by the colour gradient scale.

Influence of the Flow Rate
The rate of change in the fluid injection can significantly affect the HF path. This section assesses this parameter's influence on the HF's development pattern. Two in situ stress states are considered while changing the flow rate. First, an isotropic stress state of σ H = σ h = 6 MPa is established before the fluid injection. In the second scenario, an anisotropic stress state of σ H = 15 MPa and σ h = 5 MPa is considered. The flow rates are increased from 1 to 2 and 3 kg/s with a constant fluid kinematic viscosity of 1e-6 m 2 /s. DFNs with Fig. 11 Effect of flow rate on the HF propagation pattern under a isotropic and b anisotropic in situ stress states after 2 s of fluid injection for models with low fracture density DFNs low fracture density have no overlapping or intersection between the DFNs, while high fracture density DFNs refere to DFNs with frequent intersection. Figure 11a illustrates the developed HFs in an isotropic confining stress state of 6 MPa under various flow rates for the models with low frequency DFNs. Under the low flow rate of 1 kg/s, long HFs propagate from the injection point with similar development lengths. However, increasing the flow rate to higher values leads to the propagation of numerous short and long HFs from the injection point. In addition to the increased number of HFs under the higher flow rates, the HFs in the models with high flow rates exhibit more crossing types, while the activation and dilation + activation types dominate the interaction types under low flow rates. Also, the higher flow rates cause higher fluid pressure around the injection point, which in turn causes many short HFs and excessive rock mass breakage. Figure 11b shows the developed patterns of the HFs under anisotropic confining stresses for the models with low frequency DFNs. Under the low flow rate of 1 kg/s, HFs tend to develop in the direction of the maximum principal stress, while increasing the flow rate to higher values leads to the propagation of the HFs in all directions (but still preferably in the direction of the maximum principal stress). This trend indicates that the injection of fluids with higher flow rates minimises the influence of the confining in situ stresses. The dominant interaction type under low flow rates is observed to be crossing which is more under the influence of the confining stresses than the flow rate. However, the higher flow rates impose higher fluid pressure on the joints surfaces, which not only dilates them but also facilitates the HF propagation speed under higher stress concentration in the joints tips and thus causes more dilation + activation as the interaction type.
In order to better investigate the combined influence of the flow rate and DFN on the HF development pattern, models with high fracture density DFNs and the same dimensions as those in Fig. 11 are built. There are two sets of DFNs with the following mean and standard deviation values for dip (angle), length (m) and fracture density (m −1 ), respectively; set one 30,5; 5,1; 1,0.2 and set two 120,5; 5,1; 1,0.2. All other properties including materials and fluid properties remain unchanged. Figure 12 shows the results of the HF development pattern for these models under the flow rates of 1, 2 and 3 kg/s. For the isotopic confining stress conditions under the 1 kg/s flow rate as shown in Fig. 12a, the injected fluid creates Fig. 12 Influnece of the flow rate on the HF development pattern in models with high fracture density DFNs under a isotropic and b anisotropic confining stress state after 2 s of fluid injection some new HFs around the injection point, then fluid flows from the HFs into the DFNs, and finally, a uniform fluid pressure distribution occurs. An increase in the fluid flow rate to 2 kg/s shows more HFs' generation around the injection point and exhibits nonuniform fluid pressure distribution. This behaviour is also observed in the model with a 3 kg/s flow rate but with excessive nonuniformity in the fluid pressure distribution and the generation of more HFs around the injection point. It is observed that under the higher flow rate, the excessive fragmentation around the injection point suspends the fluid flow and leads to the high fluid pressure accumulation and shorter HFs, and thus decreased efficiency. In the anisotropic confining in situ stress regime as shown in Fig. 12b, under the low 1 kg/s flow rate, new HFs develop in the direction of the maximum principal stress, fluid then flows in a small perimeter around the injection point and a few DFNs are activated. The generation of new HFs, fluid flow distance and DFN activation for fluid flow are observed to increase once the flow rate rises to 2 kg/s. However, a high flow rate of 3 kg/s shows excessive fracturing around the injection point, leading to more DFN activation and pushing the fluid to flow further into the surrounding DFNs. Moreover, a relatively large area in the direction of the maximum principal stress is observed to have fluid pressure. Therefore, the high fracture density DFNs cause more HF arrestation, but further DFN activation and thus large fluid extension.

Effect of Fluid Kinematic Viscosity
The HF formation behaviour and pattern is largely influenced by the fracturing kinematic viscosity (μ). It is known that high fluid kinematic viscosity leads to the development of larger volume and wider fractures in the reservoirs (Nagaso et al. 2015). In a numerical study by Sarris and Papanastasiou (2015), the influence of fluid viscosity and the injection rate on a fluid-driven fracture in cohesive poroelastic and poroelastoplastic weak formations was investigated. Their result showed that under the injection of fluid with a small value, the action of the viscous fluid flow does not impose any critical effect on the fracture profiles. The fracture is wider when it is driven with highly viscous fluids. The formed aperture with the low-viscosity fluids is narrower. Also, Wrobel et al. (2021) showed that the rheological properties of fluids considerably influence the process of hydraulic fracture not only by the viscosity values but also by the range of fluid shear rates over which variation of viscosity occurs. The influence of the temperature on the fluid viscosity is not considered in this study Fig. 13 Influnece of the fluid kinematic viscosity on the HF development pattern in a low fracture density DFN model under a isotrpoic and b anistropic confining stresses after 1 s of fluid injection and a constant kinematic fluid viscosity is considered in all simulations. Thermal-to-hydraulic (T → H) processes, such as the dependency of fluid physical properties (e.g., density, viscosity) on temperature, are not accounted for in the current version of the software. Moreover, the degree of the kinematic viscosity of an injected fluid into the reservoirs alters as the temperature and the pressure change. Therefore, it is of great importance to assess this factor on the behaviour of the HF propagation. In this regard, this section analyses the effect of the fluid kinematic viscosity in the models with the two mentioned DFN sets. The isotropic and anisotropic confining in situ stresses are considered. The kinematic viscosity varies from 1 to 2, 3 and 4 mm 2 /s. The flow rate is kept at 1 kg/s. Figure 13 illustrates the development of HFs under the kinematic viscosity of 2, 3 and 4 mm 2 /s for the low fracture density DFNs under both isotropic and anisotropic confining stresses. Under the isotropic confining stress condition and the 2 mm 2 /s fluid kinematic viscosity as shown in Fig. 13a, individual HFs are observed to propagate in all directions and crack development density around the injection point is low. Increasing the fluid kinematic viscosity to 3 mm 2 /s leads to severe breakage around the injection point and then the development of a large number of HFs in all directions. The case with μ = 4 mm 2 /s exhibits severe breakage around the injection point and the development of many Hfs from the breakage zone. Moreover, the length of the HFs under higher fluid kinematic viscosity is shorter compared with the lower fluid kinematic viscosity models. The observed HF propagation behaviour in the isotropic stress conditions is also true for the models with the anisotropic confining stresses in Fig. 13b, but with the eminent directional preference in the direction of the maximum principal stress. Furthermore, under the anisotropic confining stress conditions, the length of the HFs is observed to be longer than those in the isotropic stress state due to the combined action of the fluid flow in the new HFs in the direction of σ H and the higher stress concentration in this direction. Therefore, under the higher fluid kinematic viscosities, the fluid pressure overshadows the high anisotropic in situ stresses effects, and the HFs propagate in all directions, reducing the HF treatment efficiency.
In the high fracture density DFN models under the isotropic confining stress condition, as shown in Fig. 14a, the low viscous fluid (2 mm 2 /s) flows quickly and causes rapid fluid pressure rise in the DFNs close to the borehole and also creates some new HFs, which facilitate the fluid flow through the DFNs. Under this viscosity, the fluid affected area is large. Increasing the kinematic viscosity to 3 mm 2 /s reduces the fluid flow speed which is marked by a smaller area under the fluid pressure raise. Once the kinematic viscosity is increased to 4 mm 2 /s, the fluid is unable to flow quickly through the DFNs and as seen in Fig. 14b, a tiny perimeter around the injection point is affected by the fluid flow and fluid pressure rise is limited to the vicinity DFNs. The HF patterns in the models with anisotropic confining stresses of σ H = 15 MPa and σ h = 5 MPa under various fluid viscosities are shown in Fig. 14b. The fluid injection with low viscosity of 2 mm 2 /s shows fewer new HF initiations from the injection point and the HFs preferably propagate towards the direction of σ H and also exhibit more activation and crossing types. However, the high viscous fluids (3 mm 2 /s) generate more new HFs around the injection point and severe breakage occurs at the immediate vicinity of the injection. However, the high viscous fluid is unable to flow further and fluid pressure rise in the direction of σ H is observed to be within a short distance. Therefore, low viscous fluid injection with a higher flow rate in deep hot rocks is recommended for an efficient fluid injection project.

Effect of DFN Aperture
Among the physical and mechanical properties of the DFN, the aperture plays a crucial role in the fluid path through the pre-existing fractures. It thus significantly influences the HF development pattern and efficiency. To study the effect of DFN aperture on the HF propagation pattern and extension, four numerical models with both the low and high fracture density DFN sets are built with the joint aperture varying from 0.2 to 0.3, 0.4 and 0.5 mm. The injection flow rate is fixed at 1 kg/s for all the models and the fluid kinematic viscosity is kept constant at 1 mm 2 /s. Once again, two confining stress scenarios are considered: isotropic and anisotropic confining stresses. Figure 15 shows the HF pattern for the low fracture density DFNs under both isotropic and anisotropic in situ stress states. As can be seen in Fig. 15a, under isotropic confining stresses, the HFs propagate in all directions with a similar pattern; however, the length of the HFs is considerably longer in the low aperture DFN model. Under the low DFN aperture, the fluid interaction with the DFNs decreases and fluid interaction increases with the rock matrix, resulting in  Temperature field in the models with low and high fracture densities faster rock breakage and, thus longer HFs. However, when the aperture increases, the influence of the DFNs on the fluid path increases and the fluid pressure accumulates in the DFNs, dilates them and causes more activation, crossing and offset crossing. Thus more fluid flow and energy are required to generate longer HFs. Under anisotropic confining stresses, as shown in Fig. 15b, the variation in the DFN aperture does not impose a significant change in the HF propagation pattern, and the high anisotropic confining stresses influence the fluid diffusion. Figure 16a illustrates the effect of the DFN aperture on the HF patterns under the isotropic confining stresses in the models with the high fracture density DFNs. As can be seen, the model with 0.2 mm aperture shows a larger perimeter with the high fluid pressure (marked with red and yellow contours) around the injection point. However, the fluid pressure rise is only observed up to a limited distance to the injection point (marked with blue contours). This behaviour is reversed with an increase in the DFN aperture. Once the DFN aperture increases, the area under the higher fluid pressure decreases and the area under the fluid flow extends further to a larger perimeter, as marked by blue counters of fluid pressure which is extended to all the DFN areas in the model with a 0.5 mm aperture. The reason is that when the aperture is small, the fluid can flow over a limited distance, which causes the fluid pressure to rise considerably inside the DFNs. Then, the fluid flow is stopped due to the small aperture and narrow channels. However, once the aperture increases, the DFNs can allow the fluid to flow easily inside the fractures over longer distances, decreasing fluid pressure rise in the vicinity of the injection point. Figure 16b depicts the HF patterns with varying DFN apertures under the anisotropic confining stresses. In these models, similar to the models shown in Fig. 17 under the isotropic confining stresses, the high fluid pressure builds up in the areas close to the injection point and the extension of the high fluid pressure region is larger for the low aperture DFNs. However, this extension is not uniform (different from the uniform extension of fluid pressure in Fig. 16a) and extends in the direction of σ H . Once again, the low aperture DFNs do not allow the fluid to flow over a large distance, while the DFNs with a higher aperture do allow the fluid to flow conveniently over a longer length due to the large fracture channels, which in turn avoids the high fluid pressure accumulation at the vicinity of the injection point. In Fig. 16b, it is interesting to note that under a low DFN aperture, the HFs propagate over a long distance in the direction of σ H while the high aperture DFNs case exhibits short HFs. This is because when the aperture is small, the fluid pressure rises quickly due to the difficulty in flowing in the narrow channels, and the high fluid pressure generates HFs in the direction of σ H with the high stress concentrating at the newly formed HFs, which further helps the HFs to grow faster. On the other hand, the DFNs with larger apertures provide wider flow channels where the fluid can flow easily and avoid high fluid pressure build-up in the areas close to the injection point.

Temperature Field
It is worth studying the temperature field distribution in these models. Figure 17 shows the temperature fields of some of the selected DFN models in Figs. 11,12,13,14,15,16 It can be seen that the tmeperature does not experience a significant change in the models with the DFNs. The reason is that since the matrix is impermeable and fluid is only flows through the newly formed HFs and also the DFNs, the cold fluid exchanges heat with the hot dry rock as soon as it flows into a fracture at the immidiate vicinity of the injection

Models with Blocks in a Matrix
Some reservoirs contain conglomerate rocks with blocks of different shapes and strength present in a matrix. These types of reservoirs are quite frequently formed under quick deposition in nearby origins including the slope belt of downfaulted basins (Feng et al. 2019). Based on their characteristics, generally tight conglomerate reservoirs exhibit weak porosity and permeability, but are usually associated with a high degree of matrix cementation and blocks compaction. With the rapid development of the HF stimulation methods, the conglomerate reservoirs have attracted much attention and contributed to significant production of unconventional oil and gas. These reservoirs are a supplement for other types of reservoirs. However, despite significant efforts on the characterisation of heterogeneous reservoirs, the significance of the conglomerate reservoirs has been largely ignored. In a conglomerate reservoir, the presence of blocks with different shape, strength and distributions results in intense heterogeneity, which strongly affects the HF propagation pattern and path. The complex interaction between the injecting fluid and a mixture of matrix-blocks induces an irregular and disordered HF path, which might be undesired and can reduce the gas and heat extraction rate. This behaviour can also cause unusual high fluid pressure during HF treatments, which in turn reduces the efficiency of the injection and increase the costs. In a HF treatment, the propagation of the HFs in a direction other than the desired path induces a great risk to the success of the project and, therefore, it is necessary to assess the influence of the presence of conglomerate rocks on the development of HFs. This section investigates the HF propagation in such formations.

HF Propagation in Conglomerate Reservoirs with Different Blocks
To investigate the influences of the blocks in a matrix in conglomerate reservoirs on the HF propagation, three sets of models are considered. Each model contains a single type of block. Here, three sets of blocks with various strength characteristics are simulated to investigate the effects of block strengths on the HF development pattern. The shape of the blocks is randomly built to reflect the real dimensions and shapes of a conglomerate rock mass. Table 5 presents the input data for the matrix and the blocks. It should be mentioned that the matrix is considered permeable while the blocks are considered impermeable. The permeability of the matrix is 3e-15 m 2 . Moreover, fluid can flow through the new HF. Mechanical interfaces are defined as the areas of the model where two different material properties are (or may become) in contact during the simulation and the mechanical properties of these interfaces can be explicitly defined. By default, the weaker property values of the two materials are used as the interface properties, although this choice can be overwritten by explicit assignment.
To better assess the HF development process in such rock masses, three models, each containing a single block type (block 1, block 2 or block 3, weak to strong) are built under an isotropic confining stress state of 6 MPa with a flow rate of 1 kg/s. Figure 18 illustrates the development stages of the HFs in the three models. In the model containing blocks 1 in Fig. 18a, at the early stages of the fluid injection, some blocks at the immediate vicinity of the injection point break and the new HFs cross them. Moreover, due to the low strength of the blocks, the working fluid not only breaks the blocks but also branches into the broken blocks to form more HFs. The HFs then grow further and break most of the blocks in their propagation path and, on some occasions, break the interface between the blocks and the matrix. In the models with medium strength blocks (block 2) shown in Fig. 18b, the working fluid breaks the interface between the blocks and matrix at the vicinity of the injection point and propagates further. In the propagation path, most of the fracturing is of intergranular type with some minor transgranular breakage within a few blocks. In models with hard blocks (block 3 in Fig. 18c), excessive hydraulic fracturing takes place around the injection point because the high strength blocks do not allow the working fluid to break them. As a result, the fluid pressure rises significantly in this area and breaks the matrix. In this case, no transgranular HFs are observed; all the HFs are of the intergranular type, and the interface breakage dominates the failure. Figure 19 shows the fluid pressure contours in the three models. As shown in Fig. 19a for the models with the low and medium strength blocks (block 1 and block 2), the fluid pressure contours exhibit a relatively uniform distribution around the injection point and the fluid pressure rise is limited to a small circumference. The medium strength blocks (block 2) arrest the fluid from flowing further and, therefore, fluid pressure rises around the injection point more than that of the block 1 case. When the blocks are hard as shown in Fig. 19a, c, very high fluid pressure accumulates Fig. 19 Fluid pressure contours in models containing a block 1, b block 2 and c block 3 after 1 s of fluid injection at the vicinity of the injection point and the hard blocks do not allow the fluid to cross the matrix or break them. The accumulated high fluid pressure causes excessive breakage of the matrix and the formation of intergranular HFs. As the result, the energy of the fluid is consumed in this area to break the matrix excessively and thus high fluid pressure spreads over a short distance. Figure 20 shows the contours of temperature around the injection point. For the low and medium strength blocks in Fig. 20a and b, the low-temperature fluid flows over longer distances and thermal exchange due to conduction heats up the fluid. In these models, the thermal conduction develops over a larger area because the low and medium strength blocks allow for easy fluid flow inside the blocks Fluid-matrix-block interaction types, a crossing, b interface breakage, c arresting, d interface breakage + crossing, e crossing + branching, f matrix branching, and g rock thermal shrinkage cracking and matrix, and heat exchange occurs frequently. The blocks with different thermal properties to the matrix cause some thermal cracking due to shrinkage deformation of the block, but their length is short compared to the main HFs (shown in Fig. 20d). In Fig. 20c for the hard blocks, since the hard blocks confine the fluid at the area close to the injection point, the thermal conduction occurs in a limited area, and the fluid heats up in this narrow area. In this model, the accumulation of cold water around the injection point leads to severe thermal shrinkage cracking, as shown in Fig. 20d. By comparing the fluid pressure and temperature fields in Figs. 19 and 20, the heat is only exchanged in the HFs and the low fluid pressure distribution in the matrix without HF generation does not influence the heat exchange. Therefore, the strength of the blocks plays a significant role in the heat exchange and conduction, which affects the efficiency of the HF operation in hot dry rocks.
Based on the HF paths observed in Fig. 18, the following interaction types between the fluid-block-matrix are identified and presented in Fig. 21. The crossing is the dominant interaction type observed in the models with low strength blocks (Fig. 21a). When the blocks are of high strength, interface breakage governs the interaction type (Fig. 21b).
Once the fluid has flown over a longer distance, a block may stop the further growth of a HF (Fig. 21c). In the models containing medium strength blocks (block 2), a propagating HF may cross the interface and break the block (Fig. 21d). In the blocks close to the injection point where the fluid pressure is high, the flowing fluid not only can cross a block but can also branch into more HFs inside the block (Fig. 21e). On some occasions, the fluid can branch within the matrix (Fig. 21f). The matrix branching is usually formed under two conditions: high flow pressure and the presence of high strength blocks, which do not permit the fluid to generate intergranular HFs. So fluid pressure rise in the matrix causes excessive branching. Lastly, the interaction between the cold fluid and hot rock causes rock shrinkage due to thermal shock, which results in the initiation of some short cracks from the matrix-block interfaces.

Effect of Flow Rate
The effect of flow rate on the HF propagation pattern is investigated in this section. Each model (models with block 1, block 2 and block 3) is simulated under two values of injection rate: 2 and 3 kg/s. Figure 22 shows the results of the fluid injection for each model. Figure 22a represents the HF pattern and fluid pressure in the low strength blocks (block1). To better visualize the HF patterns, an enlarged view around the injection point is cropped as shown in the yellow frame in Fig. 21a. Under the 2 kg/s flow rate, the newly formed HFs propagate uniformly in all directions and mainly break the blocks, while interface breakage occurs in a few blocks. Both the intergranular and transgranular HFs contribute almost equally. The fluid pressure distributes uniformly but with a smaller perimeter and a higher flow rate because a high flow rate consumes more energy to break more blocks and thus cannot flow longer distances. Under the high flow rate of 3 kg/s, the high pressure fluid breaks the matrix around the injection point Fig. 23 Influence of the fluid kinematic viscosity on HF pattern in the models containing a block 1, b block 2 and c block 3 after 1 s of fluid injection and excessive matrix branching is observed. Moreover, the transgranular HFs are seen to be the dominant cracking type with the blocks close to the injection point undergoing severe fracturing, reducing the HF treatment's efficiency. In the medium strength blocks shown in Fig. 22b, under both flow rates, severe matrix breakage occurs around the injection point and between some blocks, and matrix branching occurs quite frequently. However, under the injection rate of 2 kg/s, intergranular HFs dominate, while under the flow rate of 3 kg/s transgranular HFs govern the cracking type. In these models, the contours of fluid pressure show high pressure rise around the injection point with less development than those in the block 1 models. In the models with hard blocks in Fig. 22c, the fluid flow at the rate of 2 kg/s excessively branches the matrix around the injection point and intergranular HFs develop around the nearby blocks. The hard blocks do not allow the fluid to propagate over longer distances or break them. By increasing the flow rate to 3 kg/s, in addition to the severe matrix branching around the injection point and between the blocks, many intergranular HFs grow over longer distances, and a few transgranular HFs are observed due to the high energy fluid flow. In these grounds, the fluid pressure rises considerably high around the injection point, resulting in excessive matrix breakage and thus less fluid flow defuse.

Effect of Fluid Kinematic Viscosity
The three models are also examined under two different values of fluid kinematic viscosity, 2 and 3 mm 2 /s. Figure 23 presents the simulation results. In the low strength block model (block 1) shown in Fig. 23a, with the fluid of 2 mm 2 /s kinematic viscosity, the newly formed HFs propagate in all directions and break the majority of the blocks in their way with the intergranular and transgranular HFs contributing equally to create the HFs. Once the kinematic viscosity increases to 3 mm 2 /s, severe matrix breakage and branching occur around the injection point and the blocks at the vicinity of the injection point experience crossing, branching and interface breakage. In this case, the transgranular HFs occur more frequently than the intergranular one since the high viscous fluid is arrested around the injection point and is unable to flow further with sufficient energy to break the blocks; instead break the weakest part, which is the interface. In Fig. 23b for the model with medium strength blocks, the 2 mm 2 /s viscous fluid heavily breaks the matrix and some blocks around the injection point, while the 3 mm 2 /s viscous fluid excessively breaks and branches both matrix and blocks, resulting in a severe fragmentation zone. In the case of strong blocks in Fig. 23c under the 2 mm 2 /s kinematic viscosity, the fluid excessively breaks the matrix around the fluid injection point and the interface breakage occurs frequently, and the HFs propagate over a short distance. In this case, a few intergranular HFs are generated while the transgranular HFs dominate with most of the nearby blocks remaining unbroken. However, the 3 mm 2 /s viscous fluid breaks the matrix and the interfaces severely, and interestingly, many nearby blocks are also broken. In this regard, the developed distances of the HFs are short, and the high viscous fluid is unable to flow further. Furthermore, injecting the high viscous fluid causes more cold fluid accumulation due to the slow diffusion of the fluid, resulting in more heat exchange between the fluid, matrix and blocks and thus more thermal shrinkage cracking.

Effect of In Situ Stresses
It has been extensively studied and demonstrated that the in situ stresses (both direction and magnitude) have a  (Shi et al. 2022;Zhang et al. 2018). Based on the previous studies, it is evident that the HFs tend to propagate in the direction parallel to the maximum principal stress. The previous sections showed the HFs propagate in all directions with no directional preference when the confining in situ stresses are isotropic. To evaluate the effect of in situ confining stresses on the HF development behaviour in the blocky reservoirs, the three models (each containing a single type of block) are modelled under anisotropic in situ stresses of σ H = 15 MPa and σ h = 5 MPa. Figure 24 shows the HF pattern developed in the models with different blocks. In the models with low strength blocks (block 1) in Fig. 24a, relatively straight HFs propagate in the direction of σ H, and intergranular HFs dominate. In the model with block 2 shown in Fig. 24b, the direction of the HFs is observed to be almost parallel to the σ H , however, matrix branching occurs frequently in this direction due to the combined action of the fluid flow through the HFs and the high anisotropic in situ stresses. The model with block 3 in Fig. 24c shows that the high strength properties of the blocks suppress the high anisotropic in situ stresses effect, and the nearby blocks deviate the fluid to propagate in all directions with a minor tendency in the direction of σ h , with interface breakage being the dominant HF type. Moreover, once the strength of the blocks increases, the fluid accumulation occurs more in the matrix causing high fluid pressure to rise in the area. Furthermore, this fluid accumulation facilitates better heat exchange, thus, more thermal shrinkage cracking. The contours of the fluid pressure demonstrate that block strength increase causes excessive fluid pressure around the injection point (23.8 MPa for block 3 compared to 16.8 MPa for block 1) and less fluid pressure distribution over longer distances, indicating inefficiency of the hydraulic fracturing treatment. Figure 25 illustrates the temperature field for the blocky reservoirs. It should be mentioned that the models in Fig. 25 are the temperature fields of the models shown  Fig. 25a, b, the HFs propagate over long distances in the matrix containing the low to medium strength blocks (block1 and block2), and some blocks are fragmented as well. Thus the heat exchange between the fluid inside the HFs and the matrix occurs frequently. The cold fluid far from the injection point is also heated to the reservoir temperature. When the block strength is high (block 3) in Fig. 25a--b, the cold fluid is suppressed in the area close to the injection point due to the increased resistance of blocks against the fluid flow. Therefore, the low-temperature field is only limited to the area close to the injection point and the temperature in the long HFs is 200 °C due to the rapid thermal conduction and convective heat transfer. A similar temperature field is observed for the models with high kinematic viscous fluid, as shown in Fig. 25c, d. However, as the kinematic viscosity is a measure of the fluid's resistance to deformation and movement, this fluid cannot flow for long distances, meaning the fluids accumulates in the limited zone around the injection point. Therefore, the fluid pressure rises significantly, and the breakage is excessive, resulting in frequent heat exchange in this area. Thus, the fluid in longer distances is a heated fluid; therefore, the low temperature is only observed in the small perimeter around the injection point. Under the high anisotropic confining stresses shown in Fig. 25e, the strength of the blocks significantly impacts the temperature field. In particular, the low to medium strength blocks allow the cold fluid to flow over a longer distance and transfer more fluid for heat exchange. However, the high strength blocks suppress the fluid flow from creating long HFs, and thus the convective heat transfer is seen to happen in a narrow area. Overall, the temperature rise in the fluid due to the conduction between the hot blocks and the fluid is evident in the high strength blocks model, and shrinkage cracking often occurs close to the injection point. In the model with block 3, the thermal exchange between the high strength blocks and strong conduction causes frequent shrinkage cracking and less high-temperature fluid distribution over longer distances. Since the strength of block 1 is low and the fluid can break them and flow over longer distances, the fluid temperature rise is mainly through the main HFs and shrinkage cracking occurs less often.

Discussion
The distribution of the DFNs or blocks in a model and varying the controlling parameters in the previous sections showed considerable effects on the HF development pattern, fluid pressure distribution and temperature field. Under the variables mentioned above, it is worth looking at the type of the HFs: tensile or shear. It was observed from the models with DFNs that all the HFs in these models are of tensile mode and, therefore, are not shown here for briefness. Since the spatial distribution volume of the DFNs is much lower than the matrix volume, the matrix significantly dominates the hydraulic fracturing path; in this case, the HFs are all in tensile mode. This finding is confirmed by experimental tests such as those in Sun et al. (2021) and Wu et al. (2019). Moreover, we only plotted the failure modes of HFs to make comparisons. Therefore, only tensile or shear modes of HFs were presented. In the utilized software, arbitrary fracture trajectories are free to develop within the constraints imposed by the mesh topology. Depending on the local stress and relative crack wall displacements, the crack elements may yield and break under Mode I (i.e., opening mode), Mode II (i.e., sliding mode) or mixed-Mode I-II conditions. The DFNs may experience opening or slipping when interacting with the HFs, however, we did not consider them in this study since the focus was to investigate HF development patterns and failure modes. However, the blocky models show considerable failure mode change under different blocks' strengths or controlling factors. Figure 26 summarises the failure modes of the blocky models under the various conditions discussed. For the block 1 models, no shear crack is observed under any flow rate or viscosity, and the failure mode is pure tensile. Once the strength of the blocks increases to block 2, shear cracks develop mainly in the interface between the matrix and the blocks, and tensile fractures dominate the failure of the matrix and branching inside the broken blocks. The high-strength blocks (block 3) exhibit tensile failure dominated region around the injection point with transformation to shear cracking away from the centre. In these models, shear failure dominates since the HFs tend to break the interface rather than the matrix or the blocks, and the matrix breakage occurs easily under the shear mode loading. Under the high anisotropic confining stresses, as shown in Fig. 26e, the large difference between the magnitude of the in situ stresses help to break both the matrix and the interfaces; therefore, both tensile and shear crack occurs in the block 1 and block 2 models. The high strength block 3 model again shows more shear cracking due to the higher resistance against the fluid flow and the breakage of interfaces. Figures 27 and 28 show the fluid pressure at the injection point for the models with low (corresponding to Figs. 11, 13 and 15) and high-density DFNs (corresponding to Figs. 12, 14 and 16) under the mentioned controlling factors. It can be observed that increasing the flow rate or fluid kinematic viscosity causes considerable fluid pressure rise, which indicates more HF development. The fluid pressure histories for the permeability change in the low fracture density DFNs (Fig. 27c) show that a difference in the DFNs permeability does not significantly impact the fluid pressure distribution. It also offers a minor influence on the fluid pressure rise in the models with high fracture density DFNs (Fig. 28c). Figure 29 shows the fluid pressure curves at the injection point for the blocky grounds corresponding to Figs. 22, 23 and 24. As can be seen in these graphs, the strength of the blocks significantly affects the fluid pressure histories. The diagrams demonstrate that the higher the strength of the blocks is, the higher the fluid pressure rise is at the injection point. Under the high anisotropic confining stresses, as shown in Fig. 29c, block 1 and block 2 show similar fluid pressure behaviour, while block 3 exhibits a distinctive difference due to the blocks' resistance to fluid flow.
Since the study's main objective was to study the interaction of fluid-driven fractures with the sources of heterogeneity like DFNs and blocks, the simulations were terminated after a few seconds of fluid injection. This injection time was long enough to allow HFs to propagate to a proper length, so extensive interaction between the newly formed HFs and the DFNs or blocks was guaranteed. This injection time and fracture extension length showed that all possible interactions could be captured and that the fluid had enough time to exchange heat with the hot rock. Moreover, the models were large enough to avoid any boundary effects, and further extension of HFs to the boundaries could impose drawbacks which may have led to unrealistic results. Furthermore, the computational time of each model was very long, and this long run time was avoided by terminating the runs when a suitable and meaningful HF extension was achieved. Also, allowing more extension of HFs would not provide any new findings since the main interaction in the fluid path occurred close to the injection point.

Conclusion
This study provides some general recommendations for a successful HF treatment in complex EGS. The paper presents the hydraulic fracture (HF) simulation results in complex reservoir conditions containing low and high fracture density DFNs and blocks. The combined finite-discrete element method (FDEM) is employed for the thermo-hydromechanical (THM) simulations in deep hot dry rocks since this method has proven to be capable of modelling fracture propagation in rocks under complex conditions. First, the validity of the FDEM is verified against existing analytical solutions. The influence of some controlling factors is studied on the HF pattern in both ground conditions, and the following main results are drawn: 1. When the fracture density of the DFNs is low, the controlling factors such as flow rate and fluid kinematic viscosity have less effect on the HF pattern and more on the rock matrix. Higher flow rates generate more HFs regardless of the direction of the maximum principal stress and high kinematic viscosity excessively breaks the rock around the injection point. 2. In reservoirs containing high fracture density DFNs, the HFs pattern is strongly controlled by the DFNs rather than the fluid parameters. In these reservoirs, the higher flow rate causes the HFs to propagate in a specific direction that the connecting DFNs govern. The high viscous fluid is prevented from further flowing since the high fracture density DFNs act as a barrier. 3. Once the DFNs aperture increases, not only the area with a higher fluid pressure decreases but also the fluid flow affected area extends further to a larger perimeter, because the fluid can flow over a limited certain distance,causing the high fluid pressure inside the DFNs.
The fluidthen is stopped by further flow due to the small aperture. The DFNs with high aperture can allow the fluid to flow easily inside the fractures over longer distances, decreasing fluid pressure in the vicinity of the injection point. Larger aperture values promote higher flow velocities and, therefore, shorter residence times, which in turn hinder the solid-fluid heat exchange process. 4. In the conglomerate reservoirs, six interaction types are identified between the working fluid and the blocks: crossing, interface breakage, arresting, interface breakage + crossing, crossing + branching, and matrix branching. 5. In conglomerate reservoirs, the strength of the blocks strongly affects the HF development pattern and fluid pressure distribution. In the reservoirs with soft blocks, the fluid breaks the blocks at the vicinity of the injection point and flows easily over longer distances, However, when the blocks are stronger, the fluid is unable to break them and accumulates around the injection point leading to excessive break of the matrix. Also, since the further flow of the fluid is stopped, fluid pressure rise is severe and interface breakage occurs frequently. 6. In the blocky reservoirs, due to the high contrast between the thermal properties of the matrix and blocks (mainly linear thermal expansion coefficient), temperature fields follow a complex trend in which crack branching frequently occurs from the interface of the matrix and blocks. The crack branching occurs due to the excessive thermal exchange between the matrix and blocks, facilitating fluid temperature rise and better HF network generation. 7. In the reservoirss with DFNs, all the HFs are observed to be of tensile mode due to the matrix breakage being the dominant failure path. In the blocky reservoirs, the models containing low strength blocks exhibit pure tensile cracking while the ones with higher strength blocks show the development of shear cracks due to the interface breakage.
Funding Open Access funding enabled and organized by CAUL and its Member Institutions.

Data availability
The data that support the fingdings of this study are available from the corresponding author upon reasonable request.