Predicting Thermal Performance of Aquifer Thermal Energy Storage Systems in Depleted Clastic Hydrocarbon Reservoir via Machine Learning: Case Study from Hungary

1. Introduction

The end-of-life management of oil and gas wells—including plugging, abandonment, and site remediation—poses significant economic and environmental challenges. In the U.S., median decommissioning costs average USD 76,000 but can exceed USD 1 million for deep wells. In Hungary, abandonment costs typically range from 50 million to several hundred million HUF, depending on well depth and complexity. (Raimi et al., 2021; Vass, 2025). With over 3 million inactive or orphaned wells in the U.S. alone, many pose environmental risks, including methane emissions with a global-warming potential 25–34 times that of CO₂ (IPCC, 2021; Kang et al., 2019). In Pennsylvania, legacy wells contribute an estimated 5–8% of anthropogenic methane, prompting stricter plugging regulations (Kang et al., 2019; Osundare et al., 2018). Equally important is reclaiming and stabilizing well sites for land safety and ecosystem restoration. Recent advances in high-fidelity reservoir simulation paired with supervised machine-learning surrogates—such as artificial neural networks, Random Forests, and Gaussian-process emulators—now enable rapid, data-driven screening of thousands of well candidates to pinpoint those most suitable for conversion into geothermal energy assets, effectively transforming environmental liabilities into low-carbon, revenue-generating infrastructure (Jin et al., 2022; Duplyakin et al., 2022; Rohmer et al., 2023).

Amid these challenges, data-driven revitalization of depleted hydrocarbon fields has emerged as a powerful strategy to tackle renewable-energy intermittency and advance global decarbonization goals (Duggal et al., 2022; Gayayev, 2023). Recent policy roadmaps from the International Energy Agency, and the World Economic Forum (IEA, 2017; REN21, 2019; WEF, 2021) all highlight the need for adaptable, seasonal storage solutions that bolster energy security and grid flexibility (Dincer & Rosen, 2011; Van Der Roest et al., 2021). Repurposing depleted reservoirs directly supports this agenda: their well-characterized permeability, porosity, and extensive historical datasets can be leveraged to transform environmental liabilities into low-carbon, revenue-generating assets (Green et al., 2021; Lee, 2013; Li, 2016). Diverse pore-space applications—CO₂ sequestration (Qin et al., 2023), hydrogen storage (Zhu et al., 2024), subsurface electricity generation (Duggal et al., 2022), and Aquifer Thermal Energy Storage (ATES) (Matos et al., 2019; Stricker et al., 2020)—underscore their versatility in a decarbonizing energy landscape. Crucially, advances in high-resolution reservoir modelling and supervised machine-learning algorithms now allow rapid screening of thousands of well trajectories, prediction of storage performance under uncertainty, and optimisation of injection–production schemes—dramatically reducing time-to-deployment and de-risking investment (Khosravi et al., 2024; Liu et al., 2024).

With over 2,800 systems worldwide—mainly in the Netherlands supplying ~2.5 TWh/yr—ATES has proven technically mature and scalable (Dickinson et al., 2009; Fleuchaus et al., 2020; Kastner et al., 2017; van Heekeren & Bakema, 2015). It uses existing wells and data to store surplus heat for later use, helping balance energy supply and support renewable integration (Fleuchaus et al., 2020; Paksoy et al., 2000). High-Temperature ATES (HT-ATES), targeting fluid temperatures ≥ 90 °C, further extends storage potential. Recent studies now combine coupled reservoir simulations with supervised-learning surrogates—such as neural-network emulators (Jin et al., 2022), Gaussian-process metamodelling (Rohmer et al., 2023), and Random Forest regressors (Duplyakin et al., 2022)—to forecast system performance across diverse geological scenarios in seconds rather than hours. In the Upper Rhine Graben, 90% of surveyed depleted oil fields are HT-ATES-ready, with projected storage of up to 12 GWh/yr and ~82% recovery after ten years (Holstenkamp et al., 2017; Liu et al., 2024; Stricker et al., 2020). Field-scale pilots confirm feasibility: the Middenmeer project stored 85–90 °C heat in a 400 m aquifer, using real-time monitoring and MODFLOW/MT3DMS to manage issues like sand production (HEATSTORE, 2025; Oerlemans et al., 2022). In the U.S., geothermal battery projects in California and Texas repurpose oilfields for long-duration thermal storage, leveraging machine-learning surrogates to optimise injection and reduce testing time by over 70% (Khosravi et al., 2024; Liu et al., 2024; Zhu et al., 2024). These systems offer superior discharge duration over batteries and cut CO₂ emissions (Fleuchaus et al., 2020).

