Heat storage efficiency, ground surface uplift and thermo-hydro-mechanical phenomena for high-temperature aquifer thermal energy storage

Abstract

High-temperature aquifer thermal energy storage (HT-ATES) systems can help in balancing energy demand and supply for better use of infrastructures and resources. The aim of these systems is to store high amounts of heat to be reused later. HT-ATES requires addressing problems such as variations of the properties of the aquifer, thermal losses and the uplift of the surface. Coupled thermo-hydro-mechanical (THM) modelling is a good tool to analyse the viability and cost effectiveness of HT-ATES systems and to understand the interaction of processes, such as heat flux, groundwater flow and ground deformation. The main problem of this modelling is its high computational cost. We propose a dimensional and numerical analysis of the thermo-hydro-mechanical behaviour of a pilot HT-ATES. The results of this study have provided information about the dominant thermo-hydraulic fluxes, evolution of the energy efficiency of the system and the role of the hydraulic and thermal loads generated by the injection and extraction of hot water.

Introduction

Background and objectives

The injection and extractionn of thermal fluids into the underground can produce unwanted consequences such as variation of water level, alteration of temperatures, subsidence, uplift and seismicity. Examples of geoengineering facilities affected by these problems are carbon dioxide storage, radioactive waste storage, injection for enhanced oil recovery, geothermal energy and underground thermal energy storage. All these facilities may trigger coupled thermo-hydro-mechanical (THM) phenomena. An example is the German Power Plant of Landau. A leak of hot water from a well induced positive vertical displacements on the land’s surface and the shutdown of the plant in 2014. This uplift was detected by Heimlich et al. (2015).

Underground thermal energy storage (UTES) technologies store thermal energy, heat or cold, by injecting thermal energy into the underground during a period of high energy supply. The thermal energy is extracted during a period of high energy demand. Heat storage can contribute to the extension of low-carbon heat sources, reduce greenhouse gas emissions and afford flexibility in the management of supply and demand of heat (Cremer 2021).

One of the most common UTES is aquifer thermal energy storage (ATES). ATES is a bidirectional system that consists of one or more wells that inject or extract thermal energy into aquifers (Schüppler et al. 2019). ATES stands out for its high storage capacity and for being useful in both small and large facilities (Pellegrini et al. 2019; Fleuchaus et al. 2018). ATES is classified as high-temperature ATES (HT-ATES) when the temperature of the injected water is above 60 °C (Drijver et al 2012). In this temperature range, heat can be used directly (Drijver et al. 2012). HT-ATES is a clear example of a facility that can trigger THM phenomena. Because of this, a good understanding of the underground processes is necessary for obtaining safe and efficient systems.

THM modelling can help in studying the viability and cost effectiveness of geoengineering facilities dealing with geofluids. This type of modelling implies a nonlinear coupled initial boundary value problem, which needs to be solved numerically. The computational costs of THM calculations are high because of basically 2 reasons: a high number of degrees of freedom and a strong coupling among nonlinear processes (Kolditz et al. 2010). Dimensional analysis permits to verify the solution of the nonlinear equations by substituting the variables by new dimensionless variables (Hauke 2008). It is a simple way to plan scale experiments and to know what can happen in the study area.

Despite the large number of examples in the literature of THM modelling applied to geoengineering facilities (e.g., Rutqvist et al. 2005; Vilarrasa et al. 2014; Toprak et al. 2016), there are only few applications to UTES, such as that of Park et al. 2016 who studied the THM behaviour of a high-temperature thermal energy storage cavern. In HT-ATES, there are many studies on the TH behaviour of the subsurface (e.g., Jongchan et al. 2010; Mindel et al 2021) and only one simulates the HM behaviour of the underground focusing on the deformation at the ground surface (Birdsell and Saar 2021). There is a THM simulation of HT-ATES, but it does not consider uplift (Jin et al. 2022). The study of the THM behaviour of HT-ATES is important because mechanics can alter fluid flow and heat transport, modifying the performance and the efficiency of the system, and can lead to important uplift of the ground surface, affecting buildings and infrastructure.

Pilot or demonstration projects study the viability of HT-ATES projects. Pilot projects simulate real projects at a small size or with some limitations in order to see if they are technically and economically possible at real scale (Fleuchaus et al. 2018). HEATSTORE is a GEOTHERMICA ERA-NET European project which aims to reduce costs, to increase safety and performance of high-temperature underground thermal energy storage (HT-UTES) technologies and to improve the supply–demand management in heat networks (Cremer 2021). There are several pilot projects of HEATSTORE HT-UTES in Europe, of which four are HT-ATES. These HT-ATES pilot projects are located in Switzerland, Netherlands and Iceland.

The aim of this paper is to analyse and understand the mechanisms of uplift caused by heat storage by studying the THM response to a HT-ATES system. In addition, efficiency has to be quantified as it is a major issue in the heat storage and quantifies heat losses to the aquifer. For this purpose, we first carry out a dimensional analysis to understand the thermo-hydro-mechanical phenomena during hot water injection and cold water extraction coupled to deformations. Then, we present a 3D THM finite element model to simulate and analyse the subsurface response in terms of temperature, pressure and displacements distribution and evolution. We make a prediction of the long-term behaviour of the system to calculate the efficiency and the displacements over a decade of operation. Finally, we present a sensitivity analysis.

Description of the Bern project

This study is motivated by the Bern project, one of the pilot projects of HEATSTORE. The power plant of Forsthaus (Bern, Switzerland) has an excess of waste heat that is dissipated into the atmosphere. A HT-ATES system is planned close to the plant in order to store the excess of heat during summer and to reuse it during winter through a district heating network.

