A Data-Driven Sand Channel Prediction ML–Model from Well Log Data for Optimized UTES Site Selection and Thermal Breakthrough Prevention: Case Study from Hungary

1. Introduction

The accelerating global shift toward renewable energy has significantly heightened the demand for effective and scalable energy storage solutions, which are essential for addressing renewable intermittency, ensuring grid flexibility, and providing reliable seasonal energy storage [1,2,3] . While short-term battery storage technologies have matured significantly, their inherent limitations in scalability and cost-effectiveness for long-duration or seasonal storage underscore the need for alternative, large-scale solutions. Geological energy storage, particularly Underground Thermal Energy Storage (UTES), is increasingly recognized as a promising alternative due to its capacity for extensive seasonal storage and its minimal environmental footprint [4,5,6] .

UTES systems, notably Aquifer Thermal Energy Storage (ATES), are now well-established worldwide, with more than 2,800 operational installations—most prominently in the Netherlands, where they collectively provide approximately 2.5 TWh of thermal storage annually [5,7,8]. Recent technological advances, including High-Temperature ATES (HT-ATES) systems that operate at temperatures above 90 °C, further extend the potential of UTES by enabling efficient integration with renewable energy sources, industrial heat supply, and district heating networks. These developments substantially expand both the practical impact and the relevance of UTES technologies [9,10].

A critical determinant of UTES system effectiveness is the strategic placement of injection and production wells, which directly influences the efficiency, reliability, and lifespan of thermal storage operations. One of the most significant operational challenges is thermal breakthrough, where injected thermal fluids prematurely reach production wells, dramatically reducing heat-storage efficiency and overall economic viability [4,11]. Preventing or mitigating thermal breakthrough fundamentally depends on precise geological characterization—specifically, understanding the spatial distribution and geometry of permeability, porosity, and, most importantly, channelized sand bodies and low-permeability barriers that dictate subsurface fluid-flow pathways [5,12].

To robustly characterize subsurface formations, the Flow Zone Indicator (FZI) has emerged as a valuable petrophysical metric, integrating permeability and porosity into a single spatially explicit indicator that describes fluid-flow capacity within geological formations [13,14,15,16] (Amaefule et al., 1993; Jin et al., 2022). However, in complex fluvial and turbidite systems such as the Szolnok Sandstone, it is the explicit geometry, connectivity, and separation of sand-rich channels and intervening low-permeability units that most directly control fluid-flow behavior, thermal breakthrough, and UTES performance [17,18]. While FZI has been historically employed within hydrocarbon reservoir analysis for enhanced oil recovery, and has proven highly effective at delineating subsurface heterogeneity into distinct permeability zones, a more targeted modeling of channel systems offers even clearer identification of both optimal regions for fluid injection and production, as well as the natural low-permeability barriers that separate them [19,20,21].

Despite its established utility in reservoir characterization, the application of explicit channel system modeling, alongside or beyond FZI analysis, for UTES system optimization—identifying injection-production well sites and assessing barrier quality for thermal breakthrough prevention—remains comparatively underdeveloped [10]. Focusing on the three-dimensional architecture of channelized sands and their separation by low-FZI barriers represents a critical advance in targeting, designing, and optimizing UTES operations. Addressing this knowledge gap presents a valuable opportunity for advancing UTES technology by precisely targeting and optimizing subsurface geological conditions [10,15,22,23].

This research builds upon FZI analysis by developing a structured, data-driven methodology that explicitly models and predicts the three-dimensional architecture of sand-rich channel systems in the Szolnok Sandstone, with direct relevance for optimizing UTES site selection and operational planning [24]. Leveraging comprehensive petrophysical datasets—including detailed records of permeability, porosity, lithological characteristics, and historical well logs—the study employs advanced geostatistical techniques such as Ordinary Kriging and Sequential Gaussian Simulation (SGS)[25]. These methods enable the generation of high-resolution, three-dimensional maps that reflect generalized heterogeneity; further, it clearly delineates the geometry and connectivity of channelized sands and their separation by low-permeability zones [23]. While FZI-based mapping remains a valuable tool, this approach advances reservoir characterization by resolving the spatial distribution of channel bodies that most directly govern UTES performance [19,21,24].

By clearly identifying the spatial extent and continuity of sand channels suitable for injection and production wells, and mapping the intervening low-FZI barriers, this approach directly addresses the critical challenge of thermal breakthrough [24]. The explicit assessment of channel geometry, barrier thickness, continuity, and spatial arrangement allows for targeted well placement that leverages natural geological compartmentalization, enhancing thermal storage efficiency and long-term system sustainability [26]. Quantitative criteria such as barrier effectiveness and spatial distribution are thus evaluated to support robust thermal storage operations [4,15,27].

