Denise Degen

Modelling the dynamic response of reservoirs to fluid injection often relies on fully coupled thermo-hydro-mechanical models, which provide a high-fidelity representation of the governing physical processes. These models support operational and design decisions under significant geological and parametric uncertainties. However, their high computational cost severely limits their applicability in large-scale statistical analysis and thus limiting the potential to account for these uncertainties. We present Platypus, a Python package that bridges full-order process-based reservoir simulations and surrogate-based analysis. The package provides automated, reproducible workflows covering the full pipeline from parameter sampling and simulation interfacing, through ensemble curation and data management, to surrogate model construction and global sensitivity analysis. Surrogate models are built via the non-intrusive reduced-basis method to deliver computationally efficient, physics-consistent approximations of the transient reservoir response, enabling statistical analysis that require a large number of evaluations. The package is demonstrated on two geothermal applications: a synthetic reservoir benchmark and the Groß Schönebeck geothermal site. These examples illustrate the package's ability to handle time-dependent coupled hydrothermal and thermo-hydro-mechanical processes, manage large simulation ensembles, and perform global sensitivity analyses enabled by surrogates. Platypus is designed to be modular and extensible, providing a flexible tool for ensemble-based reservoir modelling workflows across a range of subsurface applications.

Reliable estimates of subsurface heat resources are essential for reducing financial risk in early-stage geothermal development. Regional assessments typically rely on the volumetric Heat In Place (HIP) method, which estimates stored thermal energy from reservoir volume, effective heat capacity, and temperature. This approach commonly assumes purely conductive heat transport and usually represents structural features such as faults indirectly through recovery factors. In reality, faults can host convective fluid flow that redistributes heat and creates localized thermal anomalies, potentially biasing resource estimates in structurally complex settings. However, this effect has not been systematically quantified. By using thermo-hydraulic simulations of a 25x25\si{\ \km\squared} fault-dominated carbonate reservoir system in North Rhine-Westphalia, Germany, we investigate how fault-hosted convection alters local HIP estimates. At elevated fault transmissivities, local HIP anomalies differ up to 60\% relative to the conductive reference case. At moderate transmissivities, the resulting perturbations are largely restricted to the kilometer-scale proximity of the fault. The simulations further demonstrate that the depth-dependent transmissivity distribution of fault damage zones strongly influences the magnitude and lateral extent of convection-induced thermal anomalies. In addition, convection self-organizes into large-scale patterns across closely spaced faults even without direct hydraulic connection between them. These results indicate that neglecting fault-controlled fluid flow can misrepresent geothermal potential in fault-dense regions. More broadly, our findings support the case for probabilistic workflows in early-stage resource assessment, particularly when fault networks are pervasive but poorly constrained.

Geomechanical-numerical modeling aims to provide a continuous description of the stress tensor within rock volumes by leveraging localized stress magnitude data for model calibration. This calibration process involves optimizing displacement boundary conditions to achieve the best possible alignment with in-situ stress measurements obtained from microhydraulic fracturing (MHF) and sleeve reopening (SR) test in boreholes providing magnitudes of the minimum and maximum horizontal stresses. However, the high cost associated with acquiring stress magnitude data often results in sparse and incomplete datasets, which can potentially hinder a robust model calibration.We investigate the relationship between calibration dataset size and stress field prediction quality by leveraging an exceptionally high-quality set of 45 in-situ stress magnitude data from MHF and SR tests performed in two boreholes in the Zürich Nordost region, a potential deep geological repository location in northern Switzerland. We present a statistical framework that analytically evaluates extensive model ensembles, where each realization uses a distinct subset of the entire stress magnitude dataset for the model calibration. This methodology quantifies the effect of data quantity on modeled stress distributions in lithological units with different elastic properties by efficiently computing stress states for multiple data combinations.By systematically comparing modeled stress magnitudes derived from incrementally larger calibration datasets, we identify the minimum number of stress magnitude data records needed to achieve a modeled stress range narrower than the predicted stress range resulting from rock stiffness variability. For Zürich Nordost, approximately 15 stress magnitude data records are sufficient to reach this threshold, beyond which additional data provides only minimal improvement to model predictions.Examination of the dataset revealed one measurement that deviated from the expected model, due to localized variations in the mechanical properties of the rock. Although this observation is geologically valid, it substantially distorts predictions when there are few calibration points available. However, its distorting effect weakens progressively as the size of the calibration dataset increases. Our analytical framework also enables the systematic detection of such anomalous records, which in turn affects the minimum number of data points required for model calibration.Our findings emphasize the importance of dataset size and composition in reducing uncertainties and help establishing a practical methodology for efficient data acquisition planning. These insights are particularly valuable for subsurface engineering applications that require reliable stress estimates, such as geothermal systems, carbon storage facilities, radioactive waste repositories and underground energy infrastructure.

