Conjugated Numerical Approach for Modelling DBHE in High Geothermal Gradient Environments

: Geothermal is a renewable energy source that can be untapped through various subsurface technologies. Closed geothermal well solutions, such as deep geothermal heat exchangers (DBHEs), consist of circulating a working fluid to recover the available heat, with less dependency on the local geological settings than conventional geothermal systems. This paper emphasizes a double numerical method to strengthen the assessment of DBHE performances. A computational fluid dynamics (CFD) commercial software and the 1D coupled wellbore-reservoir geothermal simulator T2Well have been used to investigate the heat transfer and fluid flow in a vertical DBHE in high geothermal gradient environments. The use of constant water properties to investigate the energy produced from DBHEs can lead to underestimating the overall heat transfer at high temperature and low mass flow rate. 2D axisymmetric CFD modelling improves the understanding of the return flow at the bottom of the DBHE, readjusting and better estimating the pressures losses commonly obtained with 1D modelling. This paper highlights the existence of convective cells located at the bottom of the DBHE internal tubing, with no significant effects due to the increase of injected water flow. Both codes are shown to constrain the numerical limitations to access the true potential of geothermal heat extraction from DBHEs in high geothermal gradient environments and demonstrate that they can be used for geothermal energy engineering applications.


Introduction
Geothermal is a key candidate among the renewable and low carbon emission source of energy, thanks to its great potential throughout the world. The Earth surface observed an approximate heat flux of 47 ± 2 TW [1], with a mean thermal gradient of 30 °C/km. Enhanced geothermal systems (EGS) is aiming to potentially double the worldwide power production from geothermal resources [2], but is still technically discussed due to induced seismicity, poor rock connectivity, and low recovery factors [3]. In this respect, closed-loop unconventional geothermal well designs are able to overcome these issues [4][5][6], with economic benefits and CO2 emission savings [7], notably from very deep drilling into favorable higher geothermal gradient (>45 °C/km).
The numerical modelling of such systems is a key step to establish the benefits of using complex designs for future implementation in testing facilities. Estimating production scenarios and the expected energy outcome can help decision-making.
Various closed-loop well designs modelling techniques have been reviewed in [8] for various geometries (vertical, inclined, horizontal, and u-tube) and working fluids (water, organic fluids, and CO2). A U-shaped, closed-loop concept has been assessed in [9], circulating CO2 at a flow rate of 25 kg/s, leading to an estimated recovery of over 2.5 MWth via a thermosiphon, without pumping. The hot fluid is self-flowing to the surface due to density changes.
Vertical deep borehole heat exchangers (DBHEs) have been thoroughly investigated in mediumdepth temperature ranges with various modelling methods [8]. The potential repurposing of hydrocarbon wells has also been assessed through such designs [10,11]. It represents a significant amount of heat available without drilling costs [12]. Figure 1 shows the schematics of an actual DBHE experimented in Hawaii [13], for which water flows downward in the annular space and returns to the surface in an inner tubing. This system reached a depth of 876.5 m, with a temperature of 110 °C at the bottom. The circulating water was observed to reach the surface at 98 °C, with maximum gross and net thermal outputs of 540 kWth and 370 kWth, respectively. While conduction is typically regarded as the main heat transfer mechanism in DBHEs, advective porous transport in volcanic areas close to magma is expected to contribute to additional heat recovery [4]. Similarly, a coaxial open downhole heat exchanger set in China observed a stable heat output of 275 kW [14]. The overall performance of a DBHE is limited by the heat exchange surface offered by the wellbore cylindrical area between the working fluid and the surroundings rocks. Producing energy from vertical DBHEs can be particularly efficient for supplying heat, rather than electricity, even in the case of very high geothermal gradient [15]. Investigations regarding the DBHEs performances available in the literature have mainly assumed a pure water phase, for both analytical studies [16][17][18] with the wellbore modelled as a 1D source line [19], and numerical studies through computational fluid dynamics techniques (CFD) [20].
Analytical DBHE models have shown good performances for an optimization analysis but overestimate the pressures losses at low mass flow rate and underestimate them at higher mass flow rate, compared to numerical models [21]. Numerical methods have been established as the most accurate approach, notably for modelling early transient time processes [21].
Typically, constant water properties have been applied to estimate the heat transfer in DBHEs [10], but constant fluid properties are known to lead to errors in pressure losses calculations [21]. Density variations can have an impact on the thermosiphon-effect, as it compensates pressure losses due to density variations of the ascending and injecting working fluid [22], and surface pressure from the internal tubing can be higher than the pressure at the injection point. While the assumption of constant water properties can be valid in low geothermal gradient systems (1.8% difference in the thermal resistance in [23] for 30 °C/km), this assumption becomes arguable for higher geothermal gradient and deeper DBHEs. Constant water properties tend to underestimate the produced water temperature and the heat flux for a DBHE of 6100 m [22] and overestimate those by 11% in a DBHE of 3500 m [24]. Therefore, the use of non-constant thermophysical properties such as the fluid specific heat, viscosity and thermal conductivity is critical [25].
To aid further developments of innovative closed-loop wellbore designs (two-phase closed wells [26], improved conductive wellbores [27], supercritical fluids [28]), accurate modelling of the heat transfer and fluid flow is a key requirement in the validation pathway.
Three CFD models of DBHE sections at different depths have been created to discuss the heat transfer and fluid flow under high geothermal gradients. Those models are based on the preliminary results from a DBHE system [18] using T2Well, an integrated 1D wellbore-reservoir simulator [29]. This study aims to provide a critical assessment of the 1D wellbore formulation based results compared to those generated with a two-dimensional (2D) axisymmetric CFD model with ANSYS Fluent [30], also highlighting the limits of using constant water properties compared to temperaturepressure dependent water properties (IAPWS-IF97 [31]).