Additionally, this study applies a data-driven machine learning workflow, integrating XGBoost regression with advanced spatial feature engineering and thin-plate spline correction, to explicitly predict the three-dimensional presence, continuity, and separation of sand-channel systems within the reservoir [28]. By interpolating well log-derived properties using Inverse Distance Weighting (IDW) and correcting residuals with Radial Basis Function (RBF) interpolation, the approach captures both local heterogeneity and large-scale channel connectivity [19]. This method enables the efficient evaluation of numerous well-placement scenarios, providing rigorous uncertainty quantification and detailed sensitivity analysis—key for robust strategic planning and rapid deployment of UTES systems in channelized settings [21]. Such a hybrid modeling strategy substantially reduces computational complexity while maximizing spatial accuracy and interpretability, consistent with best practices in the field [20,21,28].

Hungary represents an ideal context for this investigation due to its dense network of depleted hydrocarbon wells, moderate-to-high geothermal gradients (around 70°C initial reservoir temperatures), and close proximity to urban centers with significant district heating demands. The presence of extensive Szolnok Sandstone Formation with its channel systems and comprehensive historical well data further enhances Hungary’s suitability for channel-focused UTES modeling (Figure 1). Moreover, Hungary’s commitment to decarbonization and energy security, combined with decades of subsurface data from hydrocarbon exploration, makes it an exceptional testbed for detailed geostatistical analysis and machine-learning-driven reservoir optimization [29,30,31].

This research systematically characterizes the 3D distribution and connectivity of sand channels—supported by FZI—to identify optimal injection and production well sites for UTES. For the first time, we apply a hybrid machine-learning approach integrating XGBoost regression with spatial residual correction to explicitly predict channel architecture and barrier distribution specifically within the Szolnok Formation. Emphasizing the mapping of high-permeability channels and low-FZI barriers, the approach mitigates thermal breakthrough risks, enhancing thermal isolation. Integrating numerical simulations with supervised machine-learning techniques optimizes system design under geological uncertainty, providing a reproducible framework for scalable, efficient UTES systems aligned with global decarbonization and energy resilience goals [2,3].

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.

Channelization: The development of laterally and vertically connected, sand-rich channel bodies within sedimentary reservoirs, which strongly control subsurface fluid-flow pathways and are critical for optimizing UTES performance and preventing thermal breakthrough.

Flow Zone Indicator (FZI): A petrophysical metric combining permeability and porosity into a single value that characterizes the fluid-flow capacity of a rock interval, widely used for identifying reservoir heterogeneity and mapping high-permeability channels and low-permeability barriers.

XGBoost: An advanced, scalable machine learning algorithm based on gradient-boosted decision trees, used for regression and classification tasks in this study to predict the three-dimensional distribution of sand channel systems from well log data.

IDW (Inverse Distance Weighting): A spatial interpolation technique that estimates unknown values at unsampled points based on the values of nearby observations, with closer points weighted more heavily; used here for spatial feature engineering of well log data.

RBF (Radial Basis Function) Interpolation: A mathematical method for smoothly interpolating values across multidimensional space; in this study, thin-plate spline RBF correction was used to improve the spatial accuracy of model residuals.

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

Gaussian Process Emulator: A statistical machine learning model used to create fast, probabilistic surrogate models of complex systems, enabling uncertainty quantification and efficient scenario evaluation.

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

Thermal Breakthrough: The phenomenon where injected hot or cold fluid reaches the production well too quickly, reducing system efficiency and potentially shortening operational lifetime; a key risk in UTES that can be mitigated by natural or engineered barriers.

Hydrogeological Model: A numerical or conceptual model simulating groundwater flow and heat transport based on geological and hydraulic properties, used to inform site selection and UTES system design.

Hot Well: A well used to inject and store heated water during UTES charging phases and extract thermal energy during production phases.

Cold Well: A complementary well used for injecting cooler water during UTES extraction or extracting cold water during charging, essential for managing thermal gradients and system balance.

Permutation Importance: A machine learning method for quantifying the importance of individual input features by measuring the decrease in model performance when feature values are randomly shuffled.

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

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

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

Overpressure: Subsurface pressure exceeding hydrostatic pressure, often due to geological compaction, tectonic forces, or fluid generation, affecting reservoir management.

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.

3. Background and Regional

3.1. Geological and Hydrogeological Setting

The Pannonian Basin, located in East-Central Europe, features complex geological structures, including subsided sedimentary basins separated by horst-like uplifted blocks [24,32]. The basement primarily comprises Paleozoic metamorphic rocks, occasionally overlain by Mesozoic carbonate formations, which act as significant aquifers in various locations, providing favorable conditions for geothermal exploration and thermal storage projects [30,31,33]. Within this intricate geological setting, the Szolnok Formation in southern Hungary represents an essential Lower Pannonian stratigraphic unit, widely recognized for its considerable reservoir potential and suitability for geothermal applications [33,34].