The HT-ATES system consist of one central well surrounded by four auxilliary wells each at a distance of 40 m (Fig. 1). This distance should approximate the thermal radius, that is, the position of the temperature front after a period of heat injection (240 days in our case) considering only advection (Bloemendal et al. 2018). At the central well heat is injected in summer and heat is produced in winter. The auxiliary wells help to limit flow and have the desired pressure in the aquifer. This configuration of wells is an alternative to the well-known two well configurations. This 5-spot system permits to concentrate and keep the stored heat in the region between the auxiliary wells. It is also used in other contexts, such as, enhanced oil recovery (Hy et al. 2019).

Fig. 1

HT-ATES project of Bern in summer (A) and winter (B)

During summer, the central well injects water at 90 °C while the auxiliary wells pump water at the same rate of 25 l/s. In winter, the auxiliary wells inject water of 50 °C, while the central well extracts water. The loading lasts longer than the unloading (8 months versus 4 months). More heat is injected than extracted despite the same injection rates are considered during all the year (a constant flow rate of 25 l/s).

The aquifer is located in the Lower Freshwater Molasse (USM), which is part of the Swiss Molasse Basin. This Basin is a thick Tertiary sedimentary formation originated by detrital filling of a subsidence basin that was caused by the uplift of the Alps (Platt and Keller 1992). In the study area, the USM consists of horizontal layers of sandstones embedded within marls and claystones (Keller 1992). Quaternary unconsolidated deposits of 150 m thickness overlay the USM (Driesner 2021).

The heat will be stored in a series of confined aquifers between 200 to 500 m below the surface. The mean geothermal gradient is 3 K/100 m, which make the temperature of the aquifer to approximately 20 °C (Driesner 2021).

Mathematical formulation

Governing equations for THM modelling

The injection and extraction of hot water into an aquifer represents a case with strongly coupled THM phenomena. The injection of hot water can cause a local increase in pressure, which coupled to the increased temperature, can induce the expansion of the medium (Rutqvist et al. 2005). Heating produces thermal expansion of the rock and native porous water. The increment of pressure and temperature tends to generate uplift, vertical displacement towards the surface, which could be partly counteracted by possible loss of cementation. The effect of pore pressure and temperature may increase the zone of deformation with respect to the purely mechanical case (Mokni et al. 2013).

Heat transport in HT-ATES is a combination of advection, convection, diffusion and dispersion and is influenced by thermal, hydraulic and mechanical parameters. Heat capacity of water and solid, and flow of liquid water control advection and dispersion. Heat capacity is a function of temperature and the flow of liquid water is controlled by temperature (through water viscosity), degree of saturation (hydraulic effect) and porosity (mechanical effect). Solid density affects heat capacity. Flow of liquid water and heat capacity depend on liquid density. Both densities are functions of temperature and liquid pressure. Diffusion or heat conduction depends on thermal conductivity which is a function of porosity and temperature. Convection is a function of liquid density, thermal conductivity and water viscosity.

THM formulation for ATES system should integrate heat transport (heat conduction, heat advection, heat convection and heat dispersion), water flow (liquid advection) and mechanical behaviour (behaviour of porous materials dependent on stresses and temperature). The THM formulation presented here is based on Olivella et al. 1994.

The liquid mass balance equation in the medium is expressed as:

$$\frac{D\left({\rho }_{l}\phi \right)}{Dt}+{\rho }_{l}\phi \nabla \cdot \left(\dot{\mathbf{u}}\right)+\nabla \cdot \left({\rho }_{l}{\mathbf{q}}_{l}\right)=0,$$

(1)

where $\frac{D\left(a\right)}{Dt}$ is the material derivative with respect to solid phase velocity $\dot{\mathbf{u}}$, $\phi$ is the porosity [L³L⁻³], $\dot{\mathbf{u}}$ is the solid velocity vector [LT⁻¹] and ${\mathbf{q}}_{{\varvec{l}}}$ is the liquid advective flux [LT⁻¹].

The material derivative is defined as:

$$\frac{D\left(a\right)}{Dt}=\frac{\partial \left(a\right)}{\partial t}+\nabla \left(a\right)\cdot \dot{\mathbf{u}.}$$

(2)

The solid mass balance equation can be written as:

$$\frac{D\phi }{Dt}=\frac{\left(1-\phi \right)}{{\rho }_{s}}\frac{D{\rho }_{s}}{Dt}+\left(1-\phi \right)\nabla \cdot \dot{\mathbf{u},}$$

(3)

where $\phi$ is the porosity, $\dot{\mathbf{u}}$ is the solid velocity vector, ${\rho }_{s}$ is the solid density [ML⁻³] and $D/Dt$ represents the material time derivative of the solid. Equation (3) represents the variation of porosity caused by solid density variation and volumetric deformation.

The internal energy balance of the medium can be expressed as:

$$\frac{\partial \left(\left(\phi {c}_{l}{\rho }_{l}+\left(1-\phi \right){c}_{s}{\rho }_{s}\right)\Delta T\right)}{\partial t}+\nabla \cdot \left(-\lambda \nabla T\right)+\nabla \bullet \left({{c}_{l}\rho }_{l}\Delta T{\mathbf{q}}_{l}\right)+\nabla \cdot \left(\left(\phi {c}_{l}{\rho }_{l}+\left(1-\phi \right){c}_{s}{\rho }_{s}\right)\Delta T\dot{\mathbf{u}}\right)=0,$$

(4)

where ${c}_{l}$ is the specific heat of liquid [FLM⁻¹θ⁻¹] and $\lambda$ is the thermal conductivity [FT⁻¹θ¹].

Equation (4) is simplified by using Eqs. (1) and (4), applying the material derivative (2) and doing some algebra:

$$R{\rho }_{l}\phi {c}_{l}\frac{DT}{Dt}+\nabla \cdot \left(-\lambda \nabla T\right)+{c}_{l}{\rho }_{l}{\mathbf{q}}_{{\varvec{l}}}\cdot \nabla T=0,$$

