A novel radial jet drilling stimulation technique for enhancing heat recovery from fractured geothermal reservoirs

In this study, the application of the Radial Jet Drilling (RJD), a novel stimulating technique for enhancing productivity in the existing wells in deformable naturally fractured reservoirs was investigated using a robust three-dimensional ﬁ nite element DFM (discrete fracture-matrix) model. Results showed that the RJD laterals were more effective in enhancing injectivity/productivity in cases with lower fracture density, i.e. lower equivalent permeability, while they had no signi ﬁ cant effect on the heat production in these cases. In higher fracture density cases, the RJD laterals improved the heat production while had no signi ﬁ cant effect on the injectivity/productivity. Results also showed that in reservoirs with very low permeability matrix, the RJD laterals can be used to connect the wells to the fracture network and hence enhance the well performance. The sensitivity analysis on the average net energy production rate with respect to the length of the RJD laterals showed that in the situations where the wells were not connected directly to the fractures, the length of RJD laterals played a crucial role in enhancing the average net energy rate. However, the 100 m laterals almost removed the dependency of the average net energy production rate on the well placement for low, medium and high fracture density cases. © 2019 Elsevier Ltd. All rights reserved.


Introduction
Fossil fuels are still the main contributor in energy market around the world, however, concerns regarding the security of energy, climate change and carbon emissions have encouraged developing strategic plans towards low carbon future by relying on renewable energy resources (Bruckner et al. [52]). Geothermal energy stored in the Earth's crust is one of the promising and clean renewable energy resources in the world [1,2]. Geothermal energy has the potential to integrate with other sources of energy to improve the overall energy production and make the heat production economically attractive. For instance, the integration of geothermal energy with heavy oil production [3] is proposed. Today, geothermal energy is extracted from the earth's crust by drilling injection and production wells, and the so-called doublet system in which cold fluid (i.e. water) is injected at the injection well and the hot water or steam is produced from the production well, is very popular in geothermal sites around the world [4]. The geothermal systems rely on the injected fluid being able to flow from the injection well to the production well through the reservoir, either using the natural permeability of the formation [5,6] or fractures in the formation as flow paths [7]. In deep geothermal reservoirs, the formations are typically of very low permeability, and fractures, natural or man-made, are needed to enhance the fluid flow within these reservoirs.
In the enhanced geothermal systems (EGS), the natural permeability of the reservoir is enhanced by creating highpermeability flow paths, and many projects are carried out around the world including the Fenton Hill EGS in the United States, Soultz EGS in France, and Cooper Basin EGS in Australia [8]. The conventional method for creating the flow paths is the hydraulic fracturing (i.e. fracking), however, due to the uncertainties as well as the seismic activities associated with fracking, the application of this method has been restricted in Europe [9]. For the hydraulic fracturing to become an option for stimulation, there is a need for the industry to demonstrate to government and the community that hydraulic fracturing processes pose a minimal risk to people and the environment. Recently, a new well stimulation method, known as Radial Jet Drilling (RJD), has become an attractive technology for enhancing well productivity/injectivity in low performing geothermal wells [10]. The method has been frequently used in oil wells to bypass the damaged zone [11] or to connect the well to natural fractures as shown in Fig. 1. The RJD stimulation technique has several benefits over the conventional hydraulic fracturing technique, including negligible chance of inducing seismicity, lower amount of water needed, more controls on the direction and the length of the laterals, and most importantly, it is cost effective to perform in existing wells. Thus, the RJD is considered as a viable alternative to conventional hydraulic stimulation technologies for enhanced or engineered geothermal systems, and has triggered several research projects around the world (e.g. Refs. [10,12]). Within the geothermal fields, recently this technique was performed in a low performing injection well in Klaipeda, Lithuania, leading to an increase of injectivity of about 14% [13]. For the RJD method to be successful, a minimum porosity of about 3e4% of the formation is required [14]. The location of the jetting nozzle in the RJD technique can be estimated by monitoring the acoustic emissions generated by the jetting process as shown by Ref. [15] in a jetting experiment in a quarry.
Considerable efforts have been expended in developing thermal-hydraulic-mechanical models (THM: e.g. Refs. [7,16,17]) and together with chemical reactions (THMC: e.g. Refs. [18,19,53]) for geothermal reservoirs. Numerical models are essential tools for evaluating the performance of a potential geothermal reservoir and for studying the effects of different parameters on the productivity of the reservoir [20e23,54]. Fractures can have a significant impact on the reservoir response, and the energy production, thus, their accurate geometrical representation in the numerical model, as well as their dynamic behaviour under varying reservoir state variables (pressure, temperature and stress) is crucial in making realistic predictions [24]. In addition, fractures can contribute to the short-circuiting between the injector and producer wells, a common problem in geothermal reservoirs that affects the energy recovery from the system. In short-circuiting, a high-conductivity channel is created between the injector and producer wells (mainly through the fractures), preventing the cold fluid from accessing the bulk rock where the heat is stored. The RJD laterals can also contribute to the short-circuiting problem by connecting the well to the "wrong" fracture or mitigate the problem by creating new paths for flow. Furthermore, the local increase in the fracture (a) RJD laterals created from the main wellbore to bypass the damaged zone around the well, (b) the RJD lateral is used to connect the well to the natural fracture, (c) the RJD lateral is used to improve the access of the well to the natural fracture, (d) the cross-section of the RJD lateral drilled with water using static nozzle in chalk, (e) the cross-section of the RJD lateral drilled with water using rotating nozzle in chalk, (f) the cross-section of the RJD lateral drilled with acid using static nozzle in chalk, (g) Static nozzle: 1-nozzle body, 2-forward orifice, 3-backward orifice; (h) Rotating nozzle: 1-static inner body, 2-rotating body, 3-static outer body, 4-forward orifice, 5-rotating body orifice, 6-backward orifice. Nozzles in the pictures are in mm-scale. aperture due to the combined poroelastic and thermoelastic deformations of the rock matrix can facilitate the creation of shortcircuits [19,25e29]. Although several recent studies have investigated the utility of RJD laterals for well stimulation and well flow enhancement (e.g. [30]), to the best knowledge of authors, a comprehensive study of the application of RJD technique in fractured geothermal reservoirs is missing and the present study aims to fill the gap. In this study, the application of the RJD laterals in enhancing injectivity/productivity of the geothermal reservoirs as well as their effects on the heat production is investigated using several cases of fracture densities (low, medium, and high), fracture geometries (short and long fractures), well placements (whether the well is connected to the fractures or located away from the fractures), and the RJD laterals length (zero, 25 m, 50 m, and 100 m). The numerical analyses are performed using a robust three dimensional finite element method, capable of simulating non-isothermal flow through the wells, RJD laterals, fractures and rock matrix. The sensitivity of the net energy production rate, the production temperature, and the pumping energy to the lateral length, orientation and their connection to the fracture network is studied. The effects of matrix deformation, aperture variation, matrix permeability, and hydraulic fracture development from the RJD laterals on the production temperature and injectivity/productivity are also investigated and the results are presented. Finally, the net energy production rate from the wells is calculated using the thermal energy produced minus the pumping energy consumed. The sensitivity of the average net energy production rate to the laterals length, well placements, and the fracture densities is investigated and the results are presented.