The Szolnok Formation predominantly consists of fine-grained sand turbidites interbedded with layers of siltstone and claystone, deposited within turbidite fan environments that developed partly in prodelta regions and partly in deeper basin settings of the extensive Lake Pannon basin during the Late Miocene (Lower Pannonian) [24,35,36]. This formation exhibits substantial variations in thickness, typically ranging from a few meters up to about 50 meters, reflecting complex depositional processes involving channelized turbidity flows and lobe systems sourced from adjacent paleo-deltas [37]. These depositional patterns significantly influence reservoir heterogeneity, directly impacting the effectiveness and optimization of geothermal energy storage operations, particularly for Aquifer Thermal Energy Storage (ATES) and Underground Thermal Energy Storage (UTES) systems [10,38,39,40].

Typically, the Szolnok Formation directly overlies thick marls and claystones of the Endrőd Formation, serving as regional seals, and is subsequently capped by prodeltaic silty clay sediments of the Algyő Formation [44]. This stratigraphic stacking results in substantial overpressure conditions, significantly influencing fluid dynamics and reservoir management strategies [34,45,46,47].

3.2. Hydrodynamic Systems and Pressure Regimes

The hydrodynamic environment of the Szolnok Formation is dominated by distinct flow regimes: shallow intervals influenced by gravity-driven groundwater systems and deeper intervals characterized by significant overpressure [24,30]. This overpressure is primarily generated by tectonic compaction, sediment loading, and gas formation from organic maturation processes, creating challenging conditions for subsurface fluid management. Pressure-depth relationships in the Szolnok Formation reveal significant super-hydrostatic pressures, often exceeding hydrostatic gradients by several tens of megapascals [45,47].

3.3. Reservoir Properties and Geothermal Potential

Reservoir characterization of the Szolnok Formation indicates favorable geological properties suitable for thermal energy storage [32]. Effective porosity within the formation typically ranges between approximately 15% to 17%, with permeabilities generally varying from 10 to 270 millidarcies (mD) [34,48,49]. Hydraulic conductivity, derived from permeability estimates, extends from approximately 11 to 2500 mmD, supporting considerable fluid flow capabilities and reservoir connectivity, essential for successful geothermal operations [32,50,51].

Detailed reservoir quality analysis highlights a moderate net-to-gross sand content typically between 50% and 70%, representing channelized sand bodies interspersed with finer-grained sediments such as siltstones and claystones, which act as natural intra-reservoir barriers or baffles [35,52]. The formation’s Flow Zone Indicator (FZI) typically exceeds values of 5.0, further emphasizing the excellent fluid flow potential for geothermal injection and production operations [15,21].

Geothermal gradients in the region of the Szolnok Formation are notably high, averaging approximately 50°C/km, reflecting relatively thin crustal conditions beneath the Pannonian Basin [47,53]. Consequently, reservoir temperatures are commonly around 70°C at economically viable depths of approximately 1,500–1,800 meters, with deeper intervals at around 2,000–2,500 meters potentially reaching temperatures as high as 90–120°C [22,30,34]. Such conditions strongly favor geothermal energy storage and district heating applications [5,22].

3.4. Hydrocarbon History and Well Infrastructure

The Szolnok Formation has experienced extensive hydrocarbon exploration and production since its discovery in the mid-20th century [32,54,55]. Numerous wells penetrate this formation, providing extensive geological, reservoir, and operational data critical for subsequent geothermal assessments and thermal energy storage implementation [5,56]. Typically, wells penetrating the Szolnok Formation in the studied region have employed telescopic casing designs ranging from approximately 340 mm diameter at the surface to 140–178 mm at reservoir depths [34,57].

The existing extensive well infrastructure, predominantly from past hydrocarbon production, offers significant economic advantages for repurposing wells into geothermal and seasonal thermal storage wells [22,38,58]. Although these wells vary in mechanical integrity and may require selective workover or recompletion, their presence substantially reduces initial investment requirements, making the implementation of geothermal projects more financially viable [32,59].

3.5. Relevance to Seasonal Heat Storage

The geological, hydrodynamic, and geothermal conditions of the Szolnok Formation in southern Hungary make it particularly suitable for Aquifer Thermal Energy Storage (ATES), including advanced High-Temperature ATES (HT-ATES) applications [24]. Its combination of favorable porosity (15–17%), significant permeability (up to 270 mD), and excellent hydraulic conductivity (up to 2500 mmD) supports efficient fluid movement and heat exchange essential for effective thermal storage and recovery [24,31].