Model Parameters
The numerical study performed in this work is based on the theoretical model described in [18]. T2Well [29] uses a dedicated multiphase transient calculation coupled with the geothermal simulator TOUGH2 [32]. The calibrated T2Well model presented in [18], based on a closed DBHE system in Hawaii [13], is extended to the depth of 1990 m with a bottomhole temperature of 350 °C. The energy flow rate calculated at 1 year in this configuration does not exceed 1.4 MW [18]. Two injections of water at 10 °C are considered, at mass flow rates 2 and 10 kg/s. Table 1 shows the casing properties, the values of the constant water properties used in the models, and the radial dimension of the DBHE. The thermal conductivity of the homogenous geothermal formations is set to 1.6 W/m.K and the specific heat to 870 J/kg.K. A permeability of 10 -16 m 2 and a porosity of 1% are assumed for this study. The calculated values from two vertical adjacent cells in the DBHE section from the T2Well model are extracted and applied as boundary conditions for a 10 m vertical, axisymmetric CFD model. T2Well/EOS1 calculates the values at nodes in the center or an arbitrary position in a given cell (at the boundary between the well and the surroundings for N3 in Figure 2) Constant water properties and the IAPWS-IF97-based formulations for density, thermal conductivity and viscosity are separately applied in the CFD model. However, the specific heat capacity remains constant as it can only be function of temperature [30]. A time period of 1 year is considered, and data are extracted at depths of −500 m, −1900 m, and −1987 m. Figure 2 describes the boundary conditions obtained from T2Well and applied to three CFD models, corresponding to the depth of −500 m, −1900 m, and the bottom of the DBHE. The node located at a depth of −495 m (N1-495) provides the temperature and pressure in the annular section to be applied at the outlet of the CFD model. Similarly, the N1-505 node provides the values for the inlet of the CFD model. A constant temperature boundary is applied at the face of the internal (C1) and external (C2) casing (top and bottom), equal to the values from the nodes N2 and N4. A gradient of temperature is applied onto the side wall of the well (established from the temperature values of the nodes N5-495 and N5-505), corresponding to the contact between the cement and the external casing (C2). The same strategy is considered for nodes located at the depth of −1900 m. At −1978 m, two geothermal gradients (Grad1T −1987/−1989 and Grad2T −1989/−1990) are considered to account for the large increase of temperature in the bottomhole section. A pressure inlet and a constant temperature are applied at the bottom of the internal well and top of the annular tubing, while a pressure outlet condition is applied at the top of the internal well and bottom of the annular tubing, detailed in Table 2.