Radial jet drilling (RJD) technology
In the RJD technique, the power of a fluid jet (commonly water, sometimes mixed with acid) is utilised to drill small diameter holes from the main wellbore as shown in Fig. 1. RJD laterals can reach up to 100 m of extension [31], and in different directions from a single depth, thus, they can effectively be used to bypass the damaged zones around the wells, connect the well to the nearby fractures, and improve the efficiency by creating conductive flow paths between the reservoir and the wellbore [11]. The method is relatively fast, cost efficient, and could be used as an alternative to the conventional well stimulation techniques, especially in layered reservoirs and in formations that are close to the water contact.
The process of radial jet drilling is relatively simple: a bottomhole assembly, generally referred to as 'deflector shoe', is connected to a tubing and lowered to the target depth. For cased hole intervals, a coiled tubing conveyed milling assembly is lowered through the tubing. At the bottom, the deflector shoe deflects the milling bit towards the casing. After milling a hole into the casing, a jetting assembly is lowered through the tubing. The coiled tubing conveyed jetting assembly consists of a self-propelled jetting nozzle attached to a flexible hose [15]. Commonly, two different nozzles are used for RJD as shown in Fig. 1: static nozzle and rotating nozzle. The static nozzle is made of one body, which is connected to the high pressure hose with a thread connection. It has forward and backward orifices. The forward orifices mainly serve for eroding and breaking the surface of the formation, while the backward orifices widen the hole as well as push the nozzle forward. The diameter of each forward orifice is 0.15 mm and the diameter of the backward orifice is 0.4 mm, and 15e20 l/min flow rate provides exit velocities of about 358e478 m/s [32,33]. The rotating nozzle has three main components, the inner and outer static body and the rotating part around the static body. The rotating nozzle also has forward and backward orifices with similar functionality as in the static nozzle. The diameter of the forward orifice is 0.3 mm, while the diameter of the backward orifice is 0.5 mm, and 15e20 l/min flow rate provides exit velocities of 171e228 m/s [32]. Commonly water is used as the jetting fluid; however, in carbonate rocks such as chalks, acid-aided water is proven to be more effective [34]. The cross-sections of the RJD lateral drilled in chalk with water and acid using static and rotating nozzles are shown in Fig. 1. In some cases, water is not an option to use, for example in water-sensitive formations due to the risk of swelling. In such cases, diesel is a natural alternative to use. The diameter of the hole varies between 2 and 5 cm, depending on several factors including the size and type of nuzzle, the strength properties of formation rock, the jetting fluid, and the far field in situ stresses [34,35]. Tight chalk has shown resistance against drilling with water using the static nozzle, but some penetration has been achieved using rotating nozzles. Rotating nozzles generate the highest rate of penetration (ROP) values in both laboratory and quarry experiments [15,34].
An important aspect of the RJD technology is the sustainability of the increase in production. Having a small diameter, the RJD laterals cannot be protected by casing, thus they are susceptible to collapse under reservoir conditions. Observations from field cases show that shortly after jetting, the production increases significantly, but it suffers from the reduction in long term [36]. This may occur due to loss of strength in the formation rock; therefore it is likely that the RJD laterals have a limited life time [37,38]. After the drilling, the in situ stresses are redistributed with concentrations near the lateral's walls. If the redistributed stresses exceed the rock strength, instability issues may take place. At injection wells, the temperature difference between the fluid in the well and the surrounding formation induces thermal strains which affect both the axial and tangential stresses [39,40]. At production wells, the high inflow rates into the laterals can cause erosion and fine production, resulting in enlarged hole, hence triggering mechanical instability and eventually the hole collapse. Stability analyses have shown that small diameter RJD laterals are relatively stable under static loads i.e. immediately after drilling, but they may not be stable under dynamic reservoir conditions in the long term [40].