(5)

where $T$ is the temperature [θ], ${c}_{l}$ is the specific heat of liquid [FLM⁻¹θ⁻¹], ${c}_{s}$ is specific heat of solid [FLM⁻¹θ⁻¹], $\lambda$ is the thermal conductivity [FT⁻¹θ¹] and $R$ is the thermal retardation coefficient [−].

The characteristic thermal conductivity $\lambda$ is expressed as the sum of the bulk thermal conductivity ${\lambda }_{b}$ and the dispersive conductivity ${\lambda }_{d}$:

$$\lambda ={\lambda }_{b}+ {\lambda }_{d}.$$

(6)

The momentum balance of the porous media has to be satisfied in order to solve the mechanical part of the problem. Neglecting inertial terms, it is reduced to the equilibrium of stresses:

$$\nabla \cdot {\varvec{\upsigma}}+\mathbf{b}=0,$$

(7)

where ${\varvec{\upsigma}}$ [ML⁻¹ T⁻²] is the stress tensor and $\mathbf{b}$ [ML⁻² T⁻²] is the body forces vector.

We assume that the medium is elastic. We use linear thermoelasticity in porous media to take into account the effect of changes in fluid pressure and temperature on rock deformation. The elastic strain is based on Hooke’s law. We use the same formulation as Vilarrasa et al. (2013) for the strain, but we have added the thermal strain part:

$${\varvec{\upvarepsilon}}=\frac{{\sigma }_{m}^{^{\prime}}}{3K}\mathbf{I}+\frac{1}{2G}\left({{\varvec{\upsigma}}}^{\mathrm{^{\prime}}}-{\sigma }_{m}^{^{\prime}}\mathbf{I}\right)+ \frac{{\alpha }_{T}}{3} \mathrm{T}\mathbf{I},$$

(8)