Mathematical Model of T2well
T2Well is a coupled wellbore-reservoir simulator, solving a 1D two-phase momentum equation for the wellbore and the 3D multiphase Darcy Law in the reservoir [29]. The IAPWS-IF97 formulation [31] is considered for the equation of state of water. The well and the rock formation connection are closed in order to mimic heat transfer only between the DBHE and the surroundings, without fluid flow. Kinetic energy is considered in the wellbore, where water velocity is high, but it is neglected in the porous reservoir, where low fluid velocities are present. The conservation of mass (k = 1) and energy (k = 2) for a single phase system can be written as where V is the volume (m 3 ), n is the outward normal vector, is the surface area of the well (m 2 ) and q is the source or sink terms for the mass or energy component. M 1 is the mass accumulation term and F 1 is the mass flux, which are defined as where ρ is the density (kg/m 3 ), and u is the fluid velocity (m/s). In the wellbore, the energy flux accumulation is The term A is the cross sectional area of the wellbore, θ is the angle between the vertical direction z and the wellbore section, and z is the elevation in the well.
where U is the internal energy of water liquid. The momentum equation in the wellbore is with the perimeter of the cross section, the apparent friction defined as = 16/Re when Re < 2400, and as follows when Re > 2400: where is the roughness of the well walls and Re is the Reynolds number, defined as Re = d/ with d the diameter of the well or, corrected for an annulus via the hydraulic diameter (do − di), and μ the water viscosity. The formulations for the energy flux, accumulation and fluid velocity in the porous reservoir can be found in [29]. The artificial increase of the annulus perimeter, as mentioned in [19], is used to calibrate the T2Well DBHE model from the experimental pressure data from [13]. More details on the T2Well mathematical model used for the DBHE model are provided in [18].

Mathematical CFD Model
The governing equations are solved with the CFD code ANSYS Fluent 17.1.0, which is based on the finite volume method [30]. The k − ε realizable turbulence model is applied in the simulations. For an axisymmetric model, the continuity equation is written as The momentum equations are defined by the following equations: where u is the fluid velocity and μ its molecular viscosity. The energy equations in the fluid and solid zone are is the effective thermal conductivity + , with the turbulence thermal conductivity. is described by c /0.85 , with c and the specific heat capacity and the turbulent viscosity, respectively. The stress tensor is defined and can be retrieved from [30]. The energy E and the sensible enthalpy h are calculated via Equation (12): The standard turbulence k − ε model is based on transport equations for the turbulence kinetic energy (k) and its dissipation rate (ε) [30]. This model is fast and has proved to be efficient at high speed and suitable for numerous industrial applications. Continuum equations can be written as with the eddy viscosity , S, , described in [30]. and are respectively equal to 1.44 and 1.99 by default. The different terms present in the continuum equation for both T2Well and Fluent are reported in Table 3. and are the generation of turbulence kinetic energy, accounting for the mean velocity gradient, and due to buoyancy when the water properties are not constant, respectively. represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate [30], ∝ is the speed of sound in the fluid.
The turbulent heat transport uses the Reynold's momentum transfer such as Gravity is accounted for in the numerical simulation. A convergence criterion of 10 −8 was used for all equations of the flow solver, ensuring a good convergence of the flow solution. The simple scheme is selected for the pressure velocity coupling. A second order spatial discretization is used for the pressure. Turbulent kinetic, dissipation rate, energy, and momentum spatial discretization are set to the second order upwind. Table 2 summarizes the boundary conditions applied (pressure and temperature) in the CFD models. The roughness on the well walls is set to 2.5 × 10 −5 m.