Accurate and efficient modelling of geothermal reservoirs is important for sustainable energy production and for the reliable assessment of operational risks. Predicting thermo-hydraulic (TH) system evolution under varying injection and production scenarios remains computationally challenging, particularly when physical knowledge of the subsurface system is incomplete and observational data are sparse. High-fidelity finite-element simulators are typically used to provide physics-based predictions of coupled flow and heat transport governed by complex partial-differential equations (PDEs). Such full-order simulations are, however, often prohibitively expensive for real-time forecasting, which is essential, for instance, in the context of digital twins.Physics-based machine-learning (PBML) approaches, such as the non-intrusive reduced basis (NIRB) method address this challenge by constructing physics-consistent surrogate models that project full-order simulation outputs onto a low-dimensional subspace learned from representative snapshots. By retaining only the dominantbasis functions, the NIRB surrogate enables orders-of-magnitude speedup in parametric predictions while staying consistent with the physical transport mechanisms and structural assumptions on fracture networks encoded in the full-order model. Despite these advantages, classical NIRB surrogates are intrinsically limited to the physical regimes represented by the governing PDEs, and consequently by the training simulations. If the surrogate does not fully capture the observed system behaviour, it is important to detect and adapt to missing or misrepresented local physics revealed by observational data, such as unmodeled convective heat transport or flow channelling arising from fracture activation.To address this need, we propose a complementary residual-learning framework that augments a baseline NIRB surrogate with parameter-to-state maps of residual temperature and pressure fields learned by Kolmogorov-Arnold Networks (KANs). The residual, defined as the difference between observed data (or a synthetic reference solution) and the NIRB model prediction, is interpreted as a proxy for missing or misrepresented physics not explicitly captured by the baseline model. KANs represent mappings as sums of learned univariate functions and provide explicit access to the functional structure of parameter dependence. Thereby, KANs could act as interpretable discrepancy models by learning the residual between observations and NIRB predictions. By analysing the dominant functional families emerging in the learned residual, such as linear dependence characteristic of conduction-dominated regimes or exponential dependence associated with convection, KANs can provide diagnostic insight into missing thermo-hydraulic processes and their relevance across parameter regimes.We validate the proposed approach synthetically by comparing a conduction-only NIRB surrogate against synthetic reference observations generated with an advection–diffusion model. We expect that KAN-based residual learning both improves predictive accuracy and reveals clear functional signatures of missing convective physics, even when only pointwise information is available. As an outlook, we aim to apply this workflow to real geothermal case studies, where sparse temperature and pressure measurements are available at well locations. In such settings, functional-family learning of residuals offers a promising pathway to improve surrogate predictions and to enhance the physical interpretability of geothermal systems, ultimately supporting more reliable assessments of reservoir behaviour.

