Risk of surface movements and reservoir deformation for high-temperature aquifer thermal energy storage (HT-ATES)

Abstract

High-temperature aquifer thermal energy storage (HT-ATES) systems are designed for seasonal storage of large amounts of thermal energy to meet the demand of industrial processes or district heating systems at high temperatures (> 100 °C). The resulting high injection temperatures or pressures induce thermo- and poroelastic stress changes around the injection well. This study estimates the impact of stress changes in the reservoir on ground surface deformation and evaluates the corresponding risk. Using a simplified coupled thermo-hydraulic-mechanical (THM) model of the planned DeepStor demonstrator in the depleted Leopoldshafen oil field (Upper Rhine Graben, Germany), we show that reservoir heating is associated with stress changes of up to 6 MPa, which can cause vertical displacements at reservoir depth in the order of 10^–3 m in the immediate vicinity of the hot injection well. Both the stress changes and the resulting displacements in the reservoir are dominated by thermoelasticity, which is responsible for up to 90% of the latter. Uplift at the surface, on the contrary, is primarily controlled by poroelasticity with by two orders of magnitude attenuated displacements of << 10^–3 m. Our calculations further show that the reservoir depth, elastic modulus, and injection/production rates are the dominant controlling parameters for the uplift, showing variations of up to two order of magnitudes between shallower reservoirs with low elastic moduli and deeper and more competent reservoirs. In addition, our findings demonstrate that the cyclic operation of HT-ATES systems reduces the potential for uplift compared to the continuous injection and production of conventional geothermal doublets, hydrocarbon production, or CO₂ storage. Consequently, at realistic production and injection rates and targeting reservoirs at depths of at least several hundred meters, the risk of ground surface movement associated with HT-ATES operations in depleted oil fields in, e.g., the Upper Rhine Graben is negligible.

Introduction

In the power sector, energy storage is gaining importance due to volatility of renewable energy sources (Lee 2013; REN21 2019). In the heat sector, seasonal variation in demand leads to large quantities of excess heat in summer, which can only be stored in the underground (Lee 2013). Among other underground storage technologies (Yang et al. 2021), near-surface aquifer thermal energy storage (ATES) systems are widely used in the low-temperature sector. As of 2018, over 2800 systems were in operation worldwide, with an increasing trend (Fleuchaus et al. 2018). Process heat or water-based district heating with an inlet temperature of 80–130 °C require corresponding temperature levels and storage capacities for significantly larger amounts of thermal energy (Wesselink et al. 2018). So far, only a few HT-ATES systems are in operation, e.g., in Berlin and Neubrandenburg in northeast Germany (Holstenkamp et al. 2017) or in a demonstration phase, e.g., the Middenmeer project in the Netherlands (Dinkelman et al. 2022; Oerlemans et al. 2022).

Hydrocarbon reservoirs have favorable storage properties and ambient temperatures, as evidenced by decades of exploration, exploitation, and research in different geological settings, such as the Upper Rhine Graben (URG) in France, Germany, and Switzerland (Böcker 2015; Stricker et al. 2020), clastic aquifers in the Netherlands (van Wees et al. 2017), or the Geneva basin in Switzerland (Moscariello 2019). One challenge associated with their utilization for seasonal storage is the presence of residual hydrocarbon content. Insights can be drawn from the geothermal utilization of depleted reservoirs (Liu et al. 2015; Wang et al. 2016; Sui et al. 2019). Furthermore, compaction-induced subsidence is a well-documented phenomenon in hydrocarbon production (Geertsma 1973; Nagel 2001) and must also be considered for HT-ATES. In oil fields of the Upper Rhine Graben that are comparable to the targeted Leopoldshafen field, such as the Landau field, subsidence reached up to 7 mm a⁻¹ during production between the 1950s and the 1990s, accumulating to about 17 cm (Fuhrmann 2016). In addition, the risk of the reactivation of nearby faults and associated induced seismicity, which has been frequently observed for other utilizations of the subsurface, such as geothermal utilization or hydrocarbon production (Zang et al. 2014; Bourne et al. 2015), has to be addressed as well.

Concepts for ATES range from single-well systems (e.g., Jeon et al. 2015) to spatial patterns of multiple injection and production wells (e.g., Sommer et al. 2015). However, most ATES systems use the push–pull concept, in which excess heat is injected into the hot well of a doublet in summer and cold fluid is simultaneously produced from a second well. In winter, the system is operated in reverse (Kim et al. 2010; Jin et al. 2021). This concept is also used in the DeepStor demonstrator at the Karlsruhe Institute of Technology (Banks et al. 2021; Rudolph et al. 2022).