Computational model
In the present study, fractures are modelled as discrete surfaces in the three-dimensional matrix. The single-phase flow through fractures is defined using lubrication equation, while single-phase flow through the matrix follows the Darcy's law. The elastic deformation model is expressed satisfying the condition of equilibrium on a representative elementary volume (REV) of the porous medium as divðDε À ap m I À b s KðT m À T 0 ÞIÞ þ F ¼ 0 (1) where, D is the drained stiffness matrix, ε ¼ ðVu þ Vu T Þ=2 is the strain, u is the displacement vector, a ¼ 1 À K=K s is the Biot coefficient of poroelasticity, K is bulk modulus of rock, K s is bulk modulus of the solid grains, b s is the volumetric thermal expansion coefficient of the rock matrix, T m is the matrix temperature, T 0 is the initial temperature, I is the second-order identity tensor, F is the body force per unit volume, and p m is the fluid pressure in the matrix. Fluid flow through deformable matrix is expressed combining the mass conservation with Darcy's law for a single- where, k m is the permeability tensor of the matrix, g is the vector of gravitational acceleration, m f is the fluid viscosity, f is the porosity, c f and b f are coefficients of the fluid compressibility and volumetric thermal expansion, respectively. The governing equation for heat transfer through the rock matrix can be obtained by combining Fourier's law with an energy balance for saturated rock. It is assumed that the fluid velocity in the rock matrix is slow enough such that the solid grains and the fluid in the rock matrix are always in local thermal equilibrium. Convective, i.e., conduction and advection, heat transfer in rock matrix can be written as where, v m is the fluid velocity in matrix, l m ¼ ð1 À fÞl s þ fl f is the average thermal conductivity tensor of the matrix, l s is the thermal conductivity tensor of the rock, r m is the average density of the matrix (saturated rock), c hm is the average matrix heat capacity such that r m c hm ¼ ð1 À fÞr s c hs þ fr f c hf , r s is the rock density, c hs is the rock heat capacity, l f is the thermal conductivity tensor of the fluid, r f is the fluid density, and c hf is the fluid heat capacity.
Assuming a high aspect ratio fracture that has a lateral extent much larger than its aperture, the fluid flow through deformable fracture can be written as div where, p f is the fluid pressure in the fracture, and a f is the fracture aperture. Finally, the governing equation for heat transfer through the fluid in the fracture can be obtained by combining Fourier's law with an energy balance for the fluid. The advective heat transfer through the fluid in the fracture can be written as div where, T f is the fluid temperature in the fracture, v f is the fluid velocity in fracture. The continuity of pressure and temperature across fractures is satisfied by assuming T f ¼ T m and p f ¼ p m on the fractures. More details on the computational model can be found in Refs. [7,28]. The governing equations are approximated numerically using the finite element method. Spatial and temporal discretisations are realised using Galerkin method and finite difference techniques, respectively. These techniques are geometrically flexible so that complex intersections of fractures and the resulting complexly shaped matrix blocks can be discretised with small elements, while larger parts of the rock matrix can be meshed coarsely, (cf. [41]).
The displacement vector u, the matrix fluid pressure p m , and the rock matrix temperature T m are taken as primary variables. Using the standard Galerkin method, the primary variable X ¼ fu; p m ; p f ; T m ; T f g within an element is approximated from its nodal values as where N is the vector of shape functions and b X is the vector of nodal values. Using the finite difference technique, the time derivative of X is defined as where X tþdt and X t are the values of X at time t þ dt and t, respectively. Implicit backward Euler time discretisation technique has been used to develop the set of discretised equations in matrix form of SX ¼ F, in which S is the element's general stiffness matrix, and F is the vector of right-hand-side loadings. The components of the stiffness matrix are dependent upon the primary unknown variables, i.e., conductance, capacitance and coupling coefficients of the fracture are all dependent on the fracture aperture; therefore, a Picard iteration procedure is adopted to reach the correct solution within acceptable tolerance. For the current iteration, s þ 1 in the current step, n þ 1, the solutiondependent coefficient matrices in the stiffness matrix S are updated using weighted average solution vector X sþq nþ1 defined as where X sÀ1 nþ1 and X s nþ1 are the solution vectors of the two most recent iterations in the current timestep n þ 1, and q ¼ 2=3 is the weighing coefficient. For the first iteration s ¼ 1, the previous timestep solution is used as where X n is the solution vector from previous timestep n. The iterations are repeated until consecutive normalised values of X s nþ1 agree to within a specified tolerance ε The tolerance is set to 0.01 in this work. The discretised equations are implemented in the Complex Systems Modelling Platform (CSMPþþ, also known as CSP), an object-oriented application programme interface (API), for the simulation of complex geological processes and their interactions (formerly CSP, cf. [42]). Unstructured elements are used for spatial discretisation of surfaces (triangles) and volumes (tetrahedra). Fracture flow, and heat transfer equations are solved on the surface elements, whereas, matrix deformation, fluid flow and heat transfer equations are solved on the volume elements. The ensuing set of linear algebraic equations is solved at each timestep using the algebraic multigrid method, SAMG [43]. The proposed numerical model is validated against several analytical solutions and other numerical simulations. Validation examples are performed in other publications by the authors (e.g. Refs. [7,28].