Predictions of the three-dimensional in-situ stress state are crucial for the site selection process for deep geological repositories for nuclear waste and their long-term safety. However, the geological configurations relevant to a potential siting region, such as fault offsets, salt intrusions, and intersecting sedimentary units, create increasing structural complexities. These complexities, especially in the form of intersecting lithological boundaries and interfaces, present significant challenges for the discretization using the Finite Element Method (FEM). In the FEM, lithological boundaries and interfaces are commonly modeled using spline-based representations, which are then geometrically approximated by finite elements. This process introduces an additional layer of geometrical approximations that can lead to discretization errors, mesh distortions, and the need for repeated geometry regeneration when testing different model scenarios.This study investigates whether the Isogeometric Analysis (IGA) as a discretization method can enhance and facilitate geometric fidelity and contribute to an automated modeling workflow. Unlike FEM, IGA employs the same spline basis functions (e.g., NURBS) for both the geometrical representation and numerical approximation. This direct application of splines eliminates the need for a geometry-to-mesh approximation step, allowing for an exact representation of both lithological boundaries and structural features, such as faults. The workflow for IGA differs from traditional FEM primarily in the preprocessing, solver implementation, and postprocessing stages: geometry is handled directly through control points, spline basis functions replace conventional shape functions, and the numerical solution is stored at these control points before being mapped back to the physical domain. While IGA does not necessitate a separate meshing step, refining the spline representation may still be required.To evaluate this approach, we begin with a three-layer benchmark model previously used in sensitivity analyses and introduce a fault that offsets the lithological layers. IGA is utilized to compute stress, strain, and displacement fields, with its performance compared to that of the FEM, focusing on the impact of the geometrical approximation. The results aim to illustrate how exact geometrical representation and spline refinement influence stress predictions, particularly in areas where faults or salt contacts create sharp geometrical variations.This work represents a significant advancement toward a more automated and reliable geomechanical modeling workflow. By reducing the need for manual geometrical regeneration and directly integrating spline-based representations into the analysis, IGA can streamline model scenario exploration and support more consistent gemechanical modeling for repository-scale studies.

A key prerequisite for reliable geoscientific process simulations is the calibration of uncertain model parameters against field observations. In practice, both measurements and simulation outputs are subject to uncertainty, arising from the observational errors, limited knowledge of material properties and inexact physical models. Bayesian inference provides a framework to explicitly acknowledge multiple sources of uncertainty by encoding modelling assumptions in prior distributions and updating them against observational data through the likelihood to obtain posterior estimates. However, applying Bayesian methods remains challenging in coupled multiphysical applications, including thermo-hydro-mechanical problems, as computational costs of repeated forward evaluations grow rapidly with model complexity. To address these limitations, we develop a hierarchical simulator for Bayesian calibrations that dynamically combines fast low-fidelity surrogate models with accurate high-fidelity finite-element simulations during the sampling stage. The core of the method stems from a fidelity-selection policy embedded directly in the probabilistic model, which transparently accounts for both surrogate-induced bias and the computational cost associated with high-fidelity simulations. We provide and compare several scenarios, that represent different optimization strategies for balancing posterior accuracy and computational efficiency. The resulting hierarchical Bayesian workflow is highly modular, and it can be coupled with external high-fidelity solvers through a unified forward interface and hence applicable to a wider range of geoscientific problems.

The subsurface of North-Central Europe has been shaped by repeated glaciations, which have altered pore pressure, temperature, and groundwater flow over geological timescales. Understanding these coupled thermo-hydraulic responses is essential for groundwater management, geothermal energy utilization, and subsurface storage applications such as CO₂ sequestration.In this study, we investigate the thermo-hydraulic response of a multi-layered subsurface model of North-Central Europe over a time period the last glacial maximum up to present-day, while considering the thermal and hydraulic feedback of the glacier dynamics on the distribution in time and space of pore pressure and temperature.Independent hydraulic and thermal simulations are constructed, followed by controlled in which porosity and permeability are systematically varied within selected stratigraphic units while all other parameters are held constant. The resulting transient pore pressure and temperature fields are analyzed to assess the relative roles of conductive and advective heat transport and to identify formation-specific controls on pressure dissipation and thermal redistribution.Based on the generated simulation ensemble, we further explore the development of physics-preserving, interpretable AI-assisted surrogate models and parameter estimation. These surrogates aim to efficiently reproduce key thermo-hydraulic responses while retaining physical consistency, thereby enabling rapid sensitivity analysis and uncertainty quantification in large-scale subsurface systems.