While single-well projects and multiple well patterns are prone to temporal mass unbalancing, which can cause differential vertical deformation, doublet operations ensure overall mass balancing at any time. Nevertheless, local mass imbalances between the injection and production areas remain and are investigated in this study. Previous research mainly focused on (i) uplift caused by poroelasticity for a single injection period (Birdsell and Saar 2020), (ii) the influence of poro- and thermoelasticity on the storage efficiency (Jin et al. 2020, 2022), and (iii) the sensitivity of hydraulic parameters (e.g., the reservoir permeability) on the resulting uplift (Vidal et al. 2022).

Although these studies provide a good overview of the impact from seasonal ATES on ground surface displacements, they lack a fundamental understanding of the underlying processes and the influence of poro- and thermoelastic stress variations on the reservoir, nearby faults and resulting uplift. This is especially true for the sensitivity of other parameters, such as mechanical properties (e.g., elastic moduli), thermal properties (e.g., thermal conductivity), or the stress field. In addition, the geological and operational differences between ATES systems and other subsurface utilizations, such as deep geothermal and hydrocarbon operations, must be considered, as ATES systems are typically located at shallower depths and involve less competent rocks (Lee 2010; Fleuchaus et al. 2018).

In this study, we use a thermo-hydraulic–mechanical (THM) model to investigate the impact of seasonal HT-ATES on the subsurface stress distribution and the resulting deformation at reservoir depth and the ground surface. We particularly focus on understanding the combined poro- and thermoelastic processes associated with the operation of HT-ATES systems in shallower depths and less competent rocks than comparable deep geothermal systems. Our objective is to perform a risk assessment of uplift for the DeepStor demonstrator, which includes a sensitivity analysis of reservoir and operational parameters.

Numerical modeling

This study uses the open-source code TIGER (THMC sImulator for GEoscientific Research; Gholami Korzani et al. 2020; Egert et al. 2022), which is based on the object-oriented framework MOOSE (Permann et al. 2020; Lindsay et al. 2022). TIGER has been successfully applied to simulate, for example, fluid flow in fractures (Egert et al. 2021), solute transport in fractured EGS reservoirs (Egert et al. 2020), heat storage (Stricker et al. 2020), and coupled mechanical processes (Egert et al. 2022). In Additional file 1 to this study, the results of the validation of poroelastic stress calculations in TIGER are provided. Therein, stress calculations in TIGER were compared with the analytical solution of Rudnicki (1986), following the procedure outlined in Altmann et al. (2010). For further insights into the implementation of the underlying continuum mechanics framework in MOOSE (Tensor Mechanics Module), we refer to the comprehensive documentation by INL (2023). This study extends previous heat storage simulations by considering coupled mechanical processes.

Governing equations

The equation for mass transport used to estimate the pore pressure, p, is given by the mass balance (Eq. 1) along with the Darcy velocity, q (Eq. 2; Bear 1972):

$${{S}}_{{{m}}} \frac{{{\mathbf{\partial }}{{p}}}}{{{\mathbf{\partial }}{{t}}}} + {{\mathbf{\nabla }}}{\varvec{q}} = {{Q}}$$

(1)

$${\varvec{q}} = \frac{{\varvec{k}}}{{{\mu}}}\left( { - {{\mathbf{\nabla }}}{{p}} + {{\rho}}_{{{f}}} {\varvec{g}}} \right)$$

(2)

S_m is the mixture-specific storage coefficient of the medium, t is the time, Q is the source/sink term for injection and production, k is the permeability tensor, µ and ρ_f are the fluid dynamic viscosity and density, respectively, and g is the gravitational acceleration.

The solid and liquid phases are assumed to be in local thermodynamic equilibrium. Heat transport used to estimate temperature can be mathematically expressed using the advection–diffusion equation (Eq. 3):

$${\rho c}_{{{p}}} \frac{{{\mathbf{\partial }}{{T}}}}{{{\mathbf{\partial }}{{t}}}} - {\varvec{\lambda}}_{{{b}}} {\mathbf{\nabla }}^{2} {{T}} + \left( {{\rho c}_{{{p}}} } \right)_{{{f}}} {\varvec{q}}{\mathbf{\nabla }}{{T}} = 0$$

(3)

T is the temperature of the porous medium. ρc_p and λ_b are the heat capacity and thermal conductivity of the mixture, respectively. (ρc_p)_f represents the heat capacity of the fluid.

The deformation of a fully saturated porous medium is described by the momentum balance equation as follows for the form of effective stress (Eq. 4; Jaeger et al. 2007):