The geological, hydrodynamic, and geothermal conditions of the Szolnok Formation in southern Hungary make it particularly suitable for Aquifer Thermal Energy Storage (ATES), including advanced High-Temperature ATES (HT-ATES) applications [38,59]. Its combination of favorable porosity (15–17%), significant permeability (up to 270 mD), and excellent hydraulic conductivity (up to 2500 mmD) supports efficient fluid movement and heat exchange essential for effective thermal storage and recovery [34,48,51].

The formation’s inherent stratigraphic heterogeneity, with sand-rich channel systems surrounded by finer-grained sediments, provides naturally compartmentalized conditions [24,35,52]. These heterogeneous reservoir characteristics create ideal scenarios for strategically positioning injection and extraction wells, leveraging high-permeability channelized sands for thermal storage, while low-permeability regions serve as natural barriers preventing premature thermal breakthrough [11,12]. Optimal well-spacing and operational strategies, therefore, rely on detailed geological characterization and modeling of reservoir heterogeneity.

Furthermore, the proximity of depleted hydrocarbon wells to district heating infrastructures in southern Hungary presents advantageous conditions for efficiently integrating thermal storage solutions into existing heating networks [22,60]. However, careful management of the formation’s substantial overpressure conditions and stratigraphic complexities remains critical to ensure sustainable and efficient thermal energy storage operations [54,57].

Collectively, leveraging the Szolnok Formation’s geological insights, favorable reservoir conditions, existing well infrastructure, and high geothermal potential positions it as an outstanding candidate for developing sustainable, scalable seasonal thermal energy storage systems [31,56].

4. Materials and Methods

4.1. Methodological Framework

The methodological framework of this study began with comprehensive data curation and preparation, involving extensive modeling, cleaning, and refinement of existing hydrocarbon field datasets. Using RockWorks software, an initial Flow Zone Index (FZI) model was constructed based on reliable porosity and permeability estimates derived from detailed reservoir characterizations [61,62]. Subsequently, a machine-learning workflow was developed, leveraging well-log data to predict spatial distributions of the FZI (Figure 2). XGBoost regression was chosen for its robustness and efficiency in capturing complex geological patterns, while IDW and RBF methods effectively addressed spatial continuity and residual smoothing. The predicted FZI outputs from the initial machine-learning model served as input training data for a HistGradientBoostingRegressor [14,16,63], which was evaluated on a withheld dataset representing approximately 30% of the study area. The final predictive model demonstrated strong performance, achieving an R² of approximately 0.92, a mean absolute error (MAE) of 0.17, and a root mean squared error (RMSE) of 0.24 [14].

4.2. Data Collection and Data Preparation

The data collection phase of this study involved systematically compiling detailed subsurface geological and petrophysical data from the Szolnok Formation in southern Hungary, utilizing comprehensive historical records from 128 hydrocarbon exploration and production wells. The primary datasets included well-log measurements, core analyses, and stratigraphic interpretations, providing key information such as formation tops and bases, thickness variations, sand content distributions, effective porosity measurements, permeability data, and hydraulic conductivity values (Figure 3) [64].

4.3. Data Modelling

The data collection process for this study involved acquiring detailed subsurface data from the Szolnok Formation in southern Hungary, leveraging information derived from 128 hydrocarbon exploration and production wells. The dataset primarily comprised stratigraphic and reservoir property data, including top and base surfaces, thickness, sand content, effective porosity, effective pore volume, permeability, and reservoir quality indices (Figure 4) [32]. Stratigraphic surfaces (top and base) were estimated using kriging with a trend model [65], while sand content, porosity, and permeability were modeled using both Inverse Distance Weighted (IDW) interpolation—specifically optimized to emphasize lateral reservoir continuity—and Ordinary Kriging techniques to enhance spatial prediction reliability [66]. The effective porosity within the formation ranged between 15% and 17%, while permeability values ranged from 10 to 270 mD, corresponding to hydraulic conductivity values of approximately 11 to 2500 mmD [67]. Reservoir Quality Index (RQI) and Flow Zone Index (FZI) models were constructed by combining these spatially interpolated porosity and permeability distributions, yielding average FZI values exceeding 5.0, indicative of excellent fluid-flow capabilities within the channelized sand bodies [68,69].

4.4. Flow Zone Index Modelling

The Flow Zone Index (FZI) modelling for the Szolnok Formation involved integrating spatially interpolated porosity (ϕ) and permeability (k) datasets to quantitatively characterize and visualize reservoir fluid-flow capacities [16]. FZI was calculated using the standard reservoir engineering equation [13,70,71,72]:

$$FZI={\displaystyle \frac{RQI}{NPI}}={\displaystyle \frac{0.0314\sqrt{\frac{k}{\varphi }}}{\frac{\varphi }{1-\varphi }}}$$