Depleted reservoirs offer vast storage potential. The Carrizo-Wilcox aquifer could store 554 TWh of heat—63 TWh as electricity (Akindipe et al., 2024)—while the Upper Rhine Graben could supply 10 TWh/yr of heat, covering much of the region's demand (Stricker et al., 2020). These capacities far exceed those of pumped hydro and grid batteries (World Energy Outlook 2021 – Analysis, 2021). Repurposing wells cuts capital costs, and integration with renewables enables low levelized costs—USD 0.11/kWh for electricity, USD 0.02/kWh for heat (Anttila, 2021; Zhu et al., 2024). Though round-trip efficiency is 40–50%, geological storage fills a key seasonal balancing role (Van Der Roest et al., 2021). Although high-temperature ATES is advancing, substantial uncertainties still cloud the design and operation of dual-well systems in heterogeneous clastic reservoirs—uncertainties that brute-force deterministic simulations cannot efficiently close. Chief among these is selecting the hot-cold well spacing that suppresses premature thermal breakthrough yet preserves a high heat-recovery factor. Thermal breakthrough can erode efficiency and jeopardise long-term viability (Bloemendal & Hartog, 2018; Sommer et al., 2013). While ATES research is growing, few studies have formally linked spacing to thermal performance in repurposed hydrocarbon fields with complex stratigraphy (Kastner et al., 2017; Pellegrini et al., 2019). Cutting-edge work now couples high-fidelity thermo-hydraulic models with machine-learning surrogates—Random-Forest, artificial-neural-network, and Gaussian-process emulators—to scan thousands of spacing scenarios, propagate geological uncertainty, and pinpoint Pareto-optimal dual-well layouts in minutes rather than days (Jin et al., 2022; Rohmer et al., 2023; Duplyakin et al., 2022). Yet these data-driven methods have not been systematically deployed in depleted clastic basins (Liu et al., 2024; Khosravi et al., 2023).

This study addresses four tightly linked questions aimed at advancing seasonal thermal-energy storage in depleted clastic reservoirs. It first pinpoints the inter-well distance in a dual-well configuration that suppresses premature thermal breakthrough while maximizing heat-recovery efficiency. It then examines how key reservoir attributes in Hungarian clastic formation—porosity, permeability, and lithological anisotropy —shift the optimal spacing and shape of the overall system performance. Next, it evaluates which alternating operating schedule—summer hot-storage/winter hot-production versus summer cold-production/winter cold-storage—delivers the greatest annual energy return. Finally, it tests whether supervised machine-learning models, calibrated on suites of coupled MODFLOW-MT3DMS simulations, can predict and optimize these design variables across the full inventory of candidate wells, thereby slashing computational time and accelerating field deployment. A preliminary version of this study was presented as an abstract at the 16th European Geothermal PhD Days, held in Szeged, Hungary (H. Abdulhaq, 2025b).

Hungary is an exceptional testbed for high-temperature ATES because it (i) hosts a dense network of depleted oil-and-gas wells situated close to district-heating loads, (ii) lies within a moderate-to-high geothermal-gradient province that delivers initial reservoir temperatures of ≈ 70 °C, and (iii) operates under strong national directives for decarbonisation and energy-security gains (Nádor et al., 2022; J. D. Szanyi et al., 2025). Decades of well-log, core, and production records have generated a richly labelled subsurface dataset that can be mined with geostatistics and supervised learning (Topór et al., 2023). Leveraging these assets, we integrate ML with high-resolution MODFLOW/MT3DMS heat-transport simulations and Random-Forest surrogates—an approach shown to cut optimisation runtimes by an order of magnitude while preserving predictive accuracy (Jin et al., 2022; Duplyakin et al., 2022). The resulting workflow delivers both mechanistic insight and actionable design rules for dual-well HT-ATES implementation in Hungarian clastic reservoirs.

By coupling high-fidelity MODFLOW-MT3DMS heat-transport simulations with Random-Forest surrogates trained on thousands of synthetic well trajectories, we cut optimisation runtimes by more than 90 % while maintaining high predictive skill (Jin et al., 2022; Rohmer et al., 2023). Feature-importance analysis of these meta-models (R² ≈ 0.87 on a hold-out set) pinpoints inter-well distance, reservoir anisotropy, and cycle length as the dominant controls on Heat-Recovery Factor—insights that translate into practical spacing rules and seasonally phased operating schedules for Hungarian clastic reservoirs (Duplyakin et al., 2022; Jin et al., 2022). Validation against independent datasets from analogous Central European basins confirms the transferability of the workflow (Topór et al., 2023). Collectively, these results provide a scalable blueprint for converting depleted hydrocarbon assets into long-duration, low-carbon thermal batteries, simultaneously reducing decommissioning liabilities and advancing national energy-security and climate-mitigation goals (REN21, 2019; WEF, 2021).

2. Glossary of Terms

ATES (Aquifer Thermal Energy Storage): A technology for storing and retrieving thermal energy in aquifers, enabling seasonal energy management by injecting heat in summer and recovering it during winter.

HT-ATES (High-Temperature Aquifer Thermal Energy Storage): An advanced form of ATES designed for storage and recovery of thermal energy at fluid temperatures ≥90 °C, suitable for industrial and district heating applications.