Geomechanical modeling aims to predict the 3D in-situ stress state of the Earth’s crust and to assess the stability of subsurface rock volumes for applications such as radioactive waste disposal, energy storage, or CO₂ geo-sequestration. However, model calibration typically relies on sparse in-situ stress magnitude data which are expensive to acquire, limited in spatial coverage, and may not represent stress conditions over larger rock volumes, away from the measurement sites. Here we present a probabilistic forward-calibration framework that uses the borehole failure interpreted from routinely acquired borehole-image logs as indirect stress data and formation integrity tests (FIT) to calibrate 3D geomechanical models.Our approach integrates four types of indirect stress observations: the occurrence of borehole breakouts (BO), drilling-induced tensile fractures (DITF), formation integrity tests (FIT), and the documented absence of both BO and DITF at micro-hydraulic fracturing (MHF) stations. Although these indirect data provide only upper and lower limits on the stress state, they offer the critical advantage of scanning the entire borehole trajectory with high resolution, yielding far more extensive spatial coverage than point measurements. The absence of borehole failure provides simultaneous upper and lower bounds on horizontal stress magnitudes, addressing a key limitation in previous approaches that struggled to constrain the maximum horizontal stress magnitude. We developed a forward uncertainty quantification framework that explores hundreds of thousands of model scenarios at each observation point using linear elastic principles and compares the agreement between predicted and observed stress indicators through a probabilistic assessment.In the Zürich Nordost siting region for a potential deep geological repository for radioactive waste in northern Switzerland, we leverage an exceptional stress magnitude dataset from two deep boreholes. This dataset comprises 30 high-quality microhydraulic fracturing tests and 15 dry sleeve reopening tests, accompanied by comprehensive borehole image logs and detailed laboratory measurements of Young's modulus and rock strength. Using the stress magnitude data alone to calibrate the geomechanical model yields accurate stress predictions with well-constrained uncertainties, providing a rigorous benchmark against which to evaluate models calibrated solely with indirect stress indicators.Our results demonstrate that stress predictions based solely on indirect observations achieve comparable accuracy to those calibrated with an exceptionally large and robust dataset of in-situ stress magnitude data. For the magnitude of the minimum horizontal stress Shmin, high-agreement scenarios reproduce the reference stress predictions throughout most of the stratigraphic section, with uncertainties dominated by natural rock property variability rather than stress magnitude uncertainty. For the magnitude of the maximum horizontal stress SHmax, the approach successfully delivers constrains within physically realistic ranges, though systematic overestimation of 2–3 MPa in some formations suggests remaining model limitations. This work demonstrates that indirect stress data, readily available during routine drilling operations, can provide reliable, uncertainty-quantified stress predictions without requiring expensive in-situ stress measurement campaigns, opening new possibilities for stress field characterization in subsurface projects worldwide.

Numerical simulations are widely used as tools to aid in understanding processes or making predictions about states and their evolution in time. However, in the process of simulation setup, a multitude of choices and simplifications must be made—from the definition of the physical laws implemented, to model discretization and spatial parameterization, and the definition of initial and boundary conditions. In addition, the defined parameters are often subject to significant uncertainties. It therefore becomes paramount to investigate the extent to which a model output, such as the predicted temperature or pressure state at a location, is influenced by a certain input parameter, such as the thermal conductivity of a geological unit. Attempts to resolve this question are thus included in many research papers, and are generally referred to as sensitivity analyses. However, only a small portion of studies actually perform a sensitivity analysis in a mathematical sense. We propose that this omission occurs for two reasons: (1) a lack of knowledge about the different types of sensitivity analyses, and how they can be properly described and formulated, and (2) a lack of awareness about the tools that are readily available to perform structured sensitivity analyses. In this contribution, we aim to fill this gap through a review of the mathematical foundations of sensitivity analyses, a presentation of the different types, several intuitive examples, and a recommendation for the choice of a suitable method, with the aim of providing guidance for wider application of suitable sensitivity analyses in subsurface process simulations. Graphic