where RQI is the Reservoir Quality Index, k is permeability (mD), and ϕ is effective porosity (expressed as a percentage). Spatial interpolation was conducted using both Inverse Distance Weighted (IDW) methods, optimized specifically to emphasize lateral continuity, and Ordinary Kriging supported by directional variogram analyses (Figure 5) [73,74]. These geostatistical approaches were selected to capture regional trends and localized heterogeneities effectively, ensuring accurate and reliable delineation of subsurface flow zones.

The resulting FZI models, generated independently using both IDW and Kriging methods, consistently indicated highly favorable fluid-flow conditions, particularly within channelized sand bodies. Specifically, the IDW-derived FZI values exceeded 4.5, while Kriging-derived estimates indicated slightly higher values, generally exceeding 5.0. The consistently high FZI values reflect excellent reservoir connectivity and permeability within these turbiditic channel systems, highlighting their suitability for fluid injection and extraction operations [75]. By clearly delineating areas of optimal flow capacity, the FZI modelling results provided critical guidance for strategic site selection and well placement decisions aimed at maximizing thermal storage efficiency and preventing premature thermal breakthrough [76].

4.5. Machine Learning Setting

The training dataset was meticulously prepared for modeling the Flow Zone Index (FZI) in three-dimensional space, integrating rich petrophysical features and rigorous spatial preprocessing. Initially, interval-level data extracted from multiple well logs were processed to obtain geospatially accurate samples characterized by depth (z) and crucial petrophysical attributes, including permeability (PERM), porosity (PHI), shale volume (VSHA), sandstone volume (VSND), effective thickness (HEFF), and total thickness (thickness) (Figure 6).

An essential preprocessing step was the targeted removal of outliers from both well intervals and grid node observations, employing robust statistical techniques: Z-score (±3.0) and Interquartile Range (IQR; 1.5 × IQR) filtering. This careful cleaning significantly improved data quality, stability, and model accuracy [77].

Geometrically, the spatial domain was defined and optionally refined through the application of a buffer around core wells, effectively isolating a central region for focused analysis. Within this spatially constrained area, a systematic 80/20 train-test split was performed on the available 3D grid nodes, ensuring a reliable validation framework for predictive performance [78,79].

Subsequent feature engineering leveraged Inverse Distance Weighting (IDW) interpolation to generate robust, spatially contextualized attributes for each grid node. Specifically, petrophysical features from the nearest well intervals were interpolated, forming IDW-derived features (PERM_idw, PHI_idw, VSHA_idw, VSND_idw, HEFF_idw, and thickness_idw). Further feature augmentation included logarithmic transformations (e.g., log_PERM_idw) to normalize skewed data distributions and derived ratio features (e.g., permeability-to-porosity ratio PERM_div_PHI_idw) to capture meaningful petrophysical relationships [80].

A powerful hybrid modeling approach was employed, combining XGBoost regression (trained via 5-fold cross-validation) with a spatially explicit residual correction using Thin-Plate Spline interpolation [81]. This integrated methodology significantly enhanced prediction accuracy and spatial consistency.

Comprehensive quality control and diagnostics were systematically generated, including scatter plots, residual analyses, depth slices, feature importance rankings, and geospatial prediction maps, facilitating robust model validation and interpretability [78]. Collectively, these engineered and spatially informed features provided a robust foundation for accurate, reliable, and interpretable three-dimensional FZI predictions [82].

4.6. Machine Learning Process

The machine learning process implemented in this study integrated rigorous data preprocessing, spatial feature engineering, and robust predictive modeling to accurately forecast three-dimensional sand-channel systems indicated by the Flow Zone Index (FZI) [68]. Initially, interval-level petrophysical well data, including permeability and porosity (input features), were systematically loaded and filtered for outliers using robust Z-score statistical methods to ensure the reliability of model inputs [83]. Permeability values in the training set ranged from approximately 0.01 mD to 980 mD, while porosity values spanned from approximately 5% to 28%. Subsequently, the spatial modeling domain was defined using a polygonal boundary shapefile, optionally refined by buffering to retain core, most representative reservoir areas [24]. Advanced spatial interpolation techniques, specifically Inverse Distance Weighting (IDW), generated spatially informed predictive features at each 3-D grid node, further augmented by composite attributes and log-transformed petrophysical features designed to capture geological relationships [84]. The complete workflow, including each critical step from data loading to final predictions and evaluation, is illustrated in Figure 7.