Well model
In this study, the wells and RJD laterals are modelled as 1D line objects. Since the flow into a wellbore is mainly a radial flow, the problem of singularity arises, in which, the region where pressure gradients are the largest is closest to a well and is far smaller than the spatial size of elements. Using local mesh refinement around the well can alleviate this problem at the cost of simulation time [30,44]. The well model computes flow into the wellbore accurately using accurate well equations that allow the computation of the bottomhole pressure when a production or injection rate is given or the computation of the rate when the bottomhole pressure is known. In this study, the well flow is defined as [45].
where, k 1 and k 2 are matrix permeabilities in two orthogonal directions perpendicular to the wellbore, h is the height of the element, r e z0:2h is the equivalent radius at which the steady state flowing pressure for the actual well equals the numerically computed pressure for the well element [46], r w is the radius of the wellbore, p w is the well pressure (bottomhole pressure), and p m is the fluid pressure at the matrix. A laminar flow through the wells and laterals is assumed as Using the well model, the dependency of the pressure/flow rate on the element size has been significantly reduced as shown in Fig. 2 for a simple well test. As shown in this figure, without using the well model, results for the well pressure corresponding to a constant production rate are significantly mesh dependent. The simulation results approaches the exact (analytical) solution as the size of the elements at the well reduces. However, using the well model, the mesh dependency has been removed and the pressure results for different mesh sizes are very close to the exact solution. Using the well model is crucial for small-diameter RJD radials to capture the accurate radial pressure/flow rate and its overall contribution to the production of the main wellbore. Also, in Fig. 2, the results for the fluid temperature at production in a simple doublet system are presented and compared for two cases with and without the well-model. The application of the well-model in two cases is compared: constant pressure injection/production versus constant flow rate. When the bottom-hole-pressure (BHP) is defined, again the well flow rate will be mesh dependant, thus the fluid temperature at production can be erroneous. However, when the flow rate is defined at the wells, the error due to lack of using a well model is limited to the BHP and the pressure in vicinity of the wells, while the heat transfer and the production temperature are not affected significantly as shown in Fig. 2.