$${\mathbf{\nabla }} \cdot \left( {\varvec{\sigma^{\prime}} - {\alpha p}} \right) + {{\rho}}_{{{s}}} {\varvec{g}} = 0$$

(4)

σʹ is the effective stress tensor, α the Biot’s coefficient, and ${\rho }_{s}$ the solid rock density.

The constitutive law for stress–strain behaviour links the displacement vector u, the primary variable solved for the deformation of the porous medium, to the effective stress tensor σʹ:

$${{{\Delta}}}\varvec{\sigma^{\prime}} = {\mathbf{\mathbb{C}}}\left( {{{{\Delta}}}{\varvec{\epsilon}} - \frac{1}{3}{{\beta}}_{{{d}}} {{{\Delta}}}{{T}}} \right)$$

(5)

$${\varvec{\epsilon}} = \frac{1}{2}\left( {{\mathbf{\nabla }}{\varvec{u}} + \left( {{\mathbf{\nabla }}{\varvec{u}}} \right)^{{{T}}} } \right)$$

(6)

${\mathbf{\mathbb{C}}}$ is the elastic material tensor, ${\varvec{\epsilon}}$ the strain, β_d the volumetric drained thermal expansion coefficient, and ΔT the temperature change. Here, we consider only isotropic, non-isothermal elastic deformation. Therefore, linear elasticity can be fully described by the generalized Hooke’s law:

$${\mathbf{\mathbb{C}}} = \varvec{\lambda}{ \delta }_{{{{ij}}}} {{\delta}}_{{{{kl}}}} + 2\varvec{\mu}{ \delta }_{{{{ik}}}} {{\delta}}_{{{{jl}}}}$$

(7)

δ is the Kronecker delta, µ the shear modulus, and λ the Lamé constant.

Numerical model

As an example, we developed a simplified 3-D model of the potential DeepStor site to simulate coupled mechanical processes for high-temperature heat storage. The model extends over an area of 12.5 × 12.5 km from the ground surface to a depth of 2.2 km. It includes a 10-m-thick reservoir at a depth of about 1200 m, over- and underlain by impermeable confining layers (Fig. 1). The reservoir dips 5° to the East and resembles a typical setting within a graben block of the URG. The lateral extension of the model was chosen to avoid boundary effects, especially with regard to pressure propagation. Two vertical wells are located in the center of the model with a lateral spacing of 300 m to avoid thermal interference between them (e.g., Sommer et al. 2015). The unstructured mesh of the model consists of tetrahedral elements as well as 1D line elements to implement the borehole doublet. It was created using the Gmsh software (Geuzaine and Remacle 2009). Element sizes vary from about 3 m along the two wells to 250 m at the upper and lower boundaries of the model. Further refinement was performed, where the highest gradients in the pressure, temperature, and displacement fields are expected, particularly around the two boreholes and above and below the near-borehole-areas in the caprock. A mesh sensitivity was performed to avoid any mesh dependency in the results. In total, the model consists of 171,346 nodes connected by 1,013,332 elements.

Fig. 1

Simplified 3D mesh of the potential DeepStor demonstrator. The model consists of a 10-m-thick reservoir (red) dipping 5° to the east, embedded into two impermeable confining layers (blue). The zoom shows the refined section of the model around the hot and cold wells (HW and CW, respectively)

Hydrostatic pore pressure is applied as a Dirichlet boundary condition at the top and bottom boundaries of the model, with the appropriate initial conditions applied to the entire model. No-flux Neumann boundary conditions are applied to the lateral boundaries of the model. Injection and production flow are implemented using time-dependent mass flux functions at the intersection of the top of the reservoir and the wells (corresponding to the beginning of the open hole section). These time-dependent functions represent simplified approximations of a real pumping operation, where an instantaneous reversal of the pumping direction is assumed at the end of each cycle. The temperature distribution within the model is based on a geothermal gradient of 50 K km⁻¹ (Agemar et al. 2012) and is achieved by setting Dirichlet boundary conditions at the top and the bottom boundaries of the model and an associated initial condition applied to the entire model. To implement fluid injection at a given temperature into the two wells, nodal Dirichlet boundary conditions on top of the open-hole section are integrated and activated during the injection period of each well. No prior stress state is applied to the model, and the bottom and side boundaries of the model are fixed in normal directions by defining displacement-free Dirichlet boundary conditions (rollers). Gravity affects the pore pressure and effectives stress within the model.