In the second stage, an optimized XGBoost regression model was rigorously trained using an 80/20 train-test split of the grid nodes, and validated via a comprehensive 5-fold cross-validation to assess predictive accuracy [78]. The output predictions of the model were subsequently improved by applying spatially explicit residual corrections using thin-plate spline interpolation (RBF), significantly enhancing spatial consistency [85]. Comprehensive evaluation metrics (R², MAE, RMSE) were calculated from the independent test set, confirming the model’s high accuracy (test R² = 0.92) [86]. Ultimately, this integrated modeling workflow yielded robust, accurate, and interpretable three-dimensional predictions of sand-channel distribution, directly informing reservoir characterization and operational decision-making for optimized UTES well management [24].

4.7. Sensitivity Analysis

A sensitivity analysis was conducted by systematically evaluating combinations of key modeling parameters, specifically the number of neighbors used in inverse distance weighting (idw_k_neigh), radial basis function interpolation (rbf_neigh), the XGBoost learning rate, and the Z-score threshold (z_thresh) applied for outlier filtering [87]. A total of 50 parameter combinations were randomly sampled from a structured parameter grid to assess their impact on predictive performance, measured by cross-validation R², test R², mean absolute error (MAE), and root mean squared error (RMSE) [88]. The analysis revealed that the model was particularly sensitive to the learning rate and the spatial complexity controlled by idw_k_neigh and rbf_neigh. The optimal predictive performance (test R² = 0.971, MAE = 0.093, RMSE = 0.178) was achieved using moderate spatial complexity (idw_k_neigh=7, rbf_neigh=40), a relatively higher learning rate (0.05), and a stringent Z-score threshold (3.0) [89]. The detailed results of this optimal parameter combination and corresponding performance metrics are presented in Table 1. These findings emphasize the critical balance required in parameter selection, highlighting that intermediate values for spatial parameters combined with higher learning rates effectively improve predictive robustness while mitigating potential overfitting or underfitting scenarios.

4.8. Model Calibration and Validation

The model calibration and validation involved analyzing the relationships between the normalized predicted FZI values (FZI_norm) and key petrophysical indicators, specifically normalized sand content (VSND_norm) and normalized shale content (VSHA_norm). Correlation analysis showed a statistically significant positive correlation between FZI_norm and VSND_norm, with Pearson’s correlation coefficient (r) of 0.333 (p=0.000461) and Spearman’s rank correlation coefficient (ρ) of 0.360 (p=0.00014). Conversely, a statistically significant negative correlation was observed between FZI_norm and VSHA_norm, yielding a Pearson’s r of -0.384 (p=4.45e-05) and Spearman’s ρ of -0.389 (p=3.41e-05). These findings confirm that the predicted normalized FZI is consistent with geological expectations, positively correlating with sandy, permeable zones and negatively with shaly, lower-quality intervals (Figure 8) [90].

5. Results

5.1. FZI Prediction Result

The three-dimensional FZI prediction model demonstrated high accuracy and stability across all depth slices in the reservoir. Aggregated performance metrics reveal that the mean absolute error (MAE) averaged just 0.056, while the mean root mean squared error (RMSE) was approximately 0.17, with the majority of predictions exhibiting errors well below these values. The distribution of MAE and RMSE across individual slices is tightly clustered, with 75% of grid predictions falling under 0.055 MAE and 0.19 RMSE, and only a handful of locations reaching a maximum MAE of 0.43 and RMSE of 0.87. These statistics underscore the overall precision and consistency of the model, ensuring that both local and global FZI variations are faithfully captured (Figure 9) [91].

The predicted Flow Zone Index (FZI) values across the three-dimensional grid display a high degree of concordance with the observed data, as demonstrated by the distribution of residuals between predictions and measurements. The majority of grid points exhibit very small residuals, with the largest concentration falling in the lowest error bins—thousands of points have absolute residuals less than 0.01, indicating that the model’s predictions are extremely close to actual FZI values for most of the reservoir [91]. The gradual increase in bin counts with increasing error, followed by a sharp decline at higher error levels, further confirms that large deviations are rare and the prediction error is dominated by minor discrepancies.

5.2. Residual Distribution and Model Robustness

A closer examination of the residual histogram reveals that over 90% of the prediction errors are less than 0.2 in absolute value, with the most populated bins spanning error ranges from 10⁻⁶ to approximately 0.2 [92]. Only a small fraction of points have larger residuals, with very few exceeding 1.0. This skewed, tightly centered error distribution reflects the robustness of the machine learning workflow, which includes extensive outlier filtering, careful spatial feature engineering, and a two-step modeling process (XGBoost regression followed by RBF residual correction)[93]. The results suggest the model successfully captures both the spatial heterogeneity and the petrophysical complexity of the reservoir, yielding predictions with minimal systematic bias (Figure 10).

Overall, the prediction outcome demonstrates both accuracy and stability across the study area, with negligible evidence of overfitting or model drift. The extremely low and sparsely distributed large residuals imply that prediction failures, when they do occur, are isolated and do not compromise the reliability of the overall FZI map [16,94]. Such a result is crucial for reservoir characterization workflows, as it provides a high-resolution, spatially continuous prediction surface with quantified uncertainty, supporting confident operational decisions and further geological analysis.