MODFLOW: A modular three-dimensional finite-difference groundwater flow model developed by the U.S. Geological Survey, widely used for simulating groundwater conditions and flows.

MT3DMS: Modular Three-Dimensional Multi-Species Transport Model, used in conjunction with MODFLOW to simulate the transport of heat, solutes, or other contaminants in groundwater systems.

Random Forest: An ensemble machine learning method based on decision trees, used for regression and classification tasks, valued for its robustness and ability to model complex relationships.

Heat-Recovery Factor (HRF): The ratio of recovered thermal energy to the initially injected energy during an ATES cycle, often used as a performance metric for system efficiency.

Thermal Breakthrough: The phenomenon where injected hot or cold fluid reaches the production well too quickly, reducing system efficiency and potentially shortening operational lifetime.

UCN File (Unformatted Concentration File): A binary output file generated by MODFLOW/MT3DMS containing spatially and temporally resolved simulation results, in this study representing temperature distributions.

Residual Heat Accumulation: The progressive build-up of stored heat in the aquifer over multiple ATES cycles, typically leading to higher thermal recovery efficiencies over time.

Surrogate Model: A fast-running, data-driven model (e.g., Random Forest, neural network) trained on outputs from complex numerical simulations to predict system behavior efficiently.

Hydrogeological Model: A numerical model simulating groundwater flow based on hydraulic and geological parameters to understand subsurface water movement and storage characteristics.

Thermal Recovery Efficiency: The percentage of injected thermal energy that can be successfully recovered during production phases in a seasonal thermal energy storage system.

Overpressure: Subsurface pressure exceeding hydrostatic pressure, often due to geological compaction, tectonic forces, or fluid generation processes, which can influence reservoir behavior.

Pannonian s.l. (sensu lato):

A stratigraphic term referring broadly to the Upper Miocene sedimentary sequences in the Pannonian Basin, including formations such as Újfalu and Zagyva.

Hot Well: A well designated for the injection of heated water during storage periods and extraction during production periods in an ATES system.

Cold Well: A complementary well used to manage temperature balance in an ATES system, typically used for extracting cooler water during storage or injecting cooler water during production, depending on the operating scheme.

Stress Period: A defined time interval in a MODFLOW/MT3DMS simulation during which external stresses (e.g., injection, pumping) are assumed constant or follow a specified pattern.

Permutation Importance: A machine learning method used to assess the relative importance of input features by measuring the decrease in model performance when feature values are randomly shuffled.

3. Background and Regional

3.1. Geological and Hydrogeological Setting of the Békés Basin

The Pannonian Basin is a sedimentary basin located in East-Central Europe, characterized by a complex geological structure consisting of variously subsided basins and horst-like blocks. The basement primarily comprises metamorphic Paleozoic rocks, with Mesozoic carbonate formations present in some areas that can serve as good aquifers (Horváth et al., 2015). Within this larger geological context, the Békés Basin represents one of the two main depressions of the Southern Great Plain of Hungary, alongside the Makó Depression, with these two significant depressions divided by the Battonya Ridge (Juhász, 1991). The Békés Basin is particularly notable for its exceptional depth, reaching approximately 7,000 meters of post-Cretaceous sedimentary fill (USGS, 2023), making it one of the deepest sub-basins within the Pannonian Basin system. Figure 1 shows the location of the study area within Hungary, along with the lithological map highlighting the modelled section of the Békés Formation and a corresponding lithological cross-section.

The stratigraphic sequence of the Békés Basin follows the general pattern of the Pannonian Basin, with important variations in thickness and characteristics. At the beginning of the Lower Pannonian period, the Endrőd Marl Formation was deposited, consisting of calcareous marl and clay marl. This formation is overlain by the fine sand turbidite set of the Szolnok Formation, which reaches several hundred meters in thickness in some locations. Above the turbidites, particularly in shallower basin areas, the hemipelagic marls are covered by the thick clayey-silty layers of the Algyő Formation with a prodelta facies (Haas, 2013). A key characteristic of this formation sequence is the extremely high overpressure below and throughout the set. The sand content of the Algyő Formation increases in areas with a shallower basement, allowing the upper part of the formation to function as a water-bearing unit in certain locations. Generally, however, the Lower Pannonian formations exhibit poor water-bearing characteristics.

The Pannonian s.l. sequence, which overlies the Lower Pannonian layers, consists primarily of the Újfalu Formation and the Zagyva Formation. The Újfalu Formation, characterized by delta front and delta plain facies, represents the most hydrogeologically significant Pannonian s.l. sediment. The Zagyva Formation features deltaic background and alluvial plain facies, with dominant sediments being bed-filling and bay-mouth bar deposits that demonstrate good water-bearing properties despite their limited horizontal dimensions. These formations are hydrodynamically connected through multiple linear erosions and overlapping (Juhász, 1991). In the Békés Basin region, the bottom of the Pannonian s.l. sequences typically lie at depths of 2,000-2,500 meters from the ground surface, with the total thickness of Pannonian s.l. sediments exceeding 2,000 meters in the Békés Basin—among the thickest in the entire Pannonian region.