Simulation models and scenarios
Heat production is simulated through a combination of an injection well and a production well, i.e. a geothermal doublet. The geothermal reservoir is consisted of a single layer of thickness 100 m and lateral extent of 3 km by 3 km as shown in Fig. 3. The wells fully penetrate the layer. Two types of fractures: short and long are considered. For short fractures, three cases of fracture densities are considered: low, medium and high densities as shown in Fig. 4. The fractures extend fully through the layer thickness. Fracture network extension may follow the fractal law [47]. In this study, fractures are created using a novel method in which the propagation of fractures are simulated dynamically using the equations of motion, and mass and momentum conservations, and include the interactions between different fractures (Welch and Luthje [55]). For each short fracture set, three well placements are considered: 1) the injection and production wells are located in the rock matrix and away from the fractures, 2) the injection well is located close to the fractures and the production well intersects a fracture, 3) both wells intersect the fractures. The spacing between wells in all cases is kept constant and equal to 1017 m. For the long fracture case, the wells are located in the matrix and away from the fractures.
Four RJD laterals in both injection and production wells are considered. The laterals are positioned horizontally in the centre of the layer. For the long fracture case, the size and direction of laterals are varied to investigate the connectivity of the laterals to the fracture network. RJD laterals are a good candidate for inducing hydraulic fractures. In this case, the RJD laterals are pressurized to create hydraulic fractures, and since there are many laterals drilled in one plane, there is a good chance for the induced hydraulic fractures to intersect and create a larger plane in the desired direction. In short fracture cases with well placement 1 (wells are located away from the fractures), the effect of using a hydraulic fracture in place of the RJD laterals is investigated. It is assumed that a circular fracture of radius 100 m is induced through the RJD laterals as shown later in Fig. 9. The well placement 1 which is located on the matrix is selected for the hydraulic fracturing case, as the other cases would develop uncertainties in the induced fracture geometry due to the presence of the natural fractures connected to the induced fracture.
In all simulations, the injection and production rates are set to Q ¼ 0:02 m 3 /s, the initial reservoir pressure is set to 20 MPa, the initial temperature is set to 80 C, and the isotropic in situ stress is set to 35 MPa. Water is injected at a temperature of 20 C. The model parameters are given in Table 1. The matrix is considered permeable with permeability varying between k m ¼ 1 Â 10 À15 m 2 to k m ¼ 1 Â 10 À13 m 2 and the fractures initial aperture is set to a uniform value of 1 mm. The equivalent permeability of the low, medium and high fracture density cases with matrix permeability of k m ¼ 1 Â 10 À13 m 2 are calculated as 3.6 Â 10 À13 m 2 , 6.5 Â 10 À13 m 2 , and 1.6 Â 10 À12 m 2 , respectively. Since the fractures are distributed relatively uniform, the equivalent permeabilities in both horizontal directions are similar. For long fractures, the matrix permeability is set to k m ¼ 1 Â 10 À14 m 2 . The left and right boundaries are considered open for all cases, while other boundaries are set to no-flow boundaries.
For the deformable matrix, the classic Barton-Bandis model [48,51] is utilised to express the fracture aperture in terms of the contact stress as a f ¼ 0:011 À 7:33 Â 10 À9 s n 1 þ 6:67 Â 10 À7 s n (13) in which, s n in the normal component of the contact stress over the fracture. For short fracture cases with deformable matrix, the Young's modulus and Poisson's ratio are set to 20 GPa and 0.25, respectively [28]. While for the long fractures the Young's modulus and Poisson's ratio are set to 50 GPa and 0.25, respectively [7]. The initial aperture under initial in situ stress (35 MPa) and pressure (20 MPa) is 1 mm which is consistent with the rigid-fracture simulations [27]. For short fractures, the deformable case is considered for well placement 3 in which both wells are connected to the fractures. The classic Barton-Bandis model is chosen for this work as it is commonly employed for simulation of variable fracture aperture [27,49]. All the simulations are conducted for 30 years as the lifetime of geothermal plants is commonly designed for a maximum of 30 years [4,20].

Results and discussion
In this section, the application of RJD laterals in heat production from different geothermal reservoirs is investigated. The different cases vary in the density of the fractures, size and orientation of RJD laterals, connectivity of the laterals to the fracture, matrix permeability, as well as the size of the fractures. The simulations are performed for 30 years using an initial timestep of 1 day which increases by a factor of 1.2 every step until it reaches to a prescribed maximum timestep of 0.25 years. First the results for the short fracture sets are presented and discussed, followed by the results for long fracture set. Finally, the net production energy rate for several cases of short fractures are calculated and its sensitivity to the length of the RJD laterals are investigated.