6. Discussion

6.1. Interpretation of Predicted 3D FZI Clusters

The predicted three-dimensional FZI clusters, classified into two distinct categories—Low-FZI (<4.5) and High-FZI (≥4.5)—clearly delineate reservoir units into sand-rich channel systems and surrounding low-quality lithologies (Figure 10) [16,95,96]. The Low-FZI cluster includes 363 grid nodes and 72 wells, characterized by an average FZI value of 3.85 (std=0.67), typically indicating fine-grained or cemented lithologies with reduced permeability and limited fluid transmissivity. In contrast, the High-FZI cluster comprises 366 grid nodes and 60 wells, exhibiting a higher average FZI of 4.99 (std=0.41), with maximum values up to 8.47. These high-FZI zones strongly suggest well-connected, sandy channel systems favorable for fluid flow [68].

Wells that directly penetrate the High-FZI channels are shown in red; however, some wells located within high-FZI zones remain marked in black, indicating that although they appear spatially associated, they did not penetrate the channelized system vertically (Figure 10C). This highlights the importance of precise 3D delineation for accurately identifying reservoir connectivity and optimizing well placement [20,97].

The clear separation between the Low and High-FZI clusters emphasizes the model’s ability to accurately distinguish between sand-rich channels and surrounding lower-quality lithologies, reflecting significant petrophysical contrasts [98]. The minimal overlap and tightly constrained standard deviations within each cluster confirm the robustness of the classification, indicating internally consistent reservoir properties. The balanced spatial distribution of nodes and wells across both clusters highlights the geological complexity and heterogeneity of the reservoir [99]. Crucially, the distinct delineation of High-FZI clusters provides essential insight for reservoir management, enabling targeted placement of wells within high-permeability channel systems, thus optimizing fluid flow, enhancing energy recovery, and reducing the risk of premature thermal breakthrough [99].

6.2. Subsurface Complexity of High-FZI Regions in 3D

The horizontal (depth-mean) analysis of the predicted FZI clusters identifies distinct high-FZI regions (FZI ≥ 4.5) distributed across the reservoir, but the three-dimensional visualization highlights considerably greater complexity and detail [100]. In the 3D view, the high-FZI zones form vertically stacked, laterally extensive, channel-like features, consistent with known channel deposition characteristics of the Szolnok Formation (Figure 11 Right) [101]. Specifically, out of 60,954 analyzed grid nodes, extensive regions exhibit FZI values exceeding 4.5, closely aligning with interpreted sandy channel systems [102]. Notably, 29 out of 136 wells directly penetrate these high-FZI channels (highlighted in red), while the remaining 107 wells fall outside these zones (displayed in black). This clear delineation provides valuable spatial context for reservoir characterization, reflecting geological controls on fluid-flow pathways [100].

This detailed three-dimensional understanding of the high-FZI zones offers substantial advantages for the design and optimization of UTES (Underground Thermal Energy Storage) cold and hot well duplets [103]. The naturally channelized high-FZI zones, indicative of enhanced permeability and reservoir quality, can be strategically utilized to maximize heat transfer and storage efficiency while simultaneously mitigating the risk of thermal breakthrough. By placing injection (cold) and production (hot) wells in separate, isolated channel segments, operators can effectively use the natural geologic architecture as barriers, minimizing unintended thermal communication [104]. Consequently, leveraging the inherent channel geometry of the Szolnok Formation supports optimal placement of duplet wells and also enhances overall UTES system performance and sustainability [103].

6.3. Enhancing Decision-Making for UTES Site Selection

Direct well data-driven prediction of channel systems within the Szolnok Formation significantly enhances decision-making for UTES (Underground Thermal Energy Storage) site selection and the management of cold and hot wells in depleted hydrocarbon reservoirs [24]. By utilizing high-resolution, spatially continuous models that clearly identify channelized sand bodies based on actual well log data, operators can precisely target highly permeable zones ideal for thermal energy injection and extraction. Accurate delineation of channel geometries enables strategic placement of cold and hot well duplets, effectively utilizing natural geological boundaries to minimize the risk of premature thermal breakthrough [23]. The model’s ability to reliably capture the spatial distribution and connectivity of channelized zones reduces subsurface uncertainty, ensuring that well planning aligns closely with geological reality. Ultimately, this channel-focused predictive approach promotes more efficient, sustainable, and cost-effective thermal energy storage operations within repurposed hydrocarbon reservoirs [24].

6.4. Generalizing Channelization Prediction for the Szolnok Formation