A reference model (hereafter referred to as “reference case”) was developed to investigate the general thermo-hydraulic–mechanical coupled behavior of the potential DeepStor demonstrator. The parameterization (Table 1) is based on average reservoir properties of (former) oil fields in the URG (Stricker et al. 2020). The seasonal storage operation consists of 6 months of continuous fluid injection at 140 °C into the hot well during summer (with the cold well acting as producer) and an inverted mode during winter, when water at 70 °C is injected in the cold well (with the hot well acting as producer). This seasonal storage operation is performed over a period of 10 years. In the following, the operation in summer will be referred to as “charging period”, whereas the operation in winter will be referred to as “discharge period”.

Table 1 Parametrization of the reference model of the DeepStor site based on average reservoir properties of (former) oil fields in URG, particularly the Leopoldshafen field (Stricker et al. 2020; Banks et al. 2021)

Results

Mechanical interaction of heat storage

To improve the understanding of the coupled mechanical impact of one charging period (6 months), Fig. 2 illustrates the three-dimensional distribution of the following parameters: (a) pressure, (b) temperature, (c) vertical stress, and (d) vertical strain. As the area of interest is primarily limited to the vicinity of the wells, only an excerpt with an extent of 1000 m × 1000 m × 300 m is shown, depicting the bottom surface of the reservoir in the foreground and a vertical section in the center of the model.

Fig. 2

Perturbation of pressure (a), temperature (b), vertical stress (c), and vertical strain (d) in the center of the model near the well doublet (HW: hot well; CW: cold well) after one charging period (6 months). The figures have an extent of 1000 m × 1000 m × 300 m. a Additionally shows the location of the displayed excerpt within the full model

The pressure perturbation reaches up to 3.5 MPa at the injection and production wells, affecting a relatively large volume (Fig. 2a). It attenuates to 0.1 MPa at a lateral distance of 800 m. Above and below the reservoir in the impermeable bedrock or caprock, this pressure perturbation nearly dissipates within 25 m (Fig. 3a). In contrast, the temperature changes (Fig. 2b) are confined to the vicinity of the injection well. Overall, a cylindrical volume in the reservoir is subject to temperature changes > 1 K within about 50 m around the hot well. Similar to the pressure perturbations, the thermal perturbations reach zero about 25 m above and below the reservoir (Fig. 3b).

Fig. 3

Distribution of the perturbation of pressure and the three principal stress components (a) as well as the temperature and the three principal strain components (b) above and below the reservoir (marked in grey)

The vertical stress changes around the hot well are the superposition of injection-related poroelastic stress decrease and the thermoelastic stress increase. Figure 2c illustrates this variation as a spatially concentrated increase in the effective stress of up to ca. 1.9 MPa close to the hot well. A cylindrical volume within a distance of 50 m to the well is affected within the reservoir, resembling the observed temperature perturbations. Around this affected volume, the poroelastic stress reductions form a halo effect with extensional stresses of up to − 1.2 MPa. At greater distances to the well, only the poroelastic stress decrease is remaining. In contrast, the vertical effective stress changes around the cold well are only affected by the production-related stress increase due to a decrease in pore pressure, reaching 2.3 MPa directly at the cold well. Above and below the reservoir, the vertical stress changes dissipate quickly, but do not reach zero within the same distance as pressure and temperature (Fig. 3a). The horizontal stress changes are about four times larger in magnitude than the vertical stress changes and become negative directly above and below the reservoir before they also trend towards zero.