Mesh Independence Study
A mesh independence study is performed to make sure the results do not depend on the mesh density. The pressure inlet boundaries are artificially set to 5 bar at 209.85 °C for the internal tubing and 110 °C for the annulus. The pressure outlet is fixed at 1 bar. In this section, all boundaries are set without heat flux except for the external casing (C2) where a constant temperature value of 250 °C is applied. Three meshes of 43,500, 171,656, and 514,968 cells are built and used in a steady-state solver. Figure 3a shows the radial temperature distribution in the middle of the model, which highlights a high degree of convergence between all three meshes. Figure 3b uses a narrower range of temperature which shows the incremental increase of temperature as the distance from the axis of the well increases. As expected, the trend is smoother when using more cells, but it shows that the intermediate and fine grid results are very close.
Based on the computational cost of a simulation, the intermediate mesh (171,656 cells) has been selected. This mesh ensures a smooth increase of the temperature, preventing divergence problems, which is a problem easily faced when using the IAPWS-IF97 formulation [31], which was applied through a user defined function (UDF) called at each time step of the simulation.  A similar methodology has been applied to the third CFD model at the bottom of the DBHE (CFD 1987). This model contains 809,676 cells.

Thermal Analysis
In the two CFD models corresponding to the depths of −500 m and −1900 m (CFD 500 and CFD 1900), the pressure inlet and outlet conditions are linked to an imposed mass flow rate. Thus, the pressure field is accurately constrained in the production and injection wells compared to the T2Well values. Figure 4 shows the vertical temperature in the center of production well (inner tubing) and in the injection well (annular, taken at a 6 cm radial distance from the axis) after 1 year. In the production well, the CFD-based temperatures at −495 m are comparable (less than 1 °C difference) to the values obtained with T2Well. However, in the annulus section, the calculated steady state temperature shows increased values in the range of 1.3 °C and 8.0 °C. In all cases except when having a mass flow rate of 2 kg/s at 500 m (Figure 4f), no temperature difference could be noticed between results based on the IAPWS-IF97 formulation and those obtained with constant water properties. In Figure 4f, however, the use of constant water properties seems to underestimate the temperature distribution by 1 °C. The very low increase of temperature between the bottom and top of the production well (upwarding fluid) in Figure 4c,e,g is in agreement with the T2Well trend.
While the inlet injections in all wells (bottom for the production well and top for the injection well) match the T2Well values, the outlet velocities are higher, with a maximum amplitude of 0.2 m/s for a mass flow rate of 2 kg/s, and 0.85 m/s for 10 kg/s, respectively. Table 4 shows the error between the results obtained when using constant density, thermal conductivity, and viscosity values and those obtained when using the IAPWS-IF97-based formulations. Results based on a constant density show an error less than 5% compared to those obtained with pressure temperature dependent densities. For the water thermal conductivity, the error increases with a low mass flow rate at a greater depth (up to 14% overestimated). Pressure and temperature values are higher, showing that the assumption of constant surface water properties should not be considered. The viscosity shows the highest differences, up to 70% overestimation in the production well.   Figure 5 shows the radial temperature distribution obtained with the CFD models at the depth of 500 m and 1900 m for a mass flow rate injection of 2 kg/s and 10 kg/s with constant (a) and the IAPWS-IF97 water properties (b), respectively. The radial temperatures in the well are homogeneous in the production well (0-0.0253 m), which can validate the use of the 1D area average calculation in T2Well, for this part. However, the injection well shows a small increase of temperature (up to 5 °C) next to Casing 2. As T2Well computes the injection well temperature very close to Casing 2 for a given elevation (z or x) at the radial position of 0.0797 m, an over-estimation of the temperature in the injection well by 2-4 °C is expected when applying a 1D approach.