This model offers substantial potential for predicting channel systems in other basins containing the Szolnok Formation, particularly if comparable well log datasets are available [24]. Given the widespread occurrence and consistent depositional characteristics of the Szolnok Formation across the region, the established relationships between well log responses and channel presence are likely transferable to areas lacking detailed channelization data [105]. However, the methodology may require recalibration or adaptation in regions with significantly different depositional environments, sediment composition, or structural complexities. By applying the same preprocessing and feature-engineering workflows to well logs from the new basin, operators can rapidly generate high-resolution predictions of channel geometries and connectivity [19]. This approach enables efficient identification of sandy channel bodies ideal for targeted UTES, geothermal, or hydrocarbon development, even in locations where direct geological observations are limited [106]. Ultimately, the ability to generalize channel predictions significantly reduces uncertainty in reservoir characterization, guiding the optimal placement of wells and enhancing data-driven decisions across the broader Pannonian Basin [24].

6.5. Managing Model Uncertainty and Reproducibility

Although our model for predicting subsurface channel systems demonstrates high accuracy, some uncertainty inevitably remains due to potential differences in well log calibration, local geological variations, and the representativeness of training data—particularly when transferring the model to other basins or datasets [107,108]. While this study specifically focused on one oil field in southern Hungary, several other oil fields have historically produced from the Szolnok Formation, potentially benefiting significantly from the methodological framework developed here [109]. Applying this structured, data-driven approach to these additional sites could enhance reservoir characterization, optimize thermal energy storage, and facilitate effective site selection and management across the broader region [24]. To maximize transparency and enhance its applicability to new geological settings, I have shared the fully documented code and methodology, including the trained XGBoost model, RBF residual correction model, and feature engineering schema.[110] These files ensure consistent and reproducible predictions of channel systems by providing the complete workflow, including feature preparation, model application, and residual correction [111]. Users can leverage these saved artifacts to quantify uncertainty, verify each modeling step, and effectively adapt the channel prediction workflow to other regions with similar log data, ensuring robust and geologically meaningful interpretations [111,112].

6.6. Recommendations and Future Work

Future work should focus on improving the model’s predictive accuracy and transferability for identifying sand-rich channel systems across different geological settings by incorporating diverse training datasets from multiple basins. Developing a standardized well-log preprocessing pipeline—including consistent calibration checks and unified feature engineering—will reduce biases when applying the model to new areas. Furthermore, integrating seismic data could provide independent structural validation, while core analyses and previous hydrocarbon production histories would offer additional calibration points, substantially enhancing the resolution and confidence of channel delineations. Employing advanced uncertainty quantification methods, such as probabilistic modeling or ensemble predictions, will also provide valuable assessments of prediction reliability. Lastly, collaborative data sharing and integration of these validation sources will significantly enhance model validation and practical utility.

7. Conclusion

This research successfully demonstrated a structured, data-driven methodology for accurately predicting the three-dimensional distribution of sand-rich channel systems, indicated by high Flow Zone Index (FZI ≥ 4.5), within the Szolnok Formation in southern Hungary. Leveraging extensive petrophysical data from 128 wells, the study clearly delineated high-permeability sandy channels from surrounding lower-quality intervals characterized by low FZI. The resulting spatial models distinctly identified interconnected channelized sands separated by continuous, low-permeability zones, effectively forming natural barriers crucial for mitigating premature thermal breakthrough in Underground Thermal Energy Storage (UTES) systems.

Although this research focused on one representative oil field in southern Hungary, numerous other fields historically producing from the Szolnok Formation could similarly benefit from adopting this predictive framework. The application of this methodology across these additional sites would enhance reservoir management practices, reduce geological uncertainty, and significantly improve the operational effectiveness of UTES installations regionally.

The integrated modeling workflow, which included rigorous outlier filtering, spatial feature engineering using Inverse Distance Weighted (IDW) interpolation, and predictive modeling via XGBoost regression enhanced with residual corrections, achieved high predictive accuracy. Validation results revealed a robust model performance with an R² of 0.92, a mean absolute error (MAE) of approximately 0.056, and a root mean squared error (RMSE) of 0.17. The precision and reliability of these predictions ensure confident strategic planning and optimal well placement, significantly enhancing the operational sustainability and efficiency of thermal storage.

By explicitly identifying regions of high FZI suitable for hot and cold well placement and delineating effective low-FZI barriers, the methodology directly addresses the critical challenge of thermal breakthrough, enhancing UTES system longevity and thermal recovery efficiency. Future research should focus on further refining model accuracy and generalizability through standardized data preprocessing, integration with seismic surveys, and advanced uncertainty quantification. Ultimately, this robust, reproducible methodology substantially advances reservoir characterization, supporting the global expansion of UTES technology for sustainable, scalable, and economically viable thermal energy storage solutions.