Short fractures
The temperature distributions in the matrix after 30 years of simulation for well placement 3 with RJD laterals for different short fracture densities are shown in Fig. 5. The matrix permeability is set to k m ¼ 1 Â 10 À13 m 2 in these cases. The temperature distribution in the low-density case resembles a smooth circle around the injection well, and the fractures do not affect the temperature distribution significantly. In the high-density case, however, the fractures are the primary flow paths, thus, disturbance of the temperature is observed. The RJD laterals in the injection well are drown into the cold plume, but in the production well two laterals connect the well to the cold plume while the other two provide access to the hot matrix behind the production well. In particular, the laterals in medium and high-density cases connect the production well to a network of fractures behind the production well that provide access to the hot matrix. The pressure (injection and production) and pressure difference between the injection and production wells (DP) during 30 years of simulation for all cases of short fractures are shown in Fig. 6. At the injection well, the buildup pressure is reduced with increase in the fracture density, and the pressure drawdown at the production well is also reduced by increasing fracture density. This is evidently a sign of increase in the equivalent permeability of the reservoir due to higher fracture density. The placement of the well has a significant effect on the pressure in the low-density case with little effect on the pressure in the high-density case. In the low density case, when the placement of the wells move from the matrix to the fracture it shows reduction in the build-up and drawdown pressures at the injection and production wells, respectively. Installing RJD laterals also improves the injectivity and productivity of the wells by reducing the build-up and drawdown pressures. The effect of the RJD laterals is more dominant in the low fracture density case and less dominant in the high fracture density case.
The production temperature profiles in Fig. 7 show that high fracture density cases, while providing better injectivity/productivity, create a shorter path for the cold fluid to reach the production, therefore, the production temperature declines faster in the high fracture density cases than the low fracture density cases. Interestingly, the RJD laterals have no significant effects on the production temperature from the low fracture density cases (where the RJD laterals have the most injectivity/productivity enhancement), while the RJD laterals have positive effects on the production temperature from the high fracture density cases (where the RJD laterals have negligible or low injectivity/productivity effects). The results for the medium density cases are mixed, for instance the temperature drawdown is faster for well placement 2 (Fig. 7b), but it is slower for well placement 3 than the other two cases (low and high densities). This is due to the orientation and position of fractures around the production well as shown in Fig. 8 for well placement 2. In low density case, the production well intersects an isolated fracture, in the medium density case, the well is connected to a fracture network in front of the production well which acts unfavourably for the heat by creating a shortcut to the injection well, and in the high case, the well is connected to a fracture network behind the production well which acts favourably for the heat by accessing the hot matrix. Thus, for well placement 2, the medium density case suffers from faster temperature decline in the production as shown in Fig. 7b, while the high-density case provides better heat production. In well placement 3 in medium case, the production well is connected to a favourable fracture, so the temperature drawdown is slower as shown in Fig. 7c.
The pressure, pressure difference and temperature evolution of the hydraulically fractured cases are also included in Figs. 6a, d and 7a, respectively, shown by dotted lines. The fractured RJD laterals show better injectivity and productivity in the low and medium density fracture cases compared to the RJD laterals, but show little improvement in the highly fractured case. The effect of the hydraulic fracture on the production temperature is negligible for the low fracture density case, positive for the medium fracture density case, and negative for the high fracture density case compared to the RJD cases. By inspecting the temperature distribution given in Fig. 9, it is seen that the induced fracture connects to a few isolated fractures in the low-density case, thus it has negligible effect on the heat production. However, the induced fracture provides higher contact surfaces for improved injectivity/productivity of the system. In the medium density case, the induced fracture at the production well connects to several favourable fractures located behind the production well, thus providing access to hot matrix (i.e. larger area of the reservoir) and improving the heat production more than the RJD laterals. In the high-density case, the induced fracture at the production well connects to many fractures around the well, so it improves heat production compared to the well itself, but less efficient than the RJD laterals. Thus, the hydraulic fracturing of the RJDs in the high fracture density case perhaps is not effective as it neither improves the injectivity/productivity nor the heat production.
While the fracture apertures were constant in previous simulations, another set of simulations with deformable matrix in which the fracture apertures follow the Barton-Bandis model given in Eq. (7), are performed. The results for the pressure, pressure difference (Dp) and production temperature for deformable matrix are also included in Figs. 6c, f and 7c, respectively. Better injectivity is observed in the Fig. 6c and f for the deformable case, compared to the constant aperture cases (rigid fracture). This is due to the increase in the fracture apertures due to reduction in the contact stress around the injection well. The reduction in the contact stress comes from both poroelastic and thermoelastic deformation of the matrix and fracture. The pressure build-up at the injection well reduces the contact stress, and the reduction in the matrix temperature results in the contraction of the matrix and thus further reduction in the contact stress in the fractures. Opposite behaviour is seen at the production well where the decrease in the fluid pressure results in an increase in the contact stress and thus lower fracture aperture and lower conductivity of the fractures. However, the effect of the pressure reduction at the production well is not as significant as the effect of cooling the matrix around the injection well. The positive effect of the deformable matrix at the injection well is not as much as the effect of RJD laterals for the low fracture density case, which due to scattered fractures it is not surprising. However, for other two fracture density cases e medium and high e the effect of the deformable matrix on the injectivity is as much as the effect of RJD laterals. When it comes to the heat production, the deformablilty of matrix has negligible effects on the heat production from low fracture density case but has negative impact on the heat production from the medium to high density cases. This is due to unfavourable flow paths created in the deformable case. Therefore, for medium to high fracture densities, the effect of matrix deformability and dynamic fracture aperture variation should be considered to achieve realistic predictions for the geothermal system.
In the cases described above, the matrix had a relatively high permeability (k m ¼ 1 Â 10 À13 m 2 ) thus the heat propagated in both fractures and matrices. Fractures can play a more significant role in transferring fluids if the matrix does not have a comparable permeability. In another case of high-density fracture with well placement 3 (wells are connected to the fractures), the matrix permeability is reduced to k m ¼ 1 Â 10 À15 m 2 . As a result, the injected clod fluid mainly remains in the fractures as shown in Fig. 10, and the heat energy stored in the matrix remains unmined which results in very fast drawdown of temperature in the produced water. As shown in this figure, the cold front has penetrated towards the production well through the fractures.