Heat Transfer Analysis
The production well extracts the energy, with a surface heat flux on the casing wall higher than 100 W/m 2 . The radial Prandtl number in water at the depth of −500 m (a) and −1900 m (b) are presented in Figure 6. As the water specific heat is constant in this study, the Prandtl number reveals only the ratio of the viscosity and the thermal conductivity. This allows to identify the predominance of heat transfer convection (Pr >> 1) or conduction (Pr << 1). Figure 6 shows that Pr is higher than 1, with the lowest values obtained for a mass flow rate of 2 kg/s in Figure 6b. When constant thermal conductivity and viscosity values are applied, the Pr reaches 6.9. Using the IAPWS IF97 properties, a higher mass flow rate shows a high Pr number. The convective heat transfer is higher due to the high water velocity and larger temperature difference between the fluid and its surroundings. In the annulus at 500 m for both mass flow rates, the Pr number rises due to the high temperature drop between the annulus and the tubing shown in Figure  5.
Except for the 10 kg/s at 500 m case, the use of constant values in this study tends to overestimate the convective heat transfer. The Prandtl number shows a low decrease at the wall boundary in the injection well, outlining a higher conductive heat transfer effect in the vicinity of the well wall.

Pressure and Temperature Comparison
At the bottom of the DBHE, water returns to the surface from the annular well to the producing well. In this study, 3 m separate the bottom of the internal tubing from the bottom of the DBHE. Figure 7 shows the static pressure distribution in the bottom of the well for a mass flow rate injection of 2 kg/s (a) and 10 kg/s (b) with constant water properties, and with the IAPWS-IF97 properties (c)-(d).  Table 5 summarizes the results of pressure in the annular and temperature in the internal tubing at −1987 m obtained with the CFD model and with T2well. For the CFD-based results, using the IAPWS-IF97 properties implies a small decrease in the pressure distribution (<1000 bars) and similar temperature results with a mass flow rate of 2 kg/s. Increasing the mass flow rate with the IAPWS-IF97 properties shows a decrease of the water temperature in the internal tubing at −1987 m when compared to constant properties based results.
The pressure calculated with ANSYS Fluent is similar to the T2Well results in the case of 2 kg/s and higher by three bars at 10 kg/s. At low mass flow rate, the CFD temperatures are higher by 5 °C and in similar range at 10 kg/s. The pressure correction applied in T2Well, calibrated with the experimental data in [18], deviates at high mass flow rate, with an overestimation of 80 bars. Thus, the estimation of the circulating water power consumption in the DBHE from the CFD study, suggests a lower value than reported by the authors in [19].
With a mass flow rate of 2 kg/s, the static pressure in the fluid is lower than for 10 kg/s. The pressure decreases in the internal tubing compared to the annular well, with a local depression near the bottom of the internal tubing. This pressure difference is higher when considering a mass flow rate of 10 kg/s. Thus, the fluid is accelerated in this section of the DBHE (see Section 3.3.2).
Furthermore, the constant boundary conditions applied on the side of the CFD model might be slightly overestimated compared to the T2Well values. As the distance of the node N5 (see Figure 2) to the casing wall is not negligible, the temperature obtained with T2Well is potentially slightly higher than considered in the casing wall with the CFD model.  A lower fluid velocity is noticed for a low mass flow rate injection. When water returns to the surface, a very low velocity zone (5 cm) is observed on the external side of the internal tubing. The use of the IAPWS-IF97 properties for water leads to higher velocity values at 2 kg/s and slightly lower velocities at 10 kg/s. Figure 9 shows the turbulence kinetic energy at the bottom of the DBHE. The turbulence kinetic energy appears mainly localized in the low velocity zone (visible in Figure 8) and at the center of the annular well, where water accelerates. As expected, the turbulence kinetic energy is also higher when the injection mass flow rate is higher. For 2 kg/s and 10 kg/s, approximate pressure differences of 1 and 4 bar are obtained, respectively, between the pressure in the annular and the internal tubing. Due to higher velocities values at high mass flow rates, frictions are generated in the bottom of the well, increasing the pressure drop.

