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

Heat storage efficiency, ground surface uplift andthermo‑hydro‑mechanical phenomena forhigh‑temperature aquifer thermal energy storage Rubén Vidal1* , Sebastià Olivella1,2, Maarten W. Saaltink1,3 and François Diaz‑Maurin1 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). 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, ground‑ water 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. Keywords: HT‑ATES, Heat, Modelling, Dimensional analysis, Thermo‑hydraulic fluxes, Uplift Open Access © The Author(s) 2022. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the mate‑ rial. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. RESEARCH Vidaletal. Geothermal Energy (2022) 10:23 https://doi.org/10.1186/s40517‑022‑00233‑3 Geothermal Energy *Correspondence: [email protected] 1 Department of Civil and Environmental Engineering (DECA), School of Civil Engineering, Universitat Politècnica de Catalunya— BarcelonaTech (UPC), UPC Campus Nord, Jordi Girona 1‑3, Building D2, 08034 Barcelona, Spain 2 International Center for Numerical Methods in Engineering (CIMNE), Barcelona, Spain 3 Associated Unit: Hydrogeology Group (UPC‑CSIC), Barcelona, Spain Page 2 of 27 Vidaletal. Geothermal Energy (2022) 10:23 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 HTATES, 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. Page 3 of 27 Vidaletal. Geothermal Energy (2022) 10:23 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). 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 Fig. 1 HT‑ATES project of Bern in summer (A) and winter (B) Page 4 of 27 Vidaletal. Geothermal Energy (2022) 10:23 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: where D(a) Dt is the material derivative with respect to solid phase velocity ˙u , φ is the porosity [L3L−3], ˙u is the solid velocity vector [LT−1] and ql is the liquid advective flux [LT−1]. (1) D(ρ l φ) Dt +ρlφ∇·(˙u)+∇·(ρlql)= 0, Page 5 of 27 Vidaletal. Geothermal Energy (2022) 10:23 The material derivative is defined as: The solid mass balance equation can be written as: where φ is the porosity, ˙u is the solid velocity vector, ρs is the solid density [ML−3] 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: where cl is the specific heat of liquid [FLM−1θ−1] and  is the thermal conductivity [FT−1θ1]. Equation(4) is simplified by using Eqs.(1) and (4), applying the material derivative (2) and doing some algebra: where T is the temperature [θ], cl is the specific heat of liquid [FLM−1θ−1], cs is specific heat of solid [FLM−1θ−1],  is the thermal conductivity [FT−1θ1] and R is the thermal retardation coefficient [−]. The characteristic thermal conductivity  is expressed as the sum of the bulk thermal conductivity b and the dispersive conductivity d : 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: where σ [ML−1 T−2] is the stress tensor and b [ML−2 T−2] 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: where ε is the elastic strain tensor [−], σ ′ is the effective stress tensor [ML−1 T−2], σ ′ m is the mean effective stress [ML−1 T−2], I is the identity matrix, αT is the volumetric (3 (2) D(a) Dt = ∂(a) ∂t +∇ (a) ·˙ u . (3) Dφ Dt = (1−φ) ρs Dρ s Dt +(1−φ) ∇· ˙ u , (4) ∂ (( φc l ρ l + ( 1−φ ) c s ρ s) �T ) ∂t +∇·(−∇T)+∇•  clρl�Tql  +∇·((φclρl+(1−φ)csρs)�T˙u)= 0, (5) R ρlφcl DT Dt +∇·(−∇T)+clρlql·∇T= 0, (6) =b+d. (7) ∇· σ +b=0, (8) ε = σ ′ m 3K I+1 2G σ ′ −σ ′ mI  + αT 3 TI , Page 6 of 27 Vidaletal. Geothermal Energy (2022) 10:23 times the linear) thermal expansion coefficient of the medium [θ−1], K=E/3(2−ν) is the bulk modulus [ML−1 T−2], G=E/(2(1+2ν) is the shear modulus [ML−1 T−2] 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: where pl is the liquid pressure [ML−1 T−2]. A summary of the basic relationships for mass and momentum balance equations is given in Table1. 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: In each variable, subscript D refers to the dimensionless variable and subscript with C refers to the characteristic variable. (9) G ∇2˙u+  K−2G 3  ∇(∇·˙u)−  K+2G 3  αT ∂ ∂t ∇T−∂ ∂t ∇pl= 0, (10) t D= t tC ;TD= T−T 0 TC ;plD = p l −p 0 plC ;uD= u uC . Table 1 Basic relationships for mass and momentum balance equations Liquid advective flux (Darcy’s law) q l= −k µl (∇pl+ρlg∇z ) Heat conduction (Fourier’s law) ic=−b∇T Heat dispersion (Fourier’s law) id=−d∇T Liquid density ρl =ρ l 0exp  α lp (p l −p l 0)+α lT (T−T0)  Solid density ρ s=ρs0exp  αsp(p−ps0)+αsT (T−T0)  Thermal retardation coefficient R = φc l ρ l +(1−φ)c s ρ s φclρl Bulk thermal conductivity b=(1−φ)ns+φnl Dispersive conductivity d = clρld | ql| Page 7 of 27 Vidaletal. Geothermal Energy (2022) 10:23 Equations(5) and (9) can be written in terms of the dimensionless variables of (10) as follows: Assuming radial flow ( q l= |Q| 2πrb aq ), Eq.(11) can be written as: The characteristic length Lc 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 tC represents the characteristic time of the advection equation: 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: where Q is the injection flow rate [L3T−1] andbaq is the thickness of the injected aquifer [L]. In Eq.(15), Pecond is the the ratio between heat transfer by advection and by diffusion [−]: We define Pecond in (16) for one well. In a 5-spot scheme, we express Pecond as: Pedisp is the ratio between heat transfer by advection and by dispersion [−]: UT is a dimensionless number which expresses the ratio between the thermal strain generated by the injection of hot water and the initial total strain [−]: (11) DT D DtD −Pe −1 2 ∇ D ·(∇DTD)+|ql|πLcbaq |Q| ·∇ DTD= 0, (12) ∇ 2 D˙uD+2 3 3ν−1 2−ν ∇D(∇D·˙uD)−2 3 ν+3 2−ν UT ∂ ∂tD ∇DTD−2 3 1+2ν 2−ν UP ∂ ∂tD ∇DplD = 0. (13) 2DT D DtD −Pe−1∇D·(∇DTD)+ L C r ·∇ DTD= 0. (14) t C= L 2 cφR 2r | q l| . (15) Pe =|Q|ρlcl 2π baq ( b + d )=  Pe−1 cond +Pe−1 disp  , −1 (16) Pe cond = |Q|ρ l c l 2πb aq  b . (17) Pe cond = q l rρ l c l  b . (18) Pe disp = r dl . Page 8 of 27 Vidaletal. Geothermal Energy (2022) 10:23 UP 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: In (19) and (20), uC is the characteristic displacement: where S =  α lp φ+α sp is the storage coefficient [−] and T0 is the temperature of the medium before injection [θ]. The characteristic temperature TC [θ] is defined as: where Tinj is equal to the temperature of the injected water [θ]. The characteristic liquid pressure pC [ML−1 T−2] is equal to: where p0 is the liquid pressure in the medium previous to the injection [ML−1 T−2] and pinj is the injected liquid pressure [ML−1 T−2] defined as: where g is the gravitational acceleration [LT−2] and h is the variation in the hydraulic head defined by Thiem (1906): where rwell is the well radius [L] and k is the intrinsic permeability [L2]. We evaluate the importance of convection in the HT-ATES system with the dimensionless Rayleigh number ( Ra ). Ra is the ratio between the characteristic time for conductive transport and that for buoyant thermal transport: where �ρl is the fluid density difference [ML−3] and µl is the dynamic viscosity of liquid [ML−1 T−2]. There is a critical Rayleigh number Rac equal to 4π2 for a case of convective flow in a hexahedral domain of porous material saturated with fluid, (Lapwood 1948; Beck 1972). If Ra is higher than Rac , heat convection is the dominant process, else heat conduction is the dominant one. Rac only considers free convection as heat transport process. (19) U T= α T 3TCLc uC . (20) U P= p c L c K uC . (21) uC=p0SLc+T0(αlT φ+αsT )Lc, (22) TC=Tinj −T0, (23) pC=pinj, (24) pinj = ρ lg � h, (25) � h= Q 2π baqk ln r well Lc, (26) R a= tcond tgrav = baq2/  b Rclρlφ  baq/  �ρlkg µ l = baqRclρlφ�ρlkg µlb , Page 9 of 27 Vidaletal. Geothermal Energy (2022) 10:23 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): where Tprod is the temperature of the extracted water [θ], Tinj is the temperature of the injected water and [θ] T0 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 Ra lower than the Rac , 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). 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. (27) η = t 0ρlcl|Qext|(Text −T0)dt  t 0 ρ l c l  Qinj   Tinj −T0  dt , Page 16 of 27 Vidaletal. Geothermal Energy (2022) 10:23 time. The energy recovery factor increases to 0.89 asymptotically for a decade of operation. This value is in agreement with the relation between Ra - η proposed by Gutierrez- Neri etal. (2011). 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 UT and UP 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 UT (90°C) is higher than UP, that is, the effect of temperature dominates over the liquid pressure. UP is equal for the injection and back injection (same flow rates). In the back injection UP dominates over 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 Fig. 11 Energy recovery factor during 10 years of operation Page 17 of 27 Vidaletal. Geothermal Energy (2022) 10:23 UT (50°C). UT and UP decrease with depth due to the increment of the hydrostatic pressure. 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. 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 Fig. 12 UT and UP in the injection and in the back injection in the aquifer. UP is the same in the two different stages Fig. 13 Thermal strain (A), hydraulic strain (B) and total strain (C) in the central well (Caqu) and in the auxiliary well (Aaqu) in the aquifer Page 18 of 27 Vidaletal. Geothermal Energy (2022) 10:23 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. 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 Fig. 14 Vertical displacements with thermal expansion and without thermal expansion on the top of the central well (Csrf) and on the top of the aquifer on the central well (Caqu) Fig. 15 Surface vertical displacements at the end of the injection and at the end of the back injection. Csrf represents the top of the central well and Asrf the top of one auxiliary well Page 19 of 27 Vidaletal. Geothermal Energy (2022) 10:23 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. 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 m2), 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 (2Qi) and half (0.5Qi) that of the base case (Qi = 25 ls−1). 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). 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). 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). Page 20 of 27 Vidaletal. Geothermal Energy (2022) 10:23 Displacements Since we do not modify the temperature of the injected water in any of the sensitivity case, UP is the only strain component that changes. The UP is inversely proportional to the intrinsic permeability (Fig.21A). If we increase the flow rate in the aquifer, we observe more hydraulic strain and UP (Fig.21B). The surface vertical displacements are a consequence of events in the aquifer. Higher values of UP 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. 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) Page 21 of 27 Vidaletal. Geothermal Energy (2022) 10:23 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 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) Page 22 of 27 Vidaletal. Geothermal Energy (2022) 10:23 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π2. Only our model with the highest permeability slightly exceeds the critical Rayleigh number. 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) Page 23 of 27 Vidaletal. Geothermal Energy (2022) 10:23 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 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) Fig. 21 Up 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) 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) Page 24 of 27 Vidaletal. Geothermal Energy (2022) 10:23 fluxes), UP (effect of hydraulic strain) and UT (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. UP quantifies the effect of the hydraulic strain generated by the flow rate injected and UT 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. Abbreviations b Body forces vector [ML−2 T−2] baq Aquifer thickness [L] cl Specific heat of liquid [FLM−1θ−1] cs Specific heat of solid [FLM−1θ−1] d Heat dispersivity [L] dl Longitudinal heat dispersivity [L] D(a) Dt Material derivative with respect to solid/liquid phase velocity E Young’s modulus [ML−1 T−2] G Shear modulus [ML−1 T−2] g Gravitational acceleration [LT−2] h Variation in the hydraulic head [L] ic Heat conduction [Fθ2T−1L−1] id Heat dispersion [Fθ2T−1L−1] K Bulk modulus [ML−1 T−2] k Intrinsic permeability [L2] Lc Characteristic length [L] Pe Dimensionless Peclet number [‑] Pecond Ratio between heat transfer by advection and by diffusion [‑] Pedisp Ratio between heat transfer by advection and by dispersion [‑] plC Characteristic liquid pressure [ML−1 T−2] pinj Injected liquid pressure [ML−1 T−2] p0 Liquid pressure in the medium before the injection [ML−1 T−2] pl Liquid pressure [ML−1 T−2] Q Injection flow rate [L3T−1] 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) Page 25 of 27 Vidaletal. Geothermal Energy (2022) 10:23 ql Liquid advective flux [LT−1] R Thermal retardation coefficient [‑] Ra Rayleigh number [‑] Rac Critical Rayleigh number [‑] r Radial distance from the central well [L] rwell Radius of the well [L] S Storage coefficient [‑] T Temperature [θ] TC Characteristic temperature [θ] Tinj Temperature of the injected water a [θ] T0 Temperature in the medium before the injection [θ] Tprod Temperature of the extracted water [θ] tC Characteristic time [T] UT Ratio between the thermal strain generated by the injection of hot water and the initial total strain [‑] UP Ratio between the hydraulic strain generated by the injection of hot water and the initial total strain [‑] ˙u Solid velocity vector [LT−1] uC Characteristic displacement [L] αlp Liquid compressibility [F−1L2] αlT Volumetric (3 times the linear) thermal expansion coefficient of the liquid [θ−1] αsp Grain compressibility [F−1L2] αsT Volumetric (3 times the linear) thermal expansion coefficient of the solid [θ−1] αT Volumetric (3 times the linear) thermal expansion coefficient of the medium [θ−1] ε Elastic strain tensor [‑] η Energy recovery factor [‑]  Thermal conductivity [FT−1θ1] b Bulk thermal conductivity [FT−1θ1] d Dispersive conductivity [FT−1θ1] µl Dynamic viscosity of liquid [ML−1 T−1] ρl Liquid density [ML−3] ρl0 Liquid density for a reference temperature and pressure [ML−3] ρs Solid density [ML−3] ρs0 Solid density for a reference temperature and pressure [ML−3] σ Stress tensor [ML−1 T−2] σ ′ Effective stress tensor [ML−1 T−2] σ ′ m Mean effective stress [ML−1 T−2] ν Poisson ratio [‑] φ Porosity [L3L−3 Acknowledgements We would like to thank Peter Meier and Geo‑Energy Suisse AG for letting us participate in the Bern pilot project. We acknowledge the financial support received from the ERANET project HEATSTORE (170153‑4401). This project has been subsidized through the ERANET cofund GEOTHERMICA (Project n. 731117), from the European Commission, RVO (the Netherlands), DETEC (Switzerland), FZJ‑PTJ (Germany), ADEME (France), EUDP (Denmark), Rannis (Iceland), VEA (Belgium), FRCT (Portugal), and MINECO (Spain). We wish to thank the Department of Research and Universities of the Generalitat de Catalunya by supporting RV with a grant (2021 FI_B 00940). Author contributions RV carried out the data acquisition and processing, participated in research analysis, interpretation of results and drafting of the manuscript. SO and MWS participated in the research analysis, interpretation of results and drafting of the manu‑ script. FDM participated in the revision of manuscript. All authors read and approved the final manuscript. Funding We acknowledge the financial support received from the ERANET project HEATSTORE (170153–4401). This project has been subsidized through the ERANET cofund GEOTHERMICA (Project n. 731117), from the European Commission, RVO (the Netherlands), DETEC (Switzerland), FZJ‑PTJ (Germany), ADEME (France), EUDP (Denmark), Rannis (Iceland), VEA (Belgium), FRCT (Portugal), and MINECO (Spain). We wish to thank the Department of Research and Universities of the Generalitat de Catalunya by supporting RV with a grant (2021 FI_B 00940). Availability of data and materials They are included in the publication. Declarations Ethics approval and consent to participate This research does not involve human subjects, human material or human data.