Numerical simulation of a Deep Borehole Heat Exchanger in the Krafla geothermal system

Abstract The geothermal energy sector is facing numerous challenges related to heat recovery efficiency and economic feasibility. Research on superheated/supercritical geothermal systems is progressing in Europe, triggered by the Iceland Deep Drilling project (IDDP) and the DESCRAMBLE project in Italy. In Iceland, the IDDP-1 well, which reached a magma intrusion at a depth of 2100 m, raised new opportunities to untap the geothermal potential near magmatic intrusions. Given their highly corrosive nature, geothermal fluids weaken the wellbores integrity during conventional geothermal production. Closed-loop Deep Borehole Heat Exchangers (DBHE) that do not require fluid exchange between the subsurface and the wells represent a strategic alternative for recovering heat from these unconventional geothermal resources, while minimising the risk of in situ reservoir damage. The thermal influence and heat recovery associated with a hypothetical DBHE drilled into the IDDP geological settings are investigated via Computational Fluid Dynamics (CFD) techniques, simulating 30 years of production. Two wellbore designs are modelled, based on simplified geological properties from the IDDP-1 well description. The results show that, during the first year of production, the output temperature is function of the working fluid velocity before reaching pseudo-steady state conditions. The cooling perturbation near the bottom hole is shown to grow radially from 10 to 40 m between 1 and 10 years of production, and the calculated output power reaches up to 1.2 MWth for a single well. The heat transfer at the bottom well bore is enhanced by extending the inner well deeper into the ground. Subject to full economic analysis to be performed at field scale, the significantly lower technical risks of the closed-loop DBHE could outweigh the lower thermal output per well compared to theoretical expectations from open-loop Enhanced Geothermal Systems (EGS).