where ${\varvec{\upvarepsilon}}$ is the elastic strain tensor [−], ${{\varvec{\upsigma}}}^{\boldsymbol{^{\prime}}}$ is the effective stress tensor [ML⁻¹ T⁻²], ${\sigma }_{m}^{^{\prime}}$ is the mean effective stress [ML⁻¹ T⁻²], $\mathbf{I}$ is the identity matrix, ${\alpha }_{T}$ is the volumetric (3 times the linear) thermal expansion coefficient of the medium [θ⁻¹], $K=E/3(2-\nu )$ is the bulk modulus [ML⁻¹ T⁻²], $G=E/(2(1+2\nu )$ is the shear modulus [ML⁻¹ T⁻²] and $E$ is the Young’s modulus and ν the Poisson ratio [−].

Using (8) and taking into account (i) the compatibility equations between strains and displacements, that (ii) the volumetric strain can be defined as the divergence of the displacement vector, (iii) assuming that loads are stationary and (iv) deriving with respect to time, the momentum balance of the porous media (7) can be written as:

$$G{\nabla }^{2}\dot{\mathbf{u}}+\left(K-\frac{2G}{3}\right)\nabla \left(\nabla \cdot \dot{\mathbf{u}}\right)-\left(K+\frac{2G}{3}\right){\alpha }_{T}\frac{\partial }{\partial t}\nabla T-\frac{\partial }{\partial t}\nabla {p}_{l}=0,$$

(9)

where ${p}_{l}$ is the liquid pressure [ML⁻¹ T⁻²].

A summary of the basic relationships for mass and momentum balance equations is given in Table 1.

Table 1 Basic relationships for mass and momentum balance equations

Dimensional analysis

In this paper, we perform a dimensional analysis of the equations of internal energy balance (5) and momentum balance (9). The main purpose is to obtain dimensionless variables or dimensionless numbers that better characterize the thermo-hydro-mechanical processes.

Dimensional analysis is based on Buckingham π theorem (Buckingham 1914). This theorem establishes that for any dimensionally homogenous equation that involves n dimensional variables is possible to find an equivalent dimensionless equation with n-d dimensionless variables, where d is equal to or less than the number of independent dimensions of the problem (Hauke 2008).

The way to apply the dimensional analysis is by changing variables (Hauke 2008). The original variables of Eqs. (5) and (9) are written as a function of a characteristic variable, which is a dimensional constant, and a dimensionless variable. We used the following dimensionless variables:

$${t}_{D}=\frac{t}{{t}_{C}} ; {T}_{D}=\frac{T-{T}_{0}}{{T}_{C}} ; {p}_{lD}=\frac{{p}_{l}-{p}_{0}}{{p}_{lC}}; {\mathbf{u}}_{D}=\frac{\mathbf{u}}{{\mathrm{u}}_{C}}.$$

(10)

In each variable, subscript D refers to the dimensionless variable and subscript with C refers to the characteristic variable.

Equations (5) and (9) can be written in terms of the dimensionless variables of (10) as follows:

$$\frac{D{T}_{D}}{D{t}_{D}}-{\frac{{Pe}^{-1}}{2}\nabla }_{D}\cdot \left({\nabla }_{D}{T}_{D}\right)+\frac{\left|{\mathbf{q}}_{{\varvec{l}}}\right|\pi {L}_{c}{b}_{aq}}{\left|\mathrm{Q}\right|}\cdot {\nabla }_{D}{T}_{D}=0,$$

(11)

$${\nabla }_{D}^{2}{\dot{\mathbf{u}}}_{D}+\frac{2}{3}\left(\frac{3\nu -1}{2-\nu }\right){\nabla }_{D}\left({\nabla }_{D}\cdot {\dot{\mathbf{u}}}_{D}\right)-\frac{2}{3}\left(\frac{\nu +3}{2-\nu }\right){U}_{T}\frac{\partial }{\partial {t}_{D}}{\nabla }_{D}{T}_{D}-\frac{2}{3}\left(\frac{1+2\nu }{2-\nu }\right){U}_{P}\frac{\partial }{\partial {t}_{D}}{\nabla }_{D}{p}_{lD}=0.$$

(12)

Assuming radial flow (${q}_{l}=\frac{\left|\mathrm{Q}\right|}{2\pi r{b}_{aq}}$), Eq. (11) can be written as:

$$2\frac{D{T}_{D}}{D{t}_{D}}-{{Pe}^{-1}\nabla }_{D}\cdot \left({\nabla }_{D}{T}_{D}\right)+\frac{{L}_{C}}{r}\cdot {\nabla }_{D}{T}_{D}=0.$$

(13)

The characteristic length ${L}_{c}$ is equal to the distance between the central well and the auxiliary well [L] and $r$ is the radial distance from the central well [L]. The characteristic time ${t}_{C}$ represents the characteristic time of the advection equation:

$${t}_{C}=\frac{{L}_{c}^{2}\phi R}{2r\left|{\mathbf{q}}_{{\varvec{l}}}\right|}.$$

(14)

$Pe$ is the dimensionless Peclet number [−]. We use this number to determine which thermo-hydraulic fluxes is the dominant one (advection, conduction or dispersion) and where the Peclet number represents the ratio between heat transfer by advection and heat transfer by diffusion and dispersion. We define Pe as:

$$Pe={\frac{{\left|Q\right|\rho }_{l}{c}_{l}}{2\pi {b}_{aq}({\lambda }_{b}+{\lambda }_{d})}=\left(P{e}_{cond}^{-1}+P{e}_{disp}^{-1}\right),}^{-1}$$

(15)

where $Q$ is the injection flow rate [L³T⁻¹] $\mathrm{and }{b}_{aq}$ is the thickness of the injected aquifer [L].

In Eq. (15), $P{e}_{cond}$ is the the ratio between heat transfer by advection and by diffusion [−]:

$$P{e}_{cond}=\frac{{\left|Q\right|\rho }_{l}{c}_{l}}{2\pi {b}_{aq}{\lambda }_{b}}.$$

(16)

We define $P{e}_{cond}$ in (16) for one well. In a 5-spot scheme, we express $P{e}_{cond}$ as:

$$P{e}_{cond}=\frac{{{q}_{l}r\rho }_{l}{c}_{l}}{{\lambda }_{b}}.$$

(17)

$P{e}_{disp}$ is the ratio between heat transfer by advection and by dispersion [−]:

$$P{e}_{disp}=\frac{r}{{d}_{l}}.$$

(18)

${U}_{T}$ is a dimensionless number which expresses the ratio between the thermal strain generated by the injection of hot water and the initial total strain [−]:

$${U}_{T}=\frac{{\frac{{\alpha }_{T}}{3}{T}_{C}L}_{c}}{{u}_{C}}.$$

(19)

${U}_{P}$ is a dimensionless number which is equal to the ratio between the hydraulic strain generated by the injection of hot water and the initial total strain:

$${U}_{P} =\frac{\frac{{p}_{c}{L}_{c}}{K}}{{u}_{C}}.$$

(20)

In (19) and (20), ${u}_{C}$ is the characteristic displacement:

$${u}_{C}= {p}_{0}S{L}_{c}+ {T}_{0}\left({\alpha }_{lT}\phi +{\alpha }_{sT}\right){L}_{c},$$

(21)

where $S=\left({\alpha }_{lp}\phi +{\alpha }_{sp}\right)$ is the storage coefficient [−] and ${T}_{0}$ is the temperature of the medium before injection [θ].

The characteristic temperature ${T}_{C}$ [θ] is defined as:

$${T}_{C}={T}_{inj}-{T}_{0,}$$

(22)

where ${T}_{inj}$ is equal to the temperature of the injected water [θ].

The characteristic liquid pressure ${p}_{C}$ [ML⁻¹ T⁻²] is equal to:

$${p}_{C}={p}_{inj},$$

(23)

where ${p}_{0}$ is the liquid pressure in the medium previous to the injection [ML⁻¹ T⁻²] and ${p}_{inj}$ is the injected liquid pressure [ML⁻¹ T⁻²] defined as:

$${p}_{inj}={\rho }_{l}g\Delta h,$$

(24)

where $g$ is the gravitational acceleration [LT⁻²] and $\Delta h$ is the variation in the hydraulic head defined by Thiem (1906):

$$\Delta h=\frac{Q}{2\pi {b}_{aq}k}\mathrm{ln}\left(\frac{{r}_{well}}{{L}_{c}}\right),$$

(25)

where ${r}_{well}$ is the well radius [L] and $k$ is the intrinsic permeability [L²].

We evaluate the importance of convection in the HT-ATES system with the dimensionless Rayleigh number (${R}_{a}$). ${R}_{a}$ is the ratio between the characteristic time for conductive transport and that for buoyant thermal transport:

$${R}_{a}=\frac{{t}_{cond}}{{t}_{grav}} =\frac{{{b}_{aq}}^{2}/\left(\frac{{\lambda }_{b}}{{Rc}_{l}{\rho }_{l}\phi }\right)}{{b}_{aq}/\left(\frac{\Delta {\rho }_{l}kg}{{\mu }_{l}}\right)}=\frac{{b}_{aq}{Rc}_{l}{\rho }_{l}\phi \Delta {\rho }_{l}kg}{{\mu }_{l}{\lambda }_{b}},$$

(26)

where $\Delta {\rho }_{l}$ is the fluid density difference [ML⁻³] and ${\mu }_{l}$ is the dynamic viscosity of liquid [ML⁻¹ T⁻²].

There is a critical Rayleigh number ${R}_{ac}$ equal to 4π² for a case of convective flow in a hexahedral domain of porous material saturated with fluid, (Lapwood 1948; Beck 1972). If ${R}_{a}$ is higher than ${R}_{ac}$, heat convection is the dominant process, else heat conduction is the dominant one. ${R}_{ac}$ only considers free convection as heat transport process.

Energy recovery

The energy recovery factor permits to quantify the efficiency of the HT-ATES system from the injected and extracted heat. This factor is expressed as (Doughty et al. 1982):

$$\eta =\frac{{\int }_{0}^{t}{\rho }_{l}{c}_{l}\left|{Q}_{ext}\right|\left({T}_{ext}-{T}_{0}\right)dt}{{\int }_{0}^{t}{\rho }_{l}{c}_{l}\left|{Q}_{inj}\right|\left({T}_{inj}-{T}_{0}\right)dt},$$

(27)

where ${T}_{prod}$ is the temperature of the extracted water [θ],${T}_{inj}$ is the temperature of the injected water and [θ] ${T}_{0}$ is the initial temperature in the aquifer [θ].

According to Gutierrez-Neri et al. (2011), the energy recovery factor is related to the Rayleigh number. A HT-ATES system is more efficient when heat conduction dominates over heat convection. For a ${R}_{a}$ lower than the ${R}_{ac}$, the energy recovery factor is expected to be higher than 0.6.

Numerical model

The numerical analysis has been performed with the Finite Element Method program CODE_BRIGHT, which is a capable of solving coupled thermo-hydro-mechanical models in geological media (Olivella et al. 1996). CODE_BRIGHT can be used in combination with the pre/post-processor GiD, developed by the International Center for Numerical Methods in Engineering (CIMNE 2020). CODE_BRIGHT uses the formulation presented above: conservation Eqs. (1), (3), (4) and (9). The phenomena considered (deformations, advection, conduction and dispersion) can be summarized in view of the constitutive equations defined above, namely: linear thermo-poro-elasticity, generalized Darcy’s law (including variable density and viscosity) and Fourier`s law (including dispersion). The equations are solved by means of finite elements in space and finite differences in time. The case of study is nonlinear and for this reason, the program uses the method of Newton–Raphson. The solution is obtained monolithically. In addition to other verifications and validations which can be found in the literature (Rutqvist et al. 2005; Olivella and Gens 2005; Tamayo et. al. 2021), CODE_BRIGHT has been verified against heat storage TH problems by Mindel et al. (2021).

Geometry and mesh

We consider a HT-ATES system that is based on the Bern project described before. It assumes a 3D multi-layer system. The model represents a quarter part of the whole domain, which consists of 1 central borehole and 4 auxiliary boreholes situated at 40 m from the centre (Fig. 2A, B). So, the model is simplified to 2 vertical wells: a central well and an auxiliary well. The geological profile consists of a top (0––400 m depth), an aquifer (400–500 m depth) and a bottom rock (500–600 m depth). Top rock and bottom rock represent aquitards (Fig. 2A).

Fig. 2.

3D model (A), plant view of wells (B) and geometry of the wells used in the simulation (C). C represents the central well and A the auxiliary wells

We have modelled the well screen explicitly. The wells have an octagonal shape as an equivalent of a round well (Fig. 2C). The inner part of the wells is empty. The parts of the wells are screened at the aquifer. This is simulated by a high permeability and dispersion in the outer part of the well with respect to those of the aquifer. The parts of the wells at the top and bottom rock are closed, which is simulated by outer parts having the same material properties as the surrounding rock.

The model boundaries are far enough so that effects of injection and extraction at these boundaries can be neglected. The size in the horizontal direction is equal to 200 m and the height is 600 m. The length of the wells is 500 m. We have used a semi-structured mesh: unstructured in the horizontal and structured in the vertical direction. The elements of the mesh are hexahedral (Fig. 2A). A total of 33027 nodes and 27300 elements have been used.

Material properties

Table 2 shows a summary of the main parameters of the geological materials. It corresponds to an initial model where all the materials are isotropic. Mechanical properties are the same for all geological materials, except for the well screen at the injection zone. We have considered a linear elastic model for all materials.

Table 2 Main parameters for preliminary calculations

Boundary and initial conditions

Water is injected and extracted into/from the nodes representing the wells within the aquifer. The total flow rate (25 l/s) is divided by 4 because the model only considers a quarter part of the domain. Injection (into the central well during injection and into the auxiliary well during back injection) is simulated by a fixed flux. Extraction (from the auxiliary well during injection and into the central well during back injection) is simulated by a fixed pressure of 6 MPa. Injection in the central well is carried out at 90 °C, while injection in the 4 auxiliary wells (25/4 l/s at each one) is carried out at 50 °C. We assume that injection in the central well takes place during two-thirds of the year. Flow is reversed during the remaining one-third of the year.

We impose pressure and temperature on the top surface at atmospheric pressure (0.1 MPa) and 20 °C, which can be considered atmospheric conditions.

We have applied zero vertical displacements on the bottom surface and zero horizontal displacements at the lateral boundaries of the model. For external boundaries (non-symmetry lateral planes), the vertical displacements are also prescribed to zero. We have assumed that lateral boundaries and bottom surface are impervious and adiabatic.

The initial temperature in the entire model is 20 oC and the liquid pressure is assumed hydrostatic starting with atmospheric pressure at the top.

Time interval data

We simulate 10 years of operation of the system in order to predict the behaviour of the system for a long time. Each year of operation consists of an injection (first 8 months) and a back injection (4 last months):

• During injection water is injected in the central well with a constant flow rate of 6.25 l/s and a temperature of 90 oC. Water is extracted from the auxiliary wells at the same time.

• During the back injection thermal energy is extracted. It consists of injecting water at 50 oC with a constant flow rate equal to 6.25 l/s in the auxiliary wells (only one auxiliary well in the domain). Water is extracted from the central at the same time.

There is an initial stage, previously to the operation time, that lasts ten days. This initial stage is used to calculate the initial stresses and the initial liquid pressures in the model.

Results for the base case

Thermo-hydraulic fluxes

Before discussing the results of the numerical model, we analyse a perfectly radial case with only injection (or extraction) from a central well without extraction (or injection) elsewhere. Flow lines are radial as in the single well scheme of Fig. 4. The injection flow rate is 25 l s⁻¹ and the thickness and properties of the aquifer are equal to the studied HT-ATES system. Figure 3 shows the Peclet number ($Pe$) is decomposed in the conduction term (${Pe}_{cond}$) and the dispersion term ($P{e}_{disp}$) as defined by Eqs. (15), (16) and (17). For this radial case $P{e}_{cond}$ is constant and $P{e}_{disp}$ increases linearly with radial distance. $P{e}_{disp}$ and $P{e}_{cond}$ intersect where radial distance is r = d_l $P{e}_{cond}$. Beyond this point conduction dominates over dispersion. For $Pe$ higher than 1, advection is the most dominant process. $Pe$ is almost equal to $P{e}_{disp}$ in the first metres from the well. From the centre to a distance equal to the dispersivity (d_l), dispersion is the most dominant process. Advection dominates when r > d_l. $Pe$ tends to $P{e}_{cond}$ at larger distances. In injections with very low injection rates, $P{e}_{cond}$ would be less than 1 and advection would never predominate in the system.

Fig. 3

Perfect radial $\mathrm{Pe}$ assuming a purely radial injection

In the numerical model of the HT-ATES case, water is simultaneously injected and extracted according the five-spot scheme. As a result, flow is not perfectly radial and flow lines are curved and finish behind the auxiliary wells (Fig. 4). Moreover, the flow drops to zero at larger distances from the well because the effect of injection (or extraction) at the central well is compensated by an extraction (or injection) at the auxiliary wells. In this case, we are interested in knowing the Peclet number along a central line, passing through the auxiliary well, and a border line, passing through the middle of two auxiliary wells (Fig. 4).

Fig. 4

Flow lines for single well and five-spot scheme

The Peclet number for HT-ATES is similar to Peclet number of the perfectly radial case within the first 40 m, along the central line (Fig. 5). The first thing to note is that the Peclet number drops beyond the auxiliary wells. This is due much smaller fluxes in this zone for the five-spot scheme compared to the perfectly radial single well scheme. Figure 4 displays the distribution of fluxes qualitatively. Another difference between the two cases is that the distance r = d_l $P{e}_{cond}$ is shorter than in the perfectly radial case. This distance is quite similar in the injection and in the back injection, a few metres from the auxiliary well. In the five well scheme (HT-ATES case) a distance can be defined where conduction dominates over the other thermo-hydraulic fluxes. The dominance of conduction is reached at shorter distance during injection than during back injection. The graphs of Pe for 10 years of operation are quite similar to those for 1 year, but the r = d_l $P{e}_{cond}$ point is at larger distance in the tenth back injection (Fig. 6).

Fig. 5

Peclet number at the end of the injection (red colour) and at the end of the back injection (blue colour) in the first year of operation. Aw represents the auxiliary well. Black and grey lines represent the perfectly radial Peclet. Central well is located at 0 m and the auxiliary well is represented with A_w

Fig. 6

Peclet number at the end of the injection (red colour) and at the end of the back injection (blue colour) in the tenth year of operation. Central well is located at 0 m and the auxiliary well is represented with A_w. Black and grey lines represent the perfectly radial Peclet

The curves of Pe along the border line are not very different from those along the radial line case (Fig. 7). Conduction dominates beyond the same distance during injection and back injection of the first year of operation.

Fig. 7

Peclet number at the end of the injection (red colour) and at the end of the back injection (blue colour) in the first year of operation in the border line. Black and grey lines represent the perfectly radial Peclet. Central well is located at 0 m

We have calculated the maximum Rayleigh number for the maximum temperature increment to be expected, which is equal to the maximum temperature injected (90 °C) minus the initial temperature in the aquifer (20 °C). The increment of the liquid density is equal to the difference between the liquid densities at these two temperatures. Using Eq. 25 and parameters and properties from Tables 1 and 2, we can calculate a Rayleigh number of 8.46, lower than the critical one (4π²). Therefore, we do not expect convection to affect the efficiency a lot.

Liquid pressure

Figure 8 shows the evolution of the liquid pressures in the aquifer near the central and auxiliary wells simulated in the THM model after 10 years of operation. During the injection, the liquid pressure is higher in the central well than in the auxiliary well, which means that water flows from the centre of the system to the outside. The behaviour is inverted during the back injection. The maximums of this graph are reached at the beginning of the injection and the minimums at the beginning of the back injection.

Fig. 8

Evolution of the liquid pressures on the aquifer near central well (C_aqu) near the auxiliary well (A_aqu)

Energy storage

The evolution of the temperature near the 2 wells has different behaviour according to the stage of the operation (Fig. 9). During the injection of the first year, the temperature increases up to 90 °C in the central well (temperature of the injected water), and to 59 °C in the auxiliary well. In the back injection stage, the temperature decreases to 65 °C in the central well and it remains constant in the auxiliary well (water is injected at 50 °C in this well). During the 10 years simulated, after each injection the temperature in the auxiliary well reaches a maximum, which increases with time. The temperature in the central well at the end of each year increases from 65 °C in the first year to 77 °C in the last year.

Fig. 9

Evolution of temperature during 10 years of operation on the aquifer near central well (C_aqu) and near the auxiliary well (A_aqu)

Figure 10 shows the evolution of the isotherm in the aquifer. Heat concentrates mainly between the two wells, although part of it escapes vertically and horizontally, producing heat losses. Heat escapes only by conduction through the top and the bottom of the aquifer, ascending and descending from 30 m at the end of 1 year (first back injection) to 70 m (tenth back injection). Horizontally, the thermal front moves according to the processes explained in the Peclet graphs (Figs. 5 and 6) up to 140 m at the end of the tenth back injection.

Fig. 10

Isotherms in the aquifer in a radial vertical section from central well (Cw) to one auxiliary well (Aw) at different times: end of the injection and end of the back injection for the first year and the tenth year. Vertical dotted lines represent the wells in the aquifer

In the studied system, the energy recovery factor grows about twice as fast in the back injection as in the injection (Fig. 11). The aquifer is storing more and more heat with time. The energy recovery factor increases to 0.89 asymptotically for a decade of operation. This value is in agreement with the relation between ${R}_{a}$-$\eta$ proposed by Gutierrez-Neri et al. (2011).

Fig. 11

Energy recovery factor during 10 years of operation

Strains and displacements

The strains are a result of both hydraulic and thermal loads, which in turn are a consequence of the injection of hot water. The U_T and U_P numbers of Eqs. (19) and (20) are a preliminary approach to see the effect of these loads in the aquifer.

Figure 12 points out that during injection U_T (90 °C) is higher than U_P, that is, the effect of temperature dominates over the liquid pressure. U_P is equal for the injection and back injection (same flow rates). In the back injection U_P dominates over U_T (50 °C). U_T and U_P decrease with depth due to the increment of the hydrostatic pressure.

Fig. 12

U_T and U_P in the injection and in the back injection in the aquifer. U_P is the same in the two different stages

Figure 13 illustrates the results of the strain in the aquifer of the THM model. The variations in thermal strain are proportional to the changes in temperature (Fig. 13A). The same occurs for the hydraulic strain with the variation in liquid pressure (Fig. 13B). Both types of strains are in the same order of magnitude. The sum of both strains confirms that the total strain in the central well is higher than in the auxiliary wells (Fig. 13C). Total strain increases during the injection in both wells, but during back injection it decreases in the central well and increases in the auxiliary well.

Fig. 13

Thermal strain (A), hydraulic strain (B) and total strain (C) in the central well (C_aqu) and in the auxiliary well (A_aqu) in the aquifer

The vertical displacements of the model are higher on the surface than in the aquifer (Fig. 14). The maximum displacements are reached at the end of the year. As the heat is largely maintained within the aquifer, the displacements are mainly due to hydraulic loads on the surface. In order to see this effect, we have made a THM model where the medium is not deformed by thermal loads (Fig. 14). In this model, the volumetric thermal expansion coefficients of water and solid have been reduced by a factor of 1000. The vertical displacements of this model on the top of the central well are almost the same as in the normal case, hence the contribution of the thermal loads is very low there. Instead, in the aquifer, the vertical displacements are more than half than those of the base case, which confirms that the deformations are due to thermal and hydraulic loads in the aquifer.

Fig. 14

Vertical displacements with thermal expansion and without thermal expansion on the top of the central well (C_srf) and on the top of the aquifer on the central well (C_aqu)

Surface displacements are strongly affected by the mechanical boundary conditions of the model further laterals. Fixed horizontal displacements would give higher and almost vertical displacements. We have considered a larger model with the vertical and horizontal displacements fixed in the further laterals. The results of the model show that the vertical displacements tend to concentrate in the central well and tend to be 0 with the distance (Fig. 15). The increment of the vertical displacement is 3 times higher during injection than during back injection. For a better assessment of the vertical displacement field experiments and model calibrations are necessary. We have obtained very low horizontal displacements (lower than 1 mm) at the surface because of the symmetry of the studied problem. The vertical displacements for 10 years are approximately the same as for 1 year, because we are considering a linear elastic model for all materials.

Fig. 15

Surface vertical displacements at the end of the injection and at the end of the back injection. C_srf represents the top of the central well and A_srf the top of one auxiliary well

Sensitivity analysis

The study of underground phenomena implies a large number of uncertainties on parameters in the initial stages of a project. These uncertainties can be divided in random and knowledge-based (Les Landes et al. 2021). We can reduce the second ones with analytical and numerical models. We study the sensitivity of the intrinsic permeability of the reservoir and the injection and extraction flow rate (quantity of thermal energy) to thermo-hydraulic fluxes in the aquifer, liquid pressures, energy recovery factor, strains and surface displacements. We analyse an aquifer with half the intrinsic permeability (0.5 k) of the base case (10^–13 m²), and another one with 10 times the intrinsic permeability (10 k). We also examine two cases where the injection flow rate and the pumping flow rate are double (2Q_i) and half (0.5Q_i) that of the base case (Q_i = 25 ls⁻¹).

Thermo-hydraulic fluxes

A change in intrinsic permeability or flow rate does not produce any change in the Pe number between the central and auxiliary well (Fig. 16). The main difference is beyond the auxiliary wells (Fig. 16A, B). A higher flow rate shifts the zone of dominant conduction further away. With respect to intrinsic permeability, a more permeable aquifer makes advection more dominant beyond the auxiliary well (Fig. 17).

Fig. 16

Simulated Pe in the injection for different intrinsic flow rates of injection (0.5q, q and 2q) during the injection (A) and back injection (B)

Fig. 17

Simulated Pe in the injection for different intrinsic permeabilities (0.5 k, k and 10 k) in the aquifer during the injection (A) and back injection (B)

In the same way as in the base case, we have calculated the Rayleigh number for the different intrinsic permeabilities. Ra for 0.5 k is 4.23 and for 10 k is 84.61. For the 10 k case convection may become relevant.

Liquid pressure

Liquid pressures in the aquifer are inversely proportional to the intrinsic permeability (Fig. 18), and directly proportional to the injection flow rate (Fig. 19).

Fig. 18

Liquid pressures on the top of the aquifer for different intrinsic permeabilities (0.5 k, k and 10 k) in the aquifer: central well (A) and auxiliary well (B)

Fig. 19

Liquid pressures on the top of the aquifer for different injection flow rates (0.5q, q and 2q): central well (A) and auxiliary well (B)

Energy recovery

The energy recovery factor does not present important differences when the intrinsic permeability of the aquifer is changed for 10 years of simulation (Fig. 20 A). A higher injection/pumping rate increases the injected heat as well as the energy recovery factor (Fig. 20 B).

Fig. 20

Energy recovery factor for different intrinsic permeabilities (0.5 k, k and 10 k) of the aquifer (A) and different injection flow rates (0.5q, q and 2q) (B)

Displacements

Since we do not modify the temperature of the injected water in any of the sensitivity case, U_P is the only strain component that changes. The U_P is inversely proportional to the intrinsic permeability (Fig. 21A). If we increase the flow rate in the aquifer, we observe more hydraulic strain and U_P (Fig. 21B).

Fig. 21

U_p for different intrinsic permeabilities (0.5 k, k and 10 k) of the aquifer (A) and for different injection water flow rate (0.5q, q and 2q) (B)

The surface vertical displacements are a consequence of events in the aquifer. Higher values of U_P in the aquifer, lead to higher surface displacements (Figs. 22 and 23). In all simulated cases, these displacements do not vary much during the injection, between 3 and 4 cm (Figs. 22A and 23A). The variation is higher during the back injection, because the hydraulic loads are more important than the thermal loads in the aquifer (Figs. 22B and 23B).

Fig. 22

Surface vertical displacements for different intrinsic permeabilities (0.5 k, k and 10 k) of the aquifer in the injection (A) and in the back injection (B)

Fig. 23

Surface vertical displacements for different injection flow rates (0.5q, q and 2q) in the injection (A) and in the back injection (B)

Discussion

The main goal of this paper is to investigate the uplift, energy efficiency and THM phenomena of a HT-ATES system composed of one central well and four auxiliary wells. One year of operation is divided into two stages: injection (injection of hot water) and back injection (extraction of hot water).

High flow rates dominate in the area of the aquifer between the central and auxiliaries wells. Heat is mainly transported by dispersion in the initial metres and by advection in the rest of this area. Conduction governs beyond the auxiliaries wells. Above and below the aquifer, conduction is the dominant heat flux.

The energy recovery factor of the system decreases during injection and increases during back injection except for the first year. More importantly, it increases on a yearly basis. Because of the high initial investment in research of the site, modelling and construction of the facility (Cremer 2021), it is essential to know when the system becomes more efficient.

Both thermal loads and hydraulic loads have an important effect on the displacements within the aquifer. At the surface, the vertical displacements are only a result of the hydraulic strains generated by the injection of water into the aquifer. Nevertheless, thermal loads could impact the surface in shallow aquifers.

The intrinsic permeability is a key parameter for HT-ATES and other geoengineering facilities. The permeability of the aquifer should be in an adequate range of magnitude for an efficient system without important heat losses caused by density-driven groundwater heat flow (Driesner 2021) and for avoiding unsafe uplift. Of course, knowing the intrinsic permeability of an aquifer is not always easy. Many geothermal reservoirs consist of separated layers or, in some cases, fractured rocks. They are usually modelled as an equivalent porous media aquifer (Birkholzer et al. 2008). In our study, we consider a constant intrinsic permeability for the whole aquifer. High permeabilities of the aquifer decreases hydraulic strains and uplift and enlarges the domain where advection dominates. Although Rayleigh number is proportional to the intrinsic permeability, it does not affect much the energy efficiency. This is more or less consistent with the results of Gutierrez-Neri et al. (2011) who do not expect important effects for Rayleigh number larger than a critical value of 4π². Only our model with the highest permeability slightly exceeds the critical Rayleigh number.

Injecting more heat in the aquifer will produce more efficient systems, but it has the disadvantages of enlarging the volume of underground affected by advection, causing more hydraulic strain and more uplift.

Future work should include more data and specific details of the site of study in order to obtain a better knowledge of the THM processes in the studied case.

Conclusions

We present a study that combines a dimensional analysis with a 3D THM finite element modelling in order to study the consequences of heat storage in HT-ATES. We focus on the mechanisms of uplift, heat transport processes and energy efficiency.

The dimensional analysis has been developed from a coupled THM formulation. We obtain three dimensionless numbers that study different processes: Peclet number (heat fluxes), ${U}_{P}$ (effect of hydraulic strain) and ${U}_{T}$ (effect of thermal strain). The Peclet number, and its decomposition in conduction and dispersive terms, permits to define not only the existence of the different heat fluxes (advection, conduction and dispersion), but also the area where each heat flux dominates. ${U}_{P}$ quantifies the effect of the hydraulic strain generated by the flow rate injected and ${U}_{T}$ quantifies that of the thermal effect by the introduction of hot water into the aquifer.

The study shows that advection and dispersion dominate inside the zone defined by the auxiliary wells of the aquifer. For larger radial distances conduction is the main heat transport process. Since convection is not very significant, the energy efficiency of the system is not affected negatively. The system becomes more efficient in terms of energy recovery after longer time in operation. Both thermal and hydraulic loads have an important effect on the displacements within the aquifer. Uplift is concentrated near the central well at the surface.

The presented dimensional and numerical analysis can be extended to other applications related to the injection and extraction of hot water into/out of the underground.