3.2 Hydrodynamic Systems and Pressure Regimes

The Békés Basin, like the broader Carpathian basin, features two distinct flow regimes: an upper, gravity-driven flow system within the Pannonian s.l. sequences, and a deeper, overpressure-driven system within the Lower Pannonian formations, primarily affecting the finer deep-sea sediments and underlying formations (Mádl-Szőnyi & Tóth, 2009; Tóth & Almási, 2001). The overpressure in the deeper system is remarkably high, reaching up to 40 MPa above hydrostatic pressure. This extreme overpressure primarily results from tectonic compression of the formations, with additional contribution from gas formation during sediment maturation processes (Tóth and Almási, 2001).

In the Békés Basin region, pressure-depth profiles indicate that the dynamic pressure gradient exceeds the hydrostatic pressure by approximately 0.13 MPa (equivalent to about 13 m hydrostatic head) in Quaternary formations and by approximately 0.44 MPa (about 44 m hydrostatic head) in the Pannonian s.l. sequence. The Lower Pannonian sequence exhibits even more dramatic super-hydrostatic pressure, with the dynamic pressure gradient exceeding hydrostatic pressure by more than 60 MPa. This significant pressure differential creates complex hydrodynamic conditions that must be carefully considered when designing and implementing any subsurface fluid management system, including seasonal thermal energy storage.

3.3. Reservoir Properties and Geothermal Potential

The Pannonian s.l. sandstone reservoirs in the Békés Basin region exhibit favorable characteristics for geothermal applications, with effective porosity values typically reaching 22-25%. The permeability of these Pannonian s.l. reservoirs, which consist of highly permeable sand layers, can reach up to 2000 mD (1.97 × 10^-12 m²), corresponding to a hydraulic conductivity of 5-10 m/day(Bálint & Szanyi, 2015; Korim, 1991; J. Szanyi et al., 2015). These values represent some of the most favorable reservoir conditions in the Hungarian geothermal context.

The consolidation state of the sandstone varies depending on depth and cementation processes. The sandstone can range from consolidated to unconsolidated, with cementation typically occurring through quartz overgrowth, calcite, or kaolin precipitation. The degree of cementation significantly influences both porosity and stability, particularly during production and injection operations. Generally, sandstone induration increases with depth as cementitious material precipitates into the pores from fluid extracted during compaction. The sand bodies are typically separated by thinner fine-grained sediments, creating a complex, heterogeneous reservoir structure (Korim, 1991; Bálint and Szanyi, 2015).

The Békés Basin is characterized by an exceptionally high geothermal gradient, approximately 50°C/km, significantly above global averages due to the relatively thin crust beneath the Pannonian Basin (Lenkey et al., 2021). This elevated thermal gradient results in reservoir temperatures around 70°C at depths of approximately 1,500-1,800 meters, making these formations particularly suitable for thermal energy storage and district heating applications. At greater depths of 2,000-2,500 meters, temperatures can reach 90-120°C, offering potential for higher-temperature applications.

The combination of favorable reservoir properties (high porosity and permeability) and outstanding geothermal conditions makes the Békés Basin exceptionally suitable for geothermal energy utilization, including seasonal heat storage applications (J. Szanyi et al., 2009). The region's depleted hydrocarbon wells, many of which penetrate these favorable Pannonian s.l. sandstone formations, present valuable opportunities for repurposing as components of thermal energy storage systems.

3.4. Hydrocarbon History and Well Infrastructure

The Békés Basin has a rich history of hydrocarbon exploration and production, with extensive drilling activities dating back to the mid-20th century. These activities have resulted in a substantial inventory of wells throughout the region, many of which have now reached the end of their productive lifespan as hydrocarbon producers. The basin contains significant natural gas resources, with gases produced from multiple reservoir intervals at depths ranging from 1,800 to 2,900 meters (Survey, 2023).

These Neogene sedimentary sequences overlying the basement highs have demonstrated the best hydrocarbon reservoir characteristics in the southeastern part of Hungary (Horváth et al., 2015). The extensive exploration and production history has generated valuable geological and reservoir data, including detailed information on formation properties, temperature profiles, pressure regimes, and fluid characteristics. This wealth of data provides a significant advantage for assessing the potential of these formations for thermal energy storage applications.