The superposition of poro- and thermoelasticity is also visible in the vertical strain (Figs. 2d and 3b), causing a vertical strain of about 10^–3 due to the injection of hot water (pore pressure and temperature increase). At the cold well, negative strains of up to − 2 × 10^–4 occur due to the production-related pressure reduction. However, the thermal expansion within the reservoir and the associated positive strain leads to negative strain directly above the reservoir due to the compression of the rock. At larger vertical distances from the reservoir, the vertical strain trends towards zero. The horizontal strain has a lower magnitude compared to the vertical strain and also dissipates towards zero over short vertical distances above and below the reservoir. The strain around the wells also has an impact on the porosity and permeability distribution. By applying the exponential strain–porosity relation developed in Chen et al. (2009), a relative porosity reduction of up to 0.5% is determined around the hot well due to the occurring solid expansion, equaling an absolute change in the porosity of 0.00075 for a reservoir porosity of 0.15. These changes in porosity also have an influence on the permeability (e.g., Carman 1937), translating in relative changes of approximately 1% (equaling absolute changes of 2.5 * 10^–16 m² for a reservoir permeability of 2.5 * 10^–14 m². However, these small changes in reservoir porosity and permeability have a negligible impact on the storage operation itself as well and will thus not further be considered.

To further analyze the mechanical response after one charging period, the vertical displacements (u_z) at the top of the reservoir and the ground surface were investigated. At the top of the reservoir, an injection-related u_z of up to 5.1 mm occurs, sharply limited to a radius of about 50 m around the hot well and reduces towards the model boundaries (Fig. 4a). Only 1% of the maximum u_z remains at a distance of about 800 m from the hot well. At the ground surface, only about 0.06 mm uplift (u_z0) remains, but it shows an impact on a larger area, with the uplift trending towards zero at distances greater than 3 km to the wells (Fig. 4b). At the top of the reservoir around the cold well, a negative u_z of about − 0.6 mm occurs. This negative u_z translates into an u_z0 of ca. − 0.05 mm at the ground surface, comparable to the u_z0 above the hot well in both in magnitude and in extent.

Fig. 4

Distribution of the vertical displacement (u_z and u_z0, respectively) after one charging period at a the top of the reservoir (1195 m depth) and b the ground surface. Comparison of the displacement along the x-axis (at y = 0 m; dashed red lines in a and b) considering only poroelasticity (blue) and combining poro- and thermoelasticity (black) at the top of the reservoir (c) and the top of the model (ground surface; d). The different x-axis extents at reservoir depth (a and c) and the ground surface (b and d) correspond to the different scales of the involved physical processes at different depths

The vertical displacements at the top of the reservoir and the ground surface are caused by the superposition of poro- and thermoelasticity, reflecting the stress and strain changes described above. To distinguish between the respective influences of poro- and thermoelasticity, a second purely poroelastic model was simulated. The comparison of this model with the reference case shows that 89% of u_z at the top of the reservoir (about 4.6 mm) is related to thermoelasticity. However, its influence is limited to a radius of about 50 m around the hot well (difference between the black and the blue line in Fig. 4c), resembling the thermal front of the hot water injection after 6 months. In contrast, poroelasticity—although it influences a larger rock volume—contributes only about 0.6 mm (or 11%) of the total u_z at the top of the reservoir. Thus, at greater distances to the hot water injection well—and around the production well, where no temperature changes take place—u_z at reservoir depth is purely caused by poroelasticity. In contrast to the reservoir, the u_z0 at the ground surface is strongly dominated by the poroelastic component (contributing up to ca. 92% of u_z0), whereas the thermoelastic component is rather negligible (Fig. 4d).

This apparent contradiction between dominating thermoelasticity at reservoir depth and dominating poroelasticity at the ground surface can be explained by the different attenuation of the thermo- and poroelastic components of u_z from the top of the reservoir to the ground surface. Figure 5 illustrates that the higher thermoelastic component of u_z (red line) is caused by higher thermal vertical strain (${\epsilon}$_zz) above the reservoir due to strong thermal expansion of the rock matrix until about four meters above the top of the reservoir (dashed red line), where the highest u_z occurs. At greater vertical distances to the top of the reservoir, the rock matrix is compressed (negative ${\epsilon}$_zz), subsequently resulting in a strong attenuation of the thermoelastic u_z.

Fig. 5

Vertical displacement (solid lines; left y-axis) and strain (dashed lines; right y-axis) variation from reservoir top at 1195 m to a depth of 1000 m

In contrast, the rock (matrix and pore space) above the reservoir undergoes lower poroelastic ${\epsilon}$_zz (dashed blue line in Fig. 5), resulting in lower u_z (blue line). This subsequently leads to weaker compression of the overlying rock (negative ${\epsilon}$_zz) and thus to less pronounced attenuation of the poroelastic component of u_z. In addition, the maximum of the poroelastic component of u_z (equal to the depth, where ${\epsilon}$_zz = zero) is located further from the top of the reservoir than the maximum of the thermoelastic component (about 15 m compared to about 4 m). Because of the weaker attenuation of the poroelastic component of u_z, the poro- and thermoelastic components of u_z are equal at about 1100 m (about 90 m above the top of the reservoir), where the red and blue curve intersect. This leads to the poroelastic component being higher than the thermoelastic component in depths shallower than about 1100 m and accordingly dominating at the ground surface.

Temporal evolution

In the following, we address the influence of seasonal charging and discharge cycles on stress changes within the reservoir, considering the different contributions of poro- and thermoelasticity. Figure 6a shows the vertical stress changes (Δσ_zz) after 6 months of injection into the hot well and production from the cold well (one charging period). The figure (as well as the other parts of Fig. 6) shows an 1000 m × 1000 m excerpt at the top of the reservoir, where the most relevant changes in Δσ_zz take place. At the hot well, the superposition of poro- and thermoelastic stress changes leads to an increase in Δσ_zz of 1.8 MPa. Around the rock volume affected by the temperature change (sharply concentrated around the hot well with a radius of about 30 m), a halo effect of poroelastic stress reduction occurs, comprising a reduction in Δσ_zz of 1.1 MPa, exceeding the expected reduction in Δσ_zz for pure poroelasticity by 0.3 MPa. At greater distances from the hot well, only the pressure-induced reduction in Δσ_zz remains and slowly dissipates towards the model boundaries. As a pressure-related reduction in the effective stress of about 0.02 MPa still occurs at a lateral distance of 1500 m, the spatial influence of poroelasticity exceeds that of thermoelasticity by nearly two orders of magnitude. At the cold well, in contrast, the production of water leads to a pressure-related increase in Δσ_zz of ca. 2.0 MPa.

Fig. 6

Distribution of stress changes at the reservoir top around the two wells (CW: cold well; HW: hot well) after a one charging period, b one full charging and discharge cycle, c the tenth charging period (i.e., after 9.5 years), and d the full simulation duration of ten charging and discharge cycles

After a complete charging and discharge cycle (i.e., after 1 year), at the cold well an injection-related decrease in Δσ_zz of − 2.0 MPa occurs, equal to the production-related increase in Δσ_zz (Fig. 6b). At the hot well, the superposition of pressure reduction and residual heat leads to an increase in Δσ_zz of about 3.9 MPa, significantly exceeding the reduction in Δσ_zz around the cold well. At a distance of 35 m, still a stress increase of about 1 MPa occurs, representing the heat remaining in the reservoir after the discharge period.

After the tenth charging period (i.e., after 9.5 years), the heating of the reservoir over time leads to an increase in Δσ_zz from 1.8 MPa (after the first charging period) to 2.3 MPa, and the radius of the affected rock volume increases from 30 to 45 m. The stress reduction halo around the thermoelastically affected rock mass moves further away from the hot well (70 m compared to 45 m) and its absolute reductions in Δσ_zz compared to pure poroelasticity increase from 0.3 MPa to 0.35 MPa. No changes occur at the cold well compared to the first cycle, as it is not influenced by thermoelasticity.

After the full simulation duration of ten charging and discharge cycles, a strong effect of reservoir heating is visible. The superposition of the production-related increase in effective stress and the stress increase due to thermal expansion of the rock matrix leads to a total increase in Δσ_zz of 5.5 MPa (compared to 3.9 MPa after the first year; Fig. 6d). In addition, a growing rock mass is affected by this increase in Δσ_zz, e.g., the distance in which a stress increase of 1 MPa can be observed increases from about 35 m to about 50 m. At the cold well, no changes occur compared to the situation after 1 year.

The changes in Δσ_zz over time at the top of the reservoir caused by cyclic heat storage also impact u_z. Fig. 7a illustrates the response of u_z to the cyclic changes of Δσ_zz at the upper end of the open hole section of the hot well. The poroelastic component (blue curve in Fig. 7a) stays constant over time, because the increase and decrease in pore pressure due to injection and production do not change, resulting in equal changes in (effective) stress and strain and thus also in u_z. However, this does not apply to the thermoelastic component of u_z (red curve in Fig. 7a). It can be observed that u_z at the end of each charging period (local maxima of u_z) decreases from 4.6 mm after the first period to 4.2 mm after the tenth period. In contrast, the remaining u_z after each subsequent discharge period (local minima of u_z) increases from 1.3 mm after the first period to 2.1 mm after the tenth period.

Fig. 7

a Development of the vertical displacement at the intersection between the top of the reservoir and the hot well over time. Blue curve: poroelastic component, red curve: thermoelastic component. b Distribution of the vertical displacement (black curves) and strain (orange curves) around the reservoir (highlighted area in light grey) over time (for combined poro- and thermoelasticity)

The decrease of u_z at the end of each charging period over time can be explained by the gradual heating of the reservoir. With each charging period, this heating causes the temperature front to propagate further into the cap rock above and below. As a result, the temperature gradient between the reservoir boundaries and the cap rock above and below decreases. This leads to lower thermal strains (Fig. 7b) and consequently reduced u_z at the upper end of the open hole section of the hot well over time. The increase in u_z after each discharge period can also be explained by the gradual heating of the reservoir due to diffusive heat losses around the hot well, as water with a temperature exceeding the ambient reservoir temperature is stored (e.g.,Kim et al. 2010; Stricker et al. 2020). Figure 7b further shows that the temperature front propagating away from the reservoir (and thus also the moving thermal strain front) leads to a flattening of the strain rate (orange curves in Fig. 7b) over time. This causes an accumulation of strain and thus an increasing maximum of u_z (black curves in Fig. 7b). In addition, the maximum of u_z moves further away from the reservoir as a consequence of the moving thermal strain front.

Risk assessment of uplift at the ground surface

The results presented in chapter 3.1 have shown that the uplift at the ground surface (u_z0) is predominantly influenced by poroelasticity and subordinately by thermoelasticity. Therefore, parameters that influence either the hydraulics in the reservoir (such as the reservoir permeability or the injection/production flow rate) or the mechanical response of the rock (such as the Young’s modulus or the Poisson’s ratio) are very likely to have the strongest influence on u_z0. In the following, we focus on these parameters, complemented by the injection temperature difference, the thermal expansion coefficient, and the reservoir depth (Table 2).

Table 2 Selected ranges of identified geological and operational parameters to determine their influence on the uplift at the ground surface (u_z0)

Figure 8 shows the results of the performed parameter sensitivity analysis. The strongest influence is exerted by the reservoir depth (dashed black line), with an increase in u_z0 of about one order of magnitude for a reduction in reservoir depth from 1600 to 400 m. A decrease in the Young’s modulus (green line) from 25 to 2.5 GPa also leads to an increase in u_z0 of nearly one order of magnitude. The injection/production flow rate (red line) has a linear influence on u_z0, with a fivefold increase in u_z0 for an increase in the flow rate from 1 to 5 Ls⁻¹. The large influence of these three parameters is related to not only their relative impact (slope of the curves) but also the large potential variation of the parameters themselves. The other four parameters do not exert a significant influence on u_z0.

Fig. 8

Sensitivity analysis on the parameter influence on u_z0 after the first charging period at the hot well. The shown changes refer to the reference case

Due to the strong non-linearity of the relationship between both reservoir depth and Young’s modulus with u_z0, high values of u_z0 (superimposed by the linear influence of the flow rate) can be expected for shallow reservoirs with low Young’s moduli. To investigate this behavior more thoroughly, we simulated a range of models that vary the reservoir depth between 400 and 1600 m (shallow and deep oil reservoir in the URG, respectively; Stricker et al. 2020) and the Young’s modulus ranging between 2.5 GPa (weak claystone; Marschall and Giger 2014; Jahn et al. 2016) and 25 GPa (strong sandstone; Egert et al. 2018). All other parameters used in these simulations were the same as in the reference case. Figure 9 illustrates the strong dependency of u_z0 on the two parameters with variations between 0.02 mm for a deep reservoir with a high Young’s modulus (comparable to the parametrization of a (medium) deep geothermal system) and 1.18 mm for a shallow reservoir with a low Young’s modulus, comprising a variation in u_z0 of nearly two orders of magnitude. The curved shapes of the isolines in Fig. 9 further illustrate the non-linear relationship between the two varied parameter and u_z0.

Fig. 9

Contours of the maximum uplift after one charging period for different reservoir depths and Young’s moduli (E). The grey asterisks show the simulations, intermediate results are interpolated. The red asterisk represents the reference case

The strong variations in u_z0 observed in HT-ATES are caused by the superposition of two processes. First, the poroelastic component of u_z0 decreases with depth, because the dominant poroelastic component of u_z (ca. 0.6 mm for all reservoir depths) is attenuated towards the ground surface as the path length increases. This finding aligns with earlier results by Birdsell and Saar (2020), but contradicts the results of Vidal et al. (2022), who simulated counterintuitively increasing vertical displacements towards the ground surface due to their choice of boundary conditions. This is particularly evident in the employment of zero vertical displacement boundary conditions at the lateral boundaries of the model by Vidal et al. (2022) in addition to rollers only [as used in this study and also, for instance, in Birdsell and Saar (2020)]. In addition, the limited lateral extent of the model (200 m in contrast to approximately 6 km in this study) may influence the results. An optimized choice thereof may lead to different results in alignment with this study of the results of Birdsell and Saar (2020). Second, an increase in the Young’s modulus leads to a lower u_z in the reservoir, varying between 0.2 mm for a Young’s modulus of 25 GPa and 2.9 mm for a Young’s modulus of 2.5 GPa (both for the depth of the reference case of 1200 m). This variation in u_z at reservoir depth propagates towards the ground surface and contributes to the described variation of u_z0. Consequently, this means that heat storage operations in rather shallow and mechanically incopentent rocks pose much smaller risk for ground surface deformation than deep geothermal operations in more competent reservoir rocks, such as granite or sandstone (Heimlich et al. 2015; Békési et al. 2019; Frey et al. 2022).

Although the thermoelastic component of u_z0 increases by one order of magnitude from 0.003 to 0.025 mm for a reduction of the reservoir depth from 1600 to 400 m, this increase remains negligible compared to the increase in the poroelastic component from 0.035 to 0.30 mm. This indicates that the contribution of the poroelastic component on u_z0 is one order of magnitude higher than the thermoelastic component irrespective of the reservoir depth, matching previous observations of Vidal et al. (2022). Therefore, it can be concluded that the high injection temperatures of HT-ATES systems do not pose a significant additional risk factor regarding uplift.

Ground surface movements around a seasonal heat storage site may be monitored using the same methods as those used for oil production or geothermal utilization, such as global navigation satellite systems (GNSS) or interferometric satellite radar data (InSAR) (Heimlich et al. 2015; Fuhrmann 2016). However, even the highest simulation results of u_z0 (ca. 1 mm) for very shallow reservoirs with low Young’s moduli are approximately one order of magnitude smaller than observations in the vicinity of the Landau oilfield in the URG, where ground surface movements of up to 7 mm.a⁻¹ occurred and accumulated over several decades of production to ca. 17 cm (Fuhrmann 2016). This is because, unlike the continuous injection and/or production of deep geothermal, hydrocarbon exploitation, and CO₂ storage, seasonal heat storage involves cyclical injection and production, which prevents the accumulation of uplift over time.

In summary, the operation of HT-ATES systems in depleted oil reservoirs in the Upper Rhine Graben is expected to only have a minor influence on u_z0, with a very low risk for damage to structures at the ground surface, such as buildings. More significant u_z0 of more than 1 mm (albeit with very low differential displacements) can only occur for very shallow and mechanically weak reservoir (and surrounding) rocks, and can be avoided by selecting a suitable reservoir. In addition, the high sensitivity of u_z0 to injection and production flow rates does not pose a significant risk to HT-ATES operations, as the flow rates in these systems are usually much smaller (i.e., < 5 Ls⁻¹) than in deep geothermal systems with up to > 100 Ls⁻¹ (Evans et al. 2012).

Conclusions

As part of a risk assessment, we used coupled thermo-hydraulic–mechanical modeling to investigate the (geo-) mechanical sensitivity of HT-ATES systems to poro- and thermoelastic stress development. Our modeling results show that the displacements caused by the cyclic heat storage operation occur close to the hot well at reservoir depth, reaching up to 6 mm. Towards the ground surface, the displacements are strongly attenuated by two orders of magnitude. While the displacements in the reservoir are dominated by small-scale thermoelasticity, the related uplift and subsidence at the ground surface are predominantly related to poroelasticity due to the stronger attenuation of the thermoelastic component.

The primary factors influencing ground surface deformation in HT-ATES systems in shallow depleted oil reservoirs in settings, such as contiental rift systems (e.g., the Upper Rhine Graben; URG) are reservoir depth, Young’s modulus, and injection/production flow rates. In contrast to deep geothermal reservoirs, shallow, soft reservoirs with high injection and production flow rates are likely to cause greater ground surface deformation. Therefore, the impact of ground surface deformation can be reduced by choosing suitable reservoir conditions and operational framework. For example, HT-ATES systems typically use lower flow rates than deep geothermal systems. In addition, the cyclic operation of an HT-ATES has a smaller impact on ground surface deformation than typical continental shallow oil production or CO₂ storage in comparable depths as no significant accumulation of uplift or subsidence occurs over time. Ground surface monitoring systems can be used to investigate ground movements at a particular storage site. The highest impact is to be expected in the area above the boreholes, allowing for a relatively simple monitoring layout.

Former oil reservoirs in the URG, such as the planned scientific demonstrator project DeepStor, are considered to be ideal for use as HT-ATES systems. They are expected to be too shallow to pose a high risk of induced seismicity (due to limited hydraulically connection to larger critically stressed faults) and simultaneously deep enough to avoid a high risk of uplift or subsidence. However, future investigations of HT-ATES systems should aim to assess their specific risk of fault reactivation and induced seismicity due to the cyclically stress changes. Beyond purely elastic modeling, more sophisticated modeling approaches, such as damage or phase-field modeling may be considered, particularly to analyze the influence on changes in porosity and permeability. Further studies may also investigate the risk of damaging the cementations of the wells due to the cyclic mechanical loading.