Introduction
The need for sustainable energy supplies raises new challenges to improve low-carbon emitting technologies and decrease dependency on fossil fuels [1]. Geothermal energy already contributes to generating power and heat worldwide, having reached a total energy of 73.689 GWh in 2015. It is expected that this value will increase in the future [2,3]. In Europe, the installed geothermal electricity capacity in 2017 was about 2.8 GWe, generated from 117 power plants [4].
While Enhanced Geothermal Systems (EGS) face technological issues due their potential to induce or trigger seismicity [5,6], closed-loop single-well solutions can prevent fracture clogging and fluid losses in porous reservoirs. Borehole Heat Exchangers (BHE) or Deep Borehole Heat Exchangers (DBHE) aim to extract geothermal energy by circulating a working fluid in a well without producing geofluids or requesting injection processes, or result in thermal short-circulation between injection and production wells. Due to highly corrosive fluids at high thermal and pressure conditions, unexpected facilities failure may happen. Despite suggesting anticipated low efficiencies compared to the expectations from open-loop EGS, the risk-cost balance analysis coupled with a sustainability assessment could encourage a development of this technology in favourable locations [7,8]. Closed-well completion designs such as the DBHE can bypass the aggressive upflow and low permeability subsurface flow in very hot geothermal reservoirs by targeting exclusively the heat available. However, to date, their efficiency has not been investigated thoroughly.
Studies performed for coaxial BHE up to 1000 m deep at low geothermal gradient (1.0-3 K/100 m) demonstrate that the thermal load is proportional to the thermal gradient [9]. The heat extraction is maximum downhole, resulting in a minimal thermal sensitivity to neighbouring BHEs. Furthermore, for low temperature cases, only 1-2% of the produced energy is dedicated to circulating the working fluid and can be managed with different borehole diameters [9]. The application of DBHE worldwide would imply drilling to a maximum depth-use of 3 km for generating 0.15-2.5 MW th thermal power and 0.25-364 MW e electrical power [7,10]. Several DBHEs exist in Switzerland and Germany, up to a depth of 2.7 km for heat pump systems [11,12]. The DBHE in Weissbad (Switzerland) is 1213 m deep, see Fig. 1 and reaches a temperature T ¼ 45°C.
Assuming ideal numerical conditions, the Weissbad closed-loop well model shows a quick decrease of temperature which enables the system to be used unfavourably with loading and discharging cycles for constant energy production [13].
Apart from the injected fluid temperatures, the flow rate and the thermal gradient, the insulation of the inner pipe appears to be the main parameter ensuring an efficient thermal recovery, whereas the well diameter affects the residence time of the working fluid. Various technology derived concepts, notably horizontal BHEs, provide long-term heat power in the range 0.350-2 MW th . For this technology, three production temperature behaviours can be described as function of the length of the horizontal section [14,15]. DBHE systems were considered to retrofit abandoned oil wells [16][17][18][19][20][21] with various working fluids such as isobutane at supercritical state [22] or water [17,23]. In terms of thermal gradient and depth, several cases from 25°C/km and 50°C/km and between 2-6 km were investigated, giving variable potential output energy. When considering 50°C/km the output power starts to be efficient at depths higher than 3 km. Furthermore, direct power generation systems appear more efficient with supercritical fluids than flashing power generation systems [24]. Research on superheated/supercritical geothermal systems is progressing in Europe, notably triggered by the Iceland Deep Drilling project (IDDP-1) in Krafla [25] and the DESCRAMBLE project [26] in Italy. The IDDP-1 well raised new opportunities to untap the geothermal potential near magmatic intrusion. The well was completed at a depth of 2072 m before being terminated after a rhyolitic magma was encountered somewhere between 2092 m and 2104 m [27].
General models for volcanic geothermal systems suggest a deep circulation of fluid above and on the sides of a magmatic heat source, with an underlying thin conductive boundary layer [28]. A 30-50 m thick conductive layer was estimated above the Krafla magmatic intrusion [29].
Many geothermal numerical studies reported in the literature were performed with TOUGH2 [30,31]. However, this code, although highly efficient, does not provide detailed information on the wellbore flow and surroundings for near supercritical conditions.
To date, only few Computational Fluid Dynamics (CFD) studies have targeted geothermal energy assessments [19,[32][33][34][35]. CFD should however also be applied to predict and study the supercritical conditions present in very hot geothermal systems. CFD-based numerical techniques can provide an integrated and coupled well bore-reservoir heat transfer analysis tool for accurate geometrical designs and real working fluid compositions. As the fluid flow calculation over the whole length of a BHE is dependent on the bottom hole conditions where turbulence is present, it is expected that a CFD simulation will solve the heat and fluid flow equations more accurately than by using traditional 1D numerical or analytical models. DBHEs are conventionally investigated in low-medium temperature areas. A thorough search of the relevant literature has yielded only one related experimental study performed in Hawaii [36], where high temperature conditions of 110°C at the depth of 876.5 m have been measured.

Abbreviations
The implementation of geothermal single-well technology in high thermal gradient settings has so far been the subject of only limited analysis. The energy production from a base case singlewell design (the Weissbad BHE, described in Fig. 1) has been estimated considering the Krafla geothermal system and using field data extracted from the literature.

Material and boundary conditions
Several numerical and analytical simulations consider 1D calculations of the borehole fluid flow and heat transfer [9,37], or full 3D CFD simulations [19,32,34]. The work presented here considers an axisymmetric approach of the well bore, reducing the 3D description to a 2D model, thus significantly reducing CPU and computing time.
The modelled DBHE uses the settings of the Weissbad BHE extended to the depth of 2070 m (see Fig. 1) by keeping the proportion between the lengths of the casing sections. The cement and casing properties are also both extracted from the Weissbad BHE [13]. The flow is directed downward in the external annulus of the well and the temperature increases with depth. The circulating heated fluid then flows from the bottom of the well to the surface via the internal pipe.
Note that there is no contact between the geothermal fluids (or rocks) and the circulating fluid, as the system is fully independent. The circulating fluid chosen in the present study is liquid water, see Table 1. The choice of constant water properties is a limitation of the current model. However, this study is intended to be the first step into modelling the complex downhole geothermal environment with a CFD approach. Note that constant water properties are commonly used for similar problems, see [19,20,32] for instance. As the heat exchange surface is small (well walls) compared to conventional geothermal extraction (fluids flowing into pores and fractures), the circulation flow rate is also small to allow a sufficient residence heating time.
The computational model is represented by an axisymmetric domain extending 200 m radially (similarly to [21]) to avoid any influence of the boundary temperature on the circulating fluid. However, as it is shown later in this work (c.f. Figs. 16 and 17), any lateral distance greater than 85 m would have been suitable to prevent any influence of the lateral boundary condition on the flow inside the well.
From a simplified conceptual model [29], the Krafla volcanic system located in the vicinity of the IDDP-1 well can be described by considering three distinct zones, as shown in Fig. 2, of properties summarised in Table 1. The cylindrical area outside the DBHE is made of basalt with constant materials parameters (q = 2700 kg=m 3 ; c = 800 J/kg.K; k = 2.5 W/m K) [29]. A 30-50 m thick conductive layer (Zone 3) separates the porous Zone 2 (similar to the reservoir) from the magma region located at a depth of 2140 m. This insulated layer is set with a lower thermal conductivity value (k = 1.5 W/m K) corresponding to a solid rhyolite [29].
The wall of the well is defined with a constant 7 mm thickness whose thermal conductivity and specific heat vary with depth. These variations correspond to the different casing and cement thickness sections identified in the Weissbad BHE ( Fig. 1).
Note that in reality, the IDDP-1 well is completed at the bottom hole with a 9-5/8 00 K55 casing of thickness 13.843 mm [27]. The 7 mm wall of the Weissbad DBHE might thus not be sufficiently thick, considering the harsh environment of the Krafla well.
The CFD code solves a 1D steady heat conduction equation at the wall to compute the thermal resistance, Dx/k, where k is the effective conductivity of the wall material and Dx is the wall thickness [39]. The effective thermal conductivity of the external wall is adjusted with respect to the depth (see Table 1). The study focuses on forced heat convection in the well, simplifying the conduction heat transfer through casing and cement sections. The water enters the DBHE with an initial velocity and constant temperature and the outlet of the well is set with a pressure outlet boundary condition. A pressure inlet atmospheric conditions of 1 bar and 10°C is applied at the surface of the porous region. The bottom boundary of the model is set with a constant temperature of value discussed in Section 5.1. The sides of the porous model are considered adiabatic. A roughness height value of 0.00026 mm [24] is applied at the walls. Two different designs are investigated, based on the depth of the inner pipe. Design 1 considers the same proportions Table 1 Materials properties applied for the fluid and zones referred in Fig. 2.

Mathematical approach
An extensive review of methods for modelling a wellbore heat exchanger has recently been published [7]. Specific geothermal numerical codes such as the TOUGH2 code [30], AUTOUGH2 [40], the Complex System Modelling Platform CSMP++ [41] can be used to assess the porous heat transfer and fluid flow. 1D semi analytical calculations are often considered to be an efficient way of simulating reservoir-wellbore systems, as it is the case for T2Well [42]. CFD tools are mainly used to investigate the near-wellbore processes [19,32,34,43,44], where detailed analysis is required.
The CFD code ANSYS Fluent 17.1.0, which is based on the finite volume approximation, is used in this study to solve the governing equations of fluid. The porous medium model considers a superficial velocity porous formulation based on the volumetric flow rate. Ignoring the convective acceleration and diffusion, the pressure is proportional to the fluid velocity, reducing to the Darcy law [45].
The general form of the continuity equation in an axisymmetric geometry can be written as: In addition, the axis and radial momentum in a 2D axisymmetric domain read: The energy equation is: where J j ! is the diffusion flux of species j and s eff is the stress tensor [45]. In the porous medium, the solid and fluid region is calculated relatively to the porosity value [45]. For incompressible flows, the energy term is calculated as: where c is the heat capacity and v is the fluid velocity. The details behind Eq. (7) are available in the ANSYS Fluent Users Guide [45].
The realizable k-e turbulence model has been selected; this model modifies the calculation of the turbulent viscosity l t and the dissipation rate e [45].

Numerical model
The PISO (Pressure-Implicit with Splitting of Operators) algorithm is used to couple velocity and pressure; this model is recommended for transient flow calculations [39]. The PRESTO (PREssure STaggering Option) scheme is considered for the pressure discretization; this model is suitable for porous-based flow simulations [39]. The second order discretization scheme is applied for all convection terms. The residual convergence criteria for all equations is set to 10 À6 . The time-step is increased as necessary from 1s to 2 h, and 30 years are simulated.

Organic Rankine Cycle analysis
Several geothermal conversion technologies are applied worldwide, depending on the resource characteristics [46]. Despite a high temperature gradient investigated in this study, the electrical power generated is assumed to be produced under long-term conditions with water in the DBHE remaining as pure liquid, below the saturation point. Thus, the Organic Rankine cycle (ORC) is selected to convert the thermal energy into electricity. ORC belongs to the binary cycle geothermal power plants that use a second circulation system with a second working fluid having lower evaporating conditions than the geothermal fluids. In a heat exchanger, the latter is vaporized and directed to the steam turbine. Fig. 3 describes an ORC power plant used to convert the thermal energy extracted by the DBHE into electrical power. For the energy conversion calculation part, water entering the power plant is assumed to be a saturated liquid. Water properties are calculated with freesteam 2.1, an open source implementation of the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam.
In an ORC cycle, the heated water flows through the ORC power plant heat exchanger. This water needs to be kept at a pressure above its flash point in order to avoid the breakout of the steam [46]. The preheater (PH) brings the ORC working fluid up to the boiling point before reaching the condition of saturated vapour in the evaporator (E). The method used to obtain the thermal and electrical power of the DBHE system is extracted from [20]. The pressure-enthalpy diagram is presented in Fig. 4. The Massachusetts Institute of Technology (MIT) correlation [47] is used based on the worldwide active ORC plants to compute the thermal efficiency v th relative to the temperature extracted from the DBHE.
Following the work described in [20], the organic fluid R-C318 is chosen, with a Ph/E ratio of 47.2.
v th ¼ 0:0935T À 2:3266 The thermal H and electrical power E are calculated with:

CFD model sensitivity
The mesh, composed of hexahedral cells with local refinements from 0.1 mm to 1 m, is generated with ICEM CFD 17.1.0, and exported into the flow solver ANSYS Fluent. A mesh independence study is performed to identify a suitable mesh. Fig. 5 shows the numerical output temperatures and the radial temperature distribution in the well below the internal pipe (2008.1 m), obtained with three different meshes.
As can be seen in Fig. 5 (left), simulated output temperatures are nearly identical for all three meshes and it is difficult to conclude on the mesh independency when looking at those results only, even if a temperature difference is present initially for the 313,200 cells-based result. The radial temperature distribution obtained after 95 h at a depth of 2008.1 m, allows for more discussions. Differences are now clearly noticeable between the results obtained with the meshes comprising 313,200 and 418,000 cells, while there is no significant difference between the results with the meshes comprising 418,000 and 764,720 cells. In the interest of computational time, the coarser of these two meshes has therefore been selected for the rest of the study. Note that a mesh independence study has also been performed for Design 2, and a suitable mesh of 424,800 cells has been selected, as reported in Table 1.
The time-step is progressively increased from 1s to 2 h, to control the stability of the numerical solution and reduce the overall computing time. Fig. 6 highlights the influence of the number of iterations per time-step, which modifies the results at the beginning of the transient simulations. It appears clearly that the default number of iterations (20 for Fluent) is not high enough to obtain reliable results as significant differences can be noticed when increasing this number. The influence of the number of iterations is however reduced when applying between 60 and 100 iterations per time-step; 80 iterations per time-step are therefore applied in this work. Fig. 6 also shows that results converge and a very slow temperature decrease is present after 60 days.

Initial conditions of the Krafla subsurface model
Initial conditions near the bottom of the IDDP-1 well are discussed before evaluating the energy production in the DBHE.

Surface down to À2 km
The numerical model is simplified compared to the natural heterogeneity and complex processes above the magmatic intrusion. The initial physical conditions at the bottom of the DBHE are investigated to prevent overestimating the thermal energy available to the DBHE.
The bottomhole conditions in the IDDP-1 well were previously studied numerically [29,38,48], as temperature log measurements did not reach the bottom hole depth. Following two independent methods, the estimates appeared to be in the range 390-400°C in the first study and around 450°C in the second one [38]. The discharge test showed a transient temperature of around 450°C [49], suggesting the existence of potentially higher temperatures due to the recharge water circulation accessing the thermal energy underneath [29].
Starting from the surface down to 2 km, the temperature gradient extracted from the data log measurement is implemented in the CFD flow solver via an user-defined function (UDF). Below 2 km, pressure and temperature distributions are calculated under steady-state conditions, see Figs. 7 and 8.
It should be noted that the Krafla geothermal fluid is more likely to be a H 2 O-NaCl mixture, as investigated with CSMP++ [50]. The developed model described here uses constant water properties,   which overestimates the pressure distribution below 1000 m, where phase change and thermal pressure dependent properties affect the fluid pressure distribution.

Around the bottom hole
Numerical simulations described in the literature apply a constant magmatic temperature in the range 850-900°C [29,51,52]. Regarding the metamorphism above the intrusion, the minimal heat flow from the conductive boundary layer reaches 23 W=m 2 [53]. Fig. 9 shows the sensitivity of the initial steady-state temperature and pressure based on three different magmatic temperatures and the thickness of the conductive layer (30-50 m).
At 2 km, pressure and temperature are extracted from Fig. 7.
If the intrusion temperature exceeds 750°C, a minimum steady temperature of 520°C is obtained. The only case below 500°C is obtained with a 50 m conductive layer thickness above an intrusion at a temperature of 650°C, a simulated magmatic temperature below the published estimates at 2.1 km [29,53]. Note that the thermal conductivity of the insulated layer could be lowered, as studies invoke the existence of very low conductivity minerals, such as felsite at the edge of the magma layer in the Krafla geothermal reservoirs [53]. However, for the sake of simplicity, the configuration including a 50 m conductive layer above an intrusion at a temperature of 650°C is selected, to keep a relevant thermal distribution near the bottom hole.
The potentially recoverable heat from the simplified IDDP-1 bottom hole surroundings will now be discussed, based on these assumptions.

Results and discussions
6.1. Velocity Table 2 summarises values of downward velocity in the bottom section of the DBHE (V down ), upward velocity (V up ), mass flow rate ( _ m), produced temperature (T prod ), thermal power (W Th ) and electric power (W e ) generated after 30 years of operation. Fig. 10 shows the production water temperature over time for different injection velocities. Note that all results presented until Section 6.4 are obtained considering a 10°C injection temperature. The first weeks of recovery show a very high temperature, in agreement with the measurements in the test well [49]. During this period, pressure needs to be sufficiently high to avoid phase change in the well.  A strong temperature decrease is observed in the early stages, as the velocity is high before a transitional temperature decrease appears and reaches a pseudo-equilibrium around 10 years.
The velocity strongly constrains the heat recovery behaviour which describes three stages as detailed for a horizontal closed system [15]. The temperature distribution in the injecting annulus and internal well is shown after 10 years of simulation in Fig. 11.
The slower the water flows down, the more heat it recovers, but also the more thermal losses it faces along the inner pipe.
Water does not circulate till the bottom hole if the internal well is not deep enough, as shown in Fig. 11. All along the producing and injecting sections, there is a 2.5-3.20% long-term increased temperature for Design 2 compared to Design 1. Below the internal well (À2008.1 m and À2069 m), the temperature rises due to equilibrium with the reservoir temperature. Pumps are required to flow water inside the wells.
As reported in Table 3, the pressure differences are nearly identical for Design 1 and Design 2, with values for Design 2 being however slightly higher when increasing the inlet velocity. As expected, the pressure difference is also higher if wall roughness (r h ¼ 0:00026 m) is considered due to frictions affecting the energy consumption for pumping the fluid in the DBHE. However, no further investigation is performed as only minor differences are seen in the temperature long-term production when roughness is, or is not, accounted for.

Heat transfer analysis
The transient horizontal heat flux through the vertical walls of the well is analysed for 0.025, 0.1 and 0.5 m=s. A negative value stands for a gain of heat, whereas a positive value means a loss of heat to the neighbouring area. Fig. 12 shows the horizontal wall heat transfer for the internal wall at 1 and 10 years.
Values under 10 À3 W/m 2 are filtered out. Despite a low thermal conductivity value to simulate the thermal insulation of the inner pipe, heat is added to the annulus part. The heat flux values are higher when the injection velocity is lower. The vertical heat flux on the external well for the two designs is shown on Figs. 12-14. Up to a depth of 2 km, the heat flux values are similar for both designs. The large range of values observed at À1366.2 m is due to the sudden change of thermal conductivity between the vertical layers of the casing, as shown in Fig. 1. The heat transfer is stronger when the velocity is high but decreases slowly in time between 1 and 10 years.
Descending along the wellbore, the heat transfer value increases as the injection velocity increases. Below the internal well (À2069 m), the Design 2 configuration shows horizontal heat flux values in the range 3.5-5 MW=m 2 , stronger than for Design 1 (2.5-3.8 MW=m 2 ). Zero heat flux values are obtained right below Table 2 Vertical velocities, mass flow rates and temperature of produced water in the DBHE for an injection temperature of 10°C, thermal and electrical powers estimated from an ORC power plant after 30 years.    Table 3 Pressure difference (DP) between injection and production well at the surface (bars). À2009 m (Design 1) or À2070 m (Design 2), suggesting a thermal equilibrium between the well bore and the geothermal reservoir.

Long term thermal influence
At the bottom of the two designed wells, significant vertical heat flux is present, as shown in Fig. 15. The bottom heat transfer values for Design 1 are low (1.76-1.78 MW=m 2 ), and appear only slightly sensitive to the injection velocity. However, Design 2 shows high vertical heat fluxes associated with a high injection velocity, as previously described for the horizontal heat transfer (Figs. 13 and 14). Fig. 16 shows the radial temperature at the bottom hole position in the reservoir after 10 years, following the injection velocity. It can be noted that Design 2 only cools down the bottom hole radial environment. This cooling effect is higher if the injection velocity is high, e.g. associated with a higher heat flux. After 10 years, a steady reservoir temperature is reached at a distance  Fig. 18 shows the thermal and electric powers from an ORC plant calculated using Eqs. (9) and (10), after 30 years of operation. The thermal power calculated reaches up to 1 MWth considering a 0.5 m=s injection velocity and an injection temperature of 10 or 20°C. The overall thermal power production is lower for high injection temperatures. On the contrary, the high injection velocities and high injection temperatures lower the electric power values. The best case scenario for production electricity using an ORC cycle is for an injection velocity of 0.05 m=s.     The velocity optimization has shown that Design 2 generates more heat. However, it also generates a deeper cooling perturbation than Design 1, whereas the improved heat production is only approximately in a range of 2-3%. The choice of an optimized design has therefore to be balanced between more thermal power generated and a sustainable induced cooling effect.

Comparison of the method with other published results
Due to the lack of experimental data, the numerical results from this study could not be validated. The initial temperature and pressure gradient along the wellbore are the only robust data available here. A comparison with other simulated results from the literature for closed-loop systems is presented in Table 5. It can be seen that the numerical results from the present study agree with other published data. The present study focuses on more favourable conditions than previously investigated in the literature. The assessment performed for an output power of 134 kW e is given for a short-time calculation. A 45% decrease of efficiency is observed after five years, decreasing the power to 675 kW th and 60 kW e [20]. It can be concluded that the overall performances of a DBHE near a magmatic chamber are improved.

Discussions
Several initial and geometrical assumptions are made in this study. The homogeneous geometry and characteristics of the reservoir considered cannot mimic the real complexity of the zone, as described in [55]. The bottom temperature chosen in the model is lower than the temperature estimates of the rhyolite magma, around 850-900°C. This suggests applying a lower thermal conductivity value than used in the current work for the discontinuity (ductile felsite [55]), to better approach the underlying heat transfer.     Table 5 Numerical data from the literature [13] and from this study. The single-phase formulation for water with constant properties might produce an underestimation of the heat transfer for the reservoir conditions. The geofluids probably also have higher thermal properties than considered here, notably higher values of heat capacity due to supercritical/superheated conditions [56].
In addition, the working fluid at the bottom hole section shows temperatures above the critical point of pure water (373.94°C), suggesting a supercritical behaviour if pressure is maintained high enough. Although this study is intended to be the first step into developing a better understanding of and better characterising the physical processes taking place in unconventional geothermal wells, further investigation needs to be considered to numerically model the BHE under such high conditions with varying density and other properties [20].
Such fluid properties temperature dependency can generate a thermo-siphon-effect, which compensates the pressure losses due to density variations of the ascending and injecting working fluid, as shown in [20]. In some cases, pressure at the outlet can become greater than injected at the inlet of the system [20]. Therefore, the implementation of natural lateral geothermal fluids circulations from the reservoir recharge zones can be investigated to trigger potential large scale temperature modifications.
In the calculations of the overall output power, the pumping power consumption for compensating pressure losses are not detailed. The magmatic object encountered in Krafla is assumed infinitely supplied and the cooling extension of the magma and changing permeability/porosity of the rock in contact of the intrusion [29] are not considered. The use of a DBHE can also be considered with pulsating injections, leading to resting times for the well to be heated up.

Conclusions
Deep, open-loop geothermal doublets of the Enhanced Geothermal System type, for electricity generation are site-dependent and technology is still in the demonstration stage. DBHEs have the potential to bypass the hurdles of open-loop systems; however, to date, only a handful of DBHEs have been implemented worldwide, with mixed success and primarily for heating/cooling applications. Bringing the DBHE concept to particularly favourable geological settings such as those considered in this paper could untap significant geothermal capacity.
A geothermal BHE CFD model close to a magma intrusion has been investigated. Transient simulations have been performed, based on initial conditions from the IDDP-1 well, drilled in Iceland, considering two designs from the BHE configuration in Weissbad, Switzerland. A turbulence model has been applied in the well with a porous fluid flow formulation in the radial direction extended to 200 m. Conclusions have been drawn from this numerical study: Considering non constant well bore temperatures, the outlet temperatures from the DBHE have been constrained by the injection velocity and well design parameters. The DBHE has produced a very early high temperature, but the heat has declined rapidly to reach a pseudo-steady state for a designed operating time of 30 years. It has been shown that the use of a high insulating material in the upper part of the inner pipe could reduce the heat losses from the rising working fluid. CFD has been able to identify flow patterns in the well, which can lead to optimized geothermal well bores, based on heat transfer calculations. By deepening the internal well, the DBHE can access more heat, showing 2-3% improvement compared to the standard design. However, it also cools down laterally the reservoir at the bottom hole up to 85 m, after 30 years. The numerical description of the thermal and fluid flow processes remains a challenge; it will be necessary to improve the model in the future to better evaluate the changes of fluid and solid physical properties. In terms of electrical power production, the best case scenario for a single well has been obtained for an injection velocity of 0.025-0.05 m=s, reaching about 90 kW e after 30 years. However, the use of DBHEs, even in high thermal conditions, seems more suitable for direct-use applications. Despite low energy production compared to the short-term expected flow from supercritical/superheated resources, the lower technical risks associated with DBHEs are attractive and can justify a future economic analysis. Increasing the number of wells could efficiently compensate for the limited heat output per well compared to conventional Enhanced Geothermal Systems. Furthermore, completion technique improvements and insulating materials properties might deliver more efficient single-well technologies in the future, for accessing better overall heat extraction.
Although CFD simulations take longer than 1D numerical or analytical models, which can be computationally more efficient, current 1D approaches may carry limitations when applied to modelling unconventional geothermal well designs. More specifically, they may not adequately capture the associated unique and largely unexplored to date fluid flow and heat transfer phenomena within the downhole environment and in its immediate vicinity, e.g. the near-wellbore region. This would be particularly relevant for non-homogeneous, anisotropic near-wellbore regions. The comparison of the two approaches will be the object of future investigations.

Declaration of Competing Interest
The authors declared that there is no conflict of interest.