The existing well infrastructure, though aging, offers potential for repurposing rather than decommissioning. Many wells have been completed with telescopic designs, with casing diameters ranging from approximately 340 mm (13 3/8") at the surface to 140-178 mm (5 1/2" - 7") at reservoir depths. While some wells may require workover or partial recompletion to ensure mechanical integrity for long-term thermal storage operations, the basic infrastructure represents a valuable asset that could significantly reduce the capital costs associated with implementing seasonal heat storage systems.

3.5. Relevance to Seasonal Heat Storage

The hydrogeological characteristics of the Békés Basin make it particularly suitable for Aquifer Thermal Energy Storage (ATES) applications, especially High-Temperature ATES (HT-ATES) targeting temperatures up to 90°C or higher. The Pannonian s.l. sandstone formations, with their high porosity, good permeability, and favorable temperature conditions, provide an excellent medium for seasonal storage and retrieval of thermal energy.

The proximity of many depleted wells to population centers in the region creates opportunities for integrating seasonal heat storage with district heating systems, similar to successful implementations in other parts of Hungary, such as Szeged. The initial reservoir temperatures of approximately 70°C in the target formations are ideal for enhancement through additional heat input during summer months, with subsequent extraction during winter heating periods.

However, the complex pressure regimes, particularly the significant overpressure in deeper formations, present challenges that must be carefully managed. Additionally, the heterogeneous nature of the reservoir formations, with sand bodies separated by fine-grained sediments, creates potential for compartmentalization that could affect thermal breakthrough patterns between injection and production wells.

The dual-well system proposed for seasonal heat storage—with one well serving for summer hot storing/winter hot producing and another for summer cold producing/winter cold storing—must be carefully designed to account for these hydrogeological characteristics. The optimal spacing between wells must balance the need to prevent premature thermal breakthrough while maximizing energy recovery efficiency, taking into consideration the specific reservoir properties of the Békés Basin formations.

By leveraging the extensive geological knowledge, existing well infrastructure, and favorable reservoir conditions of the Békés Basin, seasonal heat storage systems can be optimized to provide sustainable, efficient thermal energy solutions while extending the productive life of otherwise abandoned hydrocarbon assets.

4. Materials and Methods

4.1. Methodological Framework

The methodological framework of this study began with the comprehensive curation and preparation of existing data from abandoned hydrocarbon fields, which involved extensive data modeling, cleaning, and refinement. A hydrogeological model was then constructed using MODFLOW (Harbaugh, 2005), providing the foundational flow and transport parameters needed for subsequent analyses. Building on this, a heat transport model was employed to simulate thermal performance in potential Underground Thermal Energy Storage (UTES) candidates. The resulting simulation outputs served as the training dataset for a Random Forest algorithm (Breiman, 2001; Pedregosa et al., 2011), which was designed to predict thermal performance in other areas. Finally, the predictions generated by the Random Forest were validated against actual simulation results. Figure 2 outlines the primary steps involved in this integrated approach.

4.1. Data Collection and Data Preparation

For the data collection segment of our study, we utilized well data from the southern part of the Bekes Basin in Hungary—a site formerly exploited for hydrocarbons and now recognized for its geothermal potential (Kovács & Teleki, 1994). This field encompasses two key formations: the shallow Szolnok Formation, which functions as an aquifer (Varga et al., 2019), and the deeper Bekes Formation, thereby providing a unique opportunity to evaluate heat storage capabilities. We integrated a diverse set of data—including core samples, density logs, resistivity logs, and gamma-ray measurements—using a stochastic simulation process. SGeMS was employed for the geostatistical simulation, while RockWorks facilitated the integration of simulation outcomes for gross thickness, effective porosity, and permeability. A total of 100 stochastic realizations were generated for each grid or voxel point, with the median (Md-type estimation) calculated to represent the central tendency, effectively minimizing the impact of skewed or outlier values. For the Bekes Formation, this Md-type estimation was deemed most representative of the expected geological parameters.

In a prior study by Abdulhaq et al., (2024), the southeastern section of the study area was identified as a prime candidate for energy storage due to the Bekes Formation’s average temperature of approximately 70°C, coupled with its favorable porosity and permeability characteristics. Based on these insights, a polygon delineating this promising area was selected for detailed hydrogeological and heat transport simulations. Figure 3 illustrates the distances from each well to the candidate town, with only those wells penetrating the Bekes Formation included in the analysis—wells that did not extend into the Bekes were excluded. The figure also displays the temperature distribution within the Bekes Formation. By applying a 5 km threshold for effective thermal transport and selecting wells with temperatures below 70°C as potential candidates for HT-ATES, only wells marked with red borders were retained for further simulation of their thermal performance (Fig 4).

4.3. Hydrogeological Modelling

The groundwater flow within the Bekesi Conglomerate reservoir was modeled using MODFLOW-2000, a modular finite-difference groundwater flow modeling software developed by the U.S. Geological Survey (Version 1.19.01, March 25, 2010). The groundwater flow model was developed using MODFLOW and discretized into three layers, 197 rows, and 400 columns with a cell size of 10 by 10 meters. In this configuration, the second layer represents the Bekes Formation. The model parameters were defined at either the cell-by-cell or grid scale, based on data-driven simulations and assumptions derived from field data and laboratory analyses. Table 1 summarizes the key parameters used in the MODFLOW processing model.

4.4. Heat Transport Modelling

For heat transport modeling, we employed MT3DMS within the GMS framework. Following the successful initiation and simulation of the MODFLOW model, we activated the Basic Transport Package in MT3DMS by introducing the starting temperature for each cell. To simulate the heat transport processes, we selected the advection, dispersion, source/sink mixing, and chemical reaction modules. The parameters utilized for this heat transport model are listed below.

Parameter	Value / Description	Justification
Initial Temperature	Varies spatially	Derived from drill stem tests and bottom-hole temperature data
Advection Package	Third order TVD scheme Ultimate	Selected for numerical stability and accuracy
TRPT	0.1	Assumed based on typical sedimentary conditions (Gelhar et al., 1992)
TRVT	0.01	Assumed based on typical sedimentary conditions (USGS, 2022)
DMCOEF (Effective Molecular Diffusion Coefficient)	0.01 m²/day	Literature-based estimate (ModelMuse, 2024)
longitudinal Dispersivity	Varies with lithology	Based on the Rock Type Calculation and thermal conductivity of Bekes Fm. from (Vass et al., 2018)
Sorption	Linear isotherm	Common assumption for initial reactive transport modelling
Kinetic Rate Reaction	zero order reaction	Assumed for simplification of reactive processes
Preconditioner	Jacobi	Default iterative solver preconditioner

4.5. Simulation Setting

To replicate ATES operations, the stress periods were structured into one month of system downtime, followed by five months dedicated to heat storage, another one-month break, and five months allocated for heat production, repeating this cycle over a seven-year period. A pair of wells was selected where the distance between them is more than 500 meters: one well for hot injection in summer and subsequent hot production in winter, and the other well for cold production in summer with cold injection during winter.

4.6. Training Data for Machine Learning

To ensure the integrity and consistency of our training dataset, we developed an in-house data entry module using Streamlit that integrates both manual input and CSV-based uploads. This module facilitates the incorporation of critical well parameters—including well ID, name, spatial coordinates, porosity, permeability, gamma ray measurements, thickness, distance to the cold well, and initial temperature—while also allowing for the efficient addition of temporal temperature data. For the temperature data, the module supports file uploads in CSV format, validates the presence of required columns (TimeDays, WellType, Row, Col, Temperature), and links the data to the corresponding well ID via an interactive selection box (Fig. 5). The system provides real-time previews and feedback, ensuring that all data are consistently and accurately stored in the database.

The training process involves leveraging a SQLite database that consolidates data from two tables—one containing well properties and the other recording temperature measurements. First, the data from both sources are merged to form a comprehensive DataFrame, which includes attributes such as porosity, permeability, gamma ray, thickness, distance to cold wells, initial temperatures, and the time variable, while spatial reference columns (Row, Col) are retained only for reference. After cleaning the dataset by removing any rows with missing values, the data is split into training and testing subsets. A Random Forest regressor is then trained on the assembled features (excluding the spatial reference columns) to predict well temperatures, with performance evaluated through metrics including MAE, RMSE, and R². Finally, the trained model is serialized and saved to a predefined path, ensuring reproducibility and ease of deployment within our predictive framework.

The prediction phase leverages a pre-trained Random Forest model to forecast temperature evolution for all available wells over a specified time range. In this stage, users interact with a Streamlit-based interface where they define the starting time, ending time, and time step for the predictions. The system loads the trained model and, for each well in the database, constructs an input DataFrame populated with static parameters such as porosity, permeability, gamma ray, thickness, distance to the cold well, and initial temperature, while dynamically varying the time parameter. The model then predicts the temperature for each time step, and the results are combined into a unified dataset. To facilitate analysis, the predicted curves for individual wells are plotted on a single graph—each curve clearly labeled with its corresponding WellID—providing a comprehensive visualization of thermal performance over time. The complete repository of the scripts is available online; however, the data content is not included (H. Abdulhaq, 2025a).

4.7. Model Calibration and Validation

For model calibration and validation, the hydrogeological simulation was first calibrated against historical field data—specifically, reported hydraulic head measurements—to ensure that the MODFLOW model accurately reflected the aquifer’s behavior (Anderson et al., 2015). These historical observations provided key benchmarks that guided the tuning of parameters related to flow and transport processes. Concurrently, the machine learning component was iteratively calibrated using simulation outputs as training data. As the simulation yielded more refined thermal performance data, continuous retraining of the Random Forest model progressively improved its predictive accuracy (Khosravi et al., 2024). This integrated calibration strategy, leveraging both empirical historical data and dynamically generated simulation inputs, has enhanced the overall robustness and reliability of the modeling framework.

4.8. Sensitivity Analysis

The sensitivity analysis of the machine learning model revealed that TimeDays is by far the most influential parameter affecting temperature prediction outcomes, with a permutation importance score of 1.96, substantially higher than all other features. This result underscores the dominant role of temporal evolution in determining thermal behavior within the reservoir, reflecting the accumulation and dissipation of heat over successive cycles. Other features—such as Initial Temperature (0.034), Thickness (0.029), GammaRay (0.006), and Permeability (0.006)—exhibited comparatively minor effects, suggesting that while geological and petrophysical properties contribute to system performance, their impact is secondary to the time dimension (Fisher et al., 2019). Notably, Porosity had no measurable influence, and DistanceToColdWell showed a negligible effect (0.00007), indicating that these factors may be less critical under the modeled conditions or were inadequately represented in the available dataset. These findings provide insight into which parameters should be prioritized for accurate performance prediction and model refinement.

5. Results

5.1. Heat Simulation result

The output of the heat simulation is stored in a UCN file, which encapsulates the final temperature distribution computed by the model (Ishikawa et al., 2014). To efficiently leverage this output, we developed an in-house Python module designed to load, process, and analyze UCN files (Fig 6). This module automatically identifies hot and cold wells and allows users to designate observation wells within the simulation setting. Users can load the UCN file and select specific layers of interest, while the integrated visualization sub-module generates plots of temperature versus time for defined stress periods, providing clear insights into the thermal performance of individual wells. Additionally, an animation sub-module enables dynamic playback of temperature evolution by allowing adjustable speeds and frame skips, thereby enhancing interpretability. Furthermore, a dedicated recovery efficiency sub-module computes the thermal recovery efficiency for any simulated wells. The complete suite of module scripts is available online and can be accessed independently of the dataset (H. Abdulhaq, 2025c).

In our thermal simulation studies, two sets of candidate results were obtained, both exhibiting robust injection performance with maximum temperatures consistently around 94.9°C across cycles (Fig. 7). In the first candidate, the break phases—representing reservoir conditions when the system strikes—showed a gradual increase in maximum temperatures from 79.67°C up to 84.33°C, with thermal recovery efficiency improving from 83.92% to 88.82%. In the second candidate, although the initial break phase efficiency was lower at 81.05%, a notable enhancement was observed, with the efficiency rising to a maximum of 87.93% over repeated cycles. Key performance metrics for this candidate include a ratio of last-to-first efficiency of 1.08, a percent increase of 8.50% from the initial break phase, an average efficiency of 85.46%, and a slope of 1.04 per cycle, indicating a steady improvement with each cycle. Importantly, while these reservoir-level efficiency improvements are significant, the production efficiency during winter remains even higher than these baseline values, ensuring superior operational performance. Together, these observations demonstrate that repeated cycles of heat injection and reservoir cooling not only stabilize the thermal regime but also enhance both the inherent recovery and the actual production efficiency during winter operations.

5.2. Machine Learning Result:

When initiating temperature performance predictions for candidate wells, we generated thermal efficiency results for multiple wells, predominantly located outside the originally modeled area. This indicates that our predictive approach is effective even beyond the immediate boundaries of our hydrogeological model, highlighting the potential for broader applicability (Fig. 8). Analysis of these external candidates reveals initial efficiencies ranging between approximately 80% to 84%, with a consistent improvement observed through repeated injection and recovery cycles. Specifically, efficiencies showed an average of 85.46%, and maximum values reached up to 88.84%. The ratio of final-to-initial efficiency averaged around 1.05, corresponding to an overall efficiency improvement of up to 8.01%. Furthermore, the observed positive slopes (up to 1.0 per cycle) clearly illustrate that efficiency systematically increases over successive cycles. These findings underscore the reliability and robustness of our predictive methodology in forecasting thermal performance beyond the initially calibrated region, as visualized in Figure 9, which displays a distribution map of each well and the corresponding predicted improvement in thermal performance.

6. Discussion

6.1. Alignment with Previous Studies and Theoretical Outcomes

Our simulation and machine learning outcomes align closely with prior theoretical and empirical findings related to High-Temperature Aquifer Thermal Energy Storage (HT-ATES). Previous studies indicate a typical improvement in thermal recovery efficiency with each successive injection-production cycle, attributed primarily to residual heat accumulation within the aquifer (Collignon et al., 2020; Drijver et al., 2012). Our simulation results reflect similar trends, with initial efficiencies around 80-84% gradually increasing to as high as approximately 88%. This progressive efficiency improvement closely mirrors published benchmarks from international case studies, which typically report HT-ATES efficiencies stabilizing in the range of 60–80% after several operational cycles (Collignon et al., 2020; Winterleitner et al., 2018).

The stability of our modeled temperature curves over multiple years also echoes the theoretical predictions, which suggest that a thermal equilibrium or steady state emerges after multiple cycles (Tang & Rijnaarts, 2023). However, subtle contrasts exist; for instance, our simulations may reflect idealized conditions, potentially omitting complexities such as density-driven convection or significant vertical heat migration reported in heterogeneous clastic reservoirs (Winterleitner et al., 2018). This highlights a limitation where idealized models may yield slightly optimistic efficiency and thermal stability predictions compared to more complex real-world scenarios.

6.2. Key Influencing Parameters

Our machine learning analysis identified time ('TimeDays') as the most influential parameter on temperature and efficiency predictions, surpassing all other geological and operational parameters by a substantial margin. Specifically, the permutation importance of time was approximately 1.96, whereas other parameters like initial temperature, thickness, gamma ray, and permeability demonstrated notably lower impacts (0.034, 0.029, 0.006, and 0.006, respectively). This finding emphasizes that temporal factors—specifically, cumulative cycle duration and residual heat buildup—predominantly govern thermal performance, consistent with the literature where residual thermal energy strongly influences long-term operational outcomes (Drijver et al., 2012).

6.3. Strengths and Limitations of Hydrogeological Model Calibration

Calibrating the hydrogeological model using historical head data provides several advantages, chiefly realistic and site-specific insights into subsurface dynamics. However, this method's main limitation lies in the necessity to modify boundary conditions to match historical measurements adequately, potentially introducing bias or oversimplification. Adjustments made to boundary conditions, while essential for aligning simulations with empirical data, might restrict the generalizability of the hydrogeological model to conditions significantly different from historical scenarios.

6.4. Enhancing Decision-Making for UTES Site Selection

Our predictive methodology substantially accelerates the decision-making process for selecting potential UTES sites by providing rapid, high-quality predictive outcomes derived from historical data and simulations. By transforming abandoned hydrocarbon reservoirs into a data-driven analytical framework, our approach reduces evaluation times and increases confidence in site assessments. This can significantly streamline site selection workflows, particularly valuable in regions with numerous abandoned wells and substantial historical datasets.

6.5. Implications for Scaling Geothermal Storage Projects

The successful application of this predictive approach using gamma ray logs to differentiate rock types indicates substantial potential for scaling geothermal storage projects, particularly within clastic sedimentary basins. Internationally, numerous clastic basins exist with similar sedimentological characteristics, making this methodology widely applicable. This capability facilitates rapid assessments across varied geographies, promoting efficient, scalable deployment of HT-ATES technologies globally.

6.6. Assumptions and Simplifications

Key assumptions and simplifications within our modeling framework include idealized boundary conditions, homogeneous or simplified geological heterogeneities, and consistent thermal properties across the model domain. Such assumptions likely influence the accuracy and predictive power of the model, potentially limiting its applicability to real-world scenarios with pronounced geological complexity or significantly variable hydrogeological conditions.

6.7. Uncertainties in Input Data

A major uncertainty arises from the inherent accuracy and reliability of the initial simulation data used for machine learning training. Because our predictive model relies heavily on accurate thermal simulations, any errors or oversimplifications in the initial simulations propagate through the predictions, potentially affecting reliability. This highlights the critical importance of accurate and comprehensive simulation input data.

6.8. Recommendations and Future Work

Future studies should focus on enhancing simulation fidelity, incorporating more detailed heterogeneity, and conducting diverse scenario sampling to bolster the robustness of the machine learning model. Additionally, exploring advanced physics-informed machine learning techniques could significantly improve predictive accuracy and generalizability. Extending simulation runs, including more diverse geological and operational parameters, and conducting further validation with independent field data would further strengthen confidence in our model predictions and their applicability to broader contexts.

7. Conclusion

This study set out to advance the design and deployment of high-temperature Aquifer Thermal Energy Storage (HT-ATES) in depleted clastic reservoirs through four interlinked research questions. First, by coupling high-resolution MODFLOW-MT3DMS heat-transport simulations with an in-house Python toolkit, we identified the optimal inter-well spacing that minimizes premature thermal breakthrough while maximizing cumulative heat-recovery efficiency. Our two candidate well pairs consistently achieved peak injection temperatures near 94.9 °C and demonstrated steady efficiency gains—up to 8.5 % over repeated cycles—when spaced to balance thermal front propagation and lateral heat recharge.

Second, we quantified how key reservoir properties in Hungarian clastic formations modulate this optimal spacing and overall system performance. Sensitivity analysis revealed that porosity and permeability variations shift thermal breakthrough timing and adjust peak recovery efficiencies by several percentage points, while anisotropy primarily affects the shape of the thermal plume. These insights enable tailoring well spacing to site-specific hydraulic and thermal heterogeneities.

Third, comparison of seasonal operating schedules confirmed that the conventional hot-storage/winter production cycle yields marginally higher annual energy returns than the inverse (cold-production/winter storage), owing to more effective residual-heat carryover. Specifically, winter production efficiencies exceeded baseline recovery values, underscoring the value of aligning storage and demand cycles to ambient temperature differentials.

Finally, our Random Forest surrogate models—trained on several coupled simulations—proved capable of predicting thermal recovery efficiency across a broad inventory of candidate wells outside the original calibration domain. Surrogate predictions achieved average accuracies within 2–3 % of full numerical simulations, with efficiency improvements up to 8 % over multiple cycles. This machine-learning workflow accelerates design optimization by an order of magnitude, enabling rapid screening of sites and operational schedules.

Collectively, these findings deliver a data-driven framework for HT-ATES implementation in depleted clastic reservoirs: mechanistic insights into spacing and scheduling trade-offs, parameter-specific performance adjustments, and a scalable surrogate modeling approach for design optimization. Future work should extend this framework by incorporating more complex heterogeneities, exploring physics-informed learning techniques, and validating predictions against field pilot data to further refine site-selection criteria and operational guidelines.