Long fractures
For the case of long, connected fractures, the matrix permeability is set to k m ¼ 1 Â 10 À14 m 2 and, both injection and production wells are located in the matrix and away from the fracture network. Four horizontal RJD laterals are added to each well at the centre of the layer. The effect of the RJD laterals' connectivity to the fracture network is investigated using several connected and unconnected cases as shown in Fig. 11. In this figure, the pressure distribution after the first step (1 day) for the long fracture case for three configurations of the RJD laterals: a) 100 m laterals connected to the fracture network, b) 50 m laterals unconnected to the fracture network, and c) 100 m laterals rotated 45⁰ and unconnected to the fracture network are presented. In this figure, the wells and RJD laterals are represented by the high-pressure gradients. Results show that in the first case, the connection of the laterals to the fracture network significantly improves injectivity and productivity as the pressure build-up at the injection well is just 3.81 MPa, and the pressure drawdown at the production well is just 5.45 MPa. The zones with high-pressure concentration around the injection and production wells are extended to the nearby fractures. While in the other two cases with no direct connectivity of the RJD laterals to the fracture network, the pressure build-up and drawdown varies between 13.1 and 23.0 MPa, and 14.0e23.0 MPa, respectively. Interestingly, when the laterals are connected to the fracture network, the higher fracture density around the injection well improves the injectivity more than the productivity at the production well, thus the pressure build-up (3.81 MPa) is considerably less than the pressure drawdown (5.45 MPa). But when the laterals are not connected to the fracture network, the density of the fracture network does not affect the pressure build-up and drawdown, as the overall conductivity is mainly determined by the matrix permeability.
The injection/production pressure and pressure difference between the injection and production wells (Dp) during 30 years of simulation for all cases of low permeability matrix with fracture network are shown in Fig. 12. Tremendously high-pressure buildup and drawdown is observed for the case without any laterals as well as for the cases with laterals unconnected to the fracture Table 1 The rock and fluid properties used in the simulations.

Parameter
Value Unit Matrix porosity (f) 0.25 e Matrix permeability (k m ) 1 0 À15 ,10 À14 , 10 À13 network, that makes the heat production impractical. While for the case with four laterals of length 100 m which connect the wells to the fracture network, the build-up and drawdown pressures are finite. Thus, in low permeability reservoirs it is crucial that the RJD laterals connect the well to the fracture network. The temperature drawdown profile for different cases of long, connected fractures is shown in Fig. 13. When the wells are connected to the fracture network by the RJD laterals, both injectivity/productivity and the heat production improves. Interestingly, for two cases with unconnected laterals (50 m laterals and 100 m rotated laterals) the heat production curves follow the same trend, which indicates that the matrix permeability controls the flow, thus laterals' lengths do not have significant effect on the direction of cold water movement. Included in this figure is the case with connected RJD laterals in deformable matrix. In this case, the fracture aperture varies during simulations due to the Barton-Bandis model presented in Eq. (7). Cooling down the matrix during heat production reduces the in situ stresses and increases the fracture apertures as shown in Fig. 14.
Around the injection well, the maximum aperture has increased to more than 2 mm from the initial value of 1 mm, while the aperture has reduced to 0.78 mm around the production well. Higher apertures result in higher conductivity of the fractures and higher injectivity, but at the same time creates unfavourable flow pathways towards the production well. Thus, the heat production reduces when the deformability of the matrix is considered, as shown in Fig. 13. The temperature distributions after 30 years of simulation for three cases: without RJD laterals, with 100 m laterals connected to the fracture network, and with 100 m laterals connected to the fracture network in deformable matrix are shown in Fig. 15. The laterals at the injection well divert the injected cold fluid towards areas behind the well and away from the production well. Those areas are marked by circles in Fig. 15b. Thus, laterals create access to more heat stored in the matrix and improve the heat production as shown in Fig. 13. In the deformable case, on the other hand, the higher fracture conductivity enhances the flow of the cold fluid towards the production well, leaving the heat stored in the matrix untouched in areas shown by circles in Fig. 15c.

Net energy production
The net energy production of a system is crucial as it estimates the performance of the geothermal reservoir. The net energy is determined by the energy produced (E g ) minus the energy consumed by pumps for injection (E p ) The energy produced from the production well may be written as (Vik et al. [54]) where Q is the production flow rate, T inj is the fluid temperature at injection. LT is the lifetime for when the temperature at  production reaches the minimum temperature the plant sees it economical to produce. The minimum temperature is set to 65 C in this study. The production energy is basically the area between the production temperature and the injection temperature of 20 C for the lifetime of the reservoir (until the production temperature reaches to 65 C). Energy used to pump in the injection fluid may be written as [50].
where DP is the pressure change between the injector and producer wells, and έ is the energy conversion efficiency factor, that is assumed to be 0.7 [50]. A sensitivity analysis on the length of the RJD laterals is performed to investigate its effect on the net energy production. The length of the RJD laterals is varied between zero (no RJD laterals installed), 25 m, 50 m and 100 m. Four RJD laterals are placed in the centre of each injection and production wells, for three different well placements and in three cases of short fractures: low, medium and high fracture densities. The results for the average net energy rate (net energy divided by the lifetime of the reservoir) versus laterals length are shown in Fig. 16. The average net energy is calculated as the sum of the net energy produced (positive values) divided by the time of that positive net energy production. The negative values for the net energy production are omitted. For all cases, the average net energy production increases with the length of the laterals. The reasons for this are two-fold: Firstly, the longer laterals improve injectivity/productivity of the  well hence reduce the energy needed for pumping. Secondly, the longer laterals provide access to the fracture network around the producer well, and hence provide access to the hot matrix behind the production well. For well placement 1 (in the matrix and away from the fractures), the net energy rate is more sensitive to the lateral length, as the well is placed in the matrix and away from the fractures. Therefore, higher length of laterals improves the injectivity/productivity of the wells. Higher densities of the fractures also improve the net energy rate produced as the fractures improve the injectivity/productivity of the system. For well placement 2, the net energy rate improves faster in short laterals (i.e. 25 m and 50 m) than the long laterals (100 m) as the short laterals perhaps already connect the injection well to the fractures, thus the longer laterals have lower impact on the average net energy production rate. For well placement 3, the wells are already connected to the fractures, thus the RJD laterals have minor impact on the injectivity/productivity of the wells as shown in Fig. 16. Interestingly, the 100 m laterals remove the dependency of the average net energy production rate on the well placement, and all the three well placements for each case of the fracture densities (low, medium, and high) almost develop the same average net energy production rate. The well placement 1 in high fracture density case is an exception as the long laterals improve the heat production by connecting the production well to the hot matrix behind the production well.

Conclusions
The application of the RJD laterals in enhancing the heat production from the geothermal reservoirs was investigated using several cases of low, medium and high density short fractures as well as long connected fractures with connected/unconnected RJD laterals. The key findings of this study are as follows: The RJD laterals can be very effective in enhancing injectivity/ productivity and/or heat production from the fractured reservoirs. The best performance of the RJD laterals in enhancing injectivity/productivity was observed in low fracture density cases in which the main wellbore was not directly connected to the fractures. In the high fracture density case, the RJD laterals   can be effective in heat production by providing access to the hot rock. In low permeability reservoirs that the flow mainly occurs through the fracture network, the RJD laterals should connect to the fracture network in order to be effective. The RJD laterals can be further used to induce hydraulic fractures in low fracture density cases, while in medium to high fracture density cases, the deformability of the rock matrix and its effect on the fracture aperture variation should be accounted for. The deformable matrix improved the chance of creating unfavourable pathways between the injection and production wells. Results for the net energy production rate showed higher sensitivity to the lateral lengths in low fracture density cases and in cases where the well was not connected to the fractures. The laterals can in principle remove the dependency of the net energy rate to the position of the wells with respect to the fractures by creating direct connections between the wells and the fractures.
In summary, the RJD technique can be used as an alternative well stimulation technique for geothermal reservoirs provided that the formation can be jetted properly. The extend, direction of laterals and their connectivity to the natural fractures are crucial parameters in evaluating the success of this technique.

Acknowledgments
Authors thank the European Union for partially funding this work through the EU Horizon 2020 research and innovation programme, under grant agreement No 654662. The research leading to these results has also received funding from the Danish Hydrocarbon Research and Technology Centre under the Advanced Water Flooding program.