Discussion
Despite not having considered the specific heat of water as a pressure-temperature dependent property in ANSYS Fluent, the heat transfer calculation between the CFD model and T2Well is similar, with a more detailed description of the physics when considering the CFD approach. After a simulated time period of 1 year using the IAPWS-IF97, the maximum density difference against the average value in the radial axis is less than 2% and less than 9% for the viscosity, mainly in the injection well. For the radial axis, the velocity variation reaches up to 5% in the production well and 20% in the injection well. The stagnant boundary of the water at the wall (v = 0) shows a local conductive heat transfer, impacting the convective heat transfer when the roughness of the well is high with high velocity (lower mean temperature in the well than at the contact of the well wall in Figure 4). Thus, the 1D calculation shows a slight increase of the mean temperature in a cross section of the DBHE. However, the constant geothermal gradient value used in the CFD models is subject to discussion due to the distance to the wall, mentioned in Section 3.3.1.
In the formal analysis of the energy equation, the thermal conductivity in T2well is an area averaged value considering the fluid and the casing material. As the casing thermal conductivity is higher than water, the area averaged value calculated can be higher than the water value only. This might tend to slightly overestimate the heat transfer using T2Well.
Further development can be achieved in the CFD model to deliver a pressure dependent formulation for the specific heat of pure water.
The return pressure estimated with T2Well appears highly overestimated at high mass flow rates. Despite being calibrated with an experimental study, the artificial pressure constrains used in T2Well need to be equilibrated for high mass flow rates, to avoid the overestimation of the energy and power consumed to transport water in DBHEs in high geothermal gradients. By evaluating lower pressure losses in the return flow, the CFD modelling tends to lower the pressure drop between the annular and tubing flow compared to the 1D approach. Thus, the energy consumption to impose a high pressure in the annulus of a DBHE might be subject to a re-evaluation using detail numerical models, or experimental data, if available.
The existence of independent convective cells at the bottom of the internal tubing opens new perspectives to optimize their intensity, notably with high conductive materials on the wellbore [27], to stimulate self-convection following the energy demand. Using both 1D and 2D-axisymmetric approaches enable to strengthen the heat transfer and fluid flow estimation from unconventional geothermal well designs. Due to the computational costs associated with CFD, it is advised to apply such techniques after a large set of 1D simulations is complete, or when a more detailed analysis is required.

Conclusions
This double numerical approach has aimed to refine the fluid flow and heat transfer numerical calculation in a DBHE in high geothermal gradient. A conjugated numerical use of T2Well-1D well formulation with three 2D axisymmetric CFD models from a theoretical DBHE set at 1900 m in volcanic formations with 350 °C at the depth of 2000 m, has been performed. A CFD approach refines the 1D description, considering the vertical and radial velocity distribution along with turbulence and local wall boundary effects. The difference of results highlighted in this work allows to discuss both approaches. Local wall boundary effects due to the well roughness decrease the mean temperature in the production well. The use of ANSYS Fluent and T2Well opens new perspectives to investigate and engineer DBHEs with pure liquid water. To estimate the energy produced from DBHEs, the use of constant water properties in the CFD model underestimates the overall heat transfer at high temperatures (more than 80 °C) at low mass flow rate (2 kg/s). It is advised to consider pressure-temperature dependent formulations for the water properties, to reduce uncertainties in the numerical modelling of a DBHE in high geothermal gradient areas.
Self-convective water cells are observed at the bottom of the DBHE, with no significant changes with water flow rate. The pressure losses associated with the return flow at the bottom of the DBHE are refined with CFD compared to a 1D approach, corresponding to a large pressure difference of 80 bars at a flow rate of 10 kg/s. Previous evaluations of the energy extracted and consumed in a DBHE can be discussed with lower expected return pressures. The conjugated numerical approach significantly improves the understanding of fluid flow and heat transfer in DBHEs in high geothermal gradients environments, which will allow further technological advances.