The open Delft Advanced Research Terra Simulator (open-DARTS) framework is an open-source reservoir simulation software. The open-DARTS focused on energy transition applications, such as geothermal energy production and carbon sequestration. It enables the modeling of compositional thermal flow, coupled with a geomechanical solver based on the Finite Volume discretization and adjoints method for inverse modeling. The open-DARTS supports different grid types (structured, corner-point geometry, and unstructured), discrete fracture networks, contact mechanics, and various thermal-chemical interactions. The recently proposed generic nonlinear formulation supports the most general nonlinear PDEs designed for various energy transition applications. The open-DARTS has been implemented in C++ and Python to optimize hardware utilization while ensuring flexibility. The most computationally expensive part is written in C++ and compiled into libraries, which are subsequently exposed to Python using pybind11. This allows the extension and overriding of C++ functions by user-defined Python code. For example, using only a Python interface, one can adjust a timestep strategy, nonlinear solver, or properties output. Besides, the Python interface of open-DARTS provides straightforward coupling with other Python-based numerical modeling packages, including the meshing, file storage, caching, and visualization modules. The open-DARTS core uses the advantages of C++ language, such as efficient low-level memory management, object-oriented programming, compile-time polymorphism, and parallelization with OpenMP. One of the advantages of open-DARTS is the Operator-Based Linearization (OBL) technique, which can resolve challenges associated with complex physics and reduce the computation time, especially for ensemble-based simulations. We would also like to share our experience on the project, repository, and the development workflow configuration using gitlab.com, including the build system (cmake), handling merge requests, automated testing in CI/CD pipelines, documentation management (gitlab.io), wiki utilization, and release publishing. Additionally, Python’s integration into open-DARTS offers the advantage of straightforward installation via PyPI and simplifies defining requirements for users who prefer to avoid compiling code from source files.

This paper presents and discusses the results obtained by the participants to the benchmark described in de Hoop et al, Comput. Geosci. (2024). The benchmark uses a model for CO₂ geological storage and focuses on the coupling between two-phase flow and geochemistry. Several test cases of various levels of difficulty are proposed, both in one and two spatial dimensions. Six teams participated in the benchmark, each with their own simulation code, though not all teams attempted all the cases. The codes used by the participants are described, and the results obtained on the various test cases are compared, as well as the performance of the codes. It is shown that the results obtained are widely consistent, giving a good level of confidence in the outcome of the benchmark. The general complexity of two-phase flow coupled with chemical reactions altering porous media means that some differences between the codes remain. Besides, from the convergence study, it is clear that the two-dimensional problem has a relatively high sensitivity to a spatial resolution which adds to the complexity.

Fracture networks are abundant in subsurface applications (e.g., geothermal energy production, CO₂ sequestration). Fractured reservoirs often have a very complex structure, making modeling flow and transport in such networks slow and unstable. Consequently, this limits our ability to perform uncertainty quantification and increases development costs and environmental risks. This study provides an advanced methodology for simulation based on Discrete Fracture Model approach. The preprocessing framework results in a fully conformal, uniformly distributed grid for realistic 2D fracture networks at a required level of precision. The simplified geometry and topology of the resulting network are compared with input (i.e., unchanged) data to evaluate the preprocessing influence. The resulting mesh-related parameters, such as volume distributions and orthogonality of control volume connections, are analyzed. Furthermore, changes in fluid-flow response related to preprocessing are evaluated using a high-enthalpy two-phase flow geothermal simulator. The simplified topology directly improves meshing results and, consequently, the accuracy and efficiency of numerical simulation. The main novelty of this work is the introduction of an automatic preprocessing framework allowing us to simplify the fracture network down to required level of complexity and addition of a fracture aperture correction capable of handling heterogeneous aperture distributions, low connectivity fracture networks, and sealing fractures. The graph-based framework is fully open-source and explicitly resolves small-angle intersections within the fracture network. A rigorous analysis of changes in the static and dynamic impact of the preprocessing algorithm demonstrates that explicit fracture representation can be computationally efficient, enabling their use in large-scale uncertainty quantification studies.

The positive impact that natural fractures can have on geothermal heat production from low-permeability reservoirs has become increasingly recognised and proven by subsurface case studies. In this study, we assess the potential impact of natural fractures on heat extraction from the tight Lower Buntsandstein Subgroup targeted by the recently drilled NLW-GT-01 well (West Netherlands Basin (WNB)). We integrate: (1) reservoir property characterisation using petrophysical analysis and geostatistical inversion, (2) image-log and core interpretation, (3) large-scale seismic fault extraction and characterisation, (4) Discrete Fracture Network (DFN) modelling and permeability upscaling, and (5) fluid-flow and temperature modelling. First, the results of the petrophysical analysis and geostatistical inversion indicate that the Volpriehausen has almost no intrinsic porosity or permeability in the rock volume surrounding the NLW-GT-01 well. The Detfurth and Hardegsen sandstones show better reservoir properties. Second, the image-log interpretation shows predominately NW-SE-orientated fractures, which are hydraulically conductive and show log-normal and negative-power-law behaviour for their length and aperture, respectively. Third, the faults extracted from the seismic data have four different orientations: NW-SE, N-S, NE-SW and E-W, with faults in proximity to the NLW-GT-01 having a similar strike to the observed fractures. Fourth, inspection of the reservoir-scale 2D DFNs, upscaled permeability models and fluid-flow/temperature simulations indicates that these potentially open natural fractures significantly enhance the effective permeability and heat production of the normally tight reservoir volume. However, our modelling results also show that when the natural fractures are closed, production values are negligible. Furthermore, because active well tests were not performed prior to the abandonment of the Triassic formations targeted by the NLW-GT-01, no conclusive data exist on whether the observed natural fractures are connected and hydraulically conductive under subsurface conditions. Therefore, based on the presented findings and remaining uncertainties, we propose that measures which can test the potential of fracture-enhanced permeability under subsurface conditions should become standard procedure in projects targeting deep and potentially fractured geothermal reservoirs.

The main objective of this study is to perform Uncertainty Quantification (UQ) using a detailed representation of fractured reservoirs. This is achieved by creating a simplified representation of the fracture network while preserving the main characteristics of the high-fidelity model. We include information at different scales in the UQ workflow which allows for a large reduction in the computational time while converging to the high-fidelity full ensemble statistics. Ultimately, it allows us to make accurate predictions on geothermal energy production in highly heterogeneous fractured porous media. The numerical reservoir simulation tool we use in this work is the Delft Advanced Research Terra Simulator (DARTS). It is based on Finite Volume approximation in space, fully coupled explicit approximation in time, and uses the novel linearization technique called Operator-Based Linearization (OBL) for the system of discretized nonlinear governing equations. We use a fracture network generation algorithm that uses distributions for length, angles, size of fracture sets, and connectivity as its main input. This allows us to generate a large number of complex fracture networks at the reservoir scale. We developed a pre-processing algorithm to simplify the fracture network and use graph theory to analyze the connectivity before and after pre-processing. Furthermore, we use metric space modeling methods for statistical analysis. A robust coarsening method for the Discrete Fracture Matrix model (DFM) is developed which allows for great control over the degree of simplification of the network’s topology and connectivity. We apply the framework to modeling of geothermal energy extraction. The pre-processing algorithm allows for a hierarchical representation of the fracture network, which in turn is utilized in the reduced UQ methodology. The reduced UQ workflow uses the coarser representation of the fracture networks to partition/rank the high-fidelity parameter space. Then a small subset of high-fidelity models is chosen to represent the full ensemble statistics. Hereby, the computational time of the UQ is reduced by two/three orders of magnitude, while converging to similar statistics as the high-fidelity full ensemble statistics. The methods developed in this study are part of a larger project on a prediction of energy production from carboniferous carbonates. The final goal is to perform UQ in pre-salt reservoirs which are characterized by complex reservoir architecture (i.e., large karstification and fracture networks). The UQ of fractured reservoirs is typically done in the petroleum industry using effective media models. We present a methodology that can efficiently handle a large ensemble of DFM models, which represent complex fracture networks and allow for making decisions under uncertainty using more detailed high-resolution numerical models.

Multiphase mass and heat transfer are ubiquitous in the subsurface within manifold applications. The presence of fractures over several scales and complex geometry magnifies the uncertainty of the heat transfer phenomena, which will significantly impact, or even dominate, the dynamic transport process. Capturing the details of fluid and heat transport within the fractured system is beneficial to the subsurface operations. However, accurate modeling methodologies for thermal high-enthalpy multiphase flow within fractured reservoirs are quite limited. In this work, multiphase flow in fractured geothermal reservoirs is numerically investigated. A discrete-fracture model is utilized to describe the fractured system. To characterize the thermal transport process accurately and efficiently, the resolution of discretization is necessarily optimized. A synthetic fracture model is firstly selected to run on different levels of discretization with different initial thermodynamic conditions. A comprehensive analysis is conducted to compare the convergence and computational efficiency of simulations. The numerical scheme is implemented within the Delft Advanced Research Terra Simulator (DARTS), which can provide fast and robust simulation to energy applications in the subsurface. Based on the converged numerical solutions, a thermal Péclet number is defined to characterize the interplay between thermal convection and conduction, which are the two governing mechanisms in geothermal development. Different heat transfer stages are recognized on the Péclet curve in conjunction with production regimes of the synthetic fractured reservoir. A fracture network, sketched and scaled up from a digital map of a realistic outcrop, is then utilized to perform a sensitivity analysis of the key parameters influencing the heat and mass transfer. Thermal propagation and Péclet number are found to be sensitive to flow rate and thermal parameters (e.g., rock heat conductivity and heat capacity). This paper presents a numerical simulation framework for fractured geothermal reservoirs, which provides the necessary procedures for practical investigations regarding geothermal developments with uncertainties.

A coupled description of flow and thermal-reactive transport is spanning a wide range of scales in space and time, which often introduces a significant complexity for the modelling of such processes. Subsurface reservoir heterogeneity with complex multi-scale features increases the modelling complexity even further. Traditional multiscale techniques are usually focused on the accuracy of the pressure solution and often ignore the transport. Improving the transport solution can however be quite significant for the performance of the simulation, especially in complex applications related to thermal-compositional flow. The use of an Adaptive Mesh Refinement enables the grid to adapt dynamically during the simulation, which facilitates the efficient use of computational resources. This is especially important in applications with thermal flow and transport where the region requires high-resolution calculations as often localized in space. In this work, the aim is to develop an Adaptive Mesh Refinement framework for geothermal reservoir simulation. The approach uses a multi-level connection list and can be applied to fully unstructured grids. The adaptivity of the grid in the developed framework is based on a hierarchical connectivity list. First, the fine-scale model is constructed, which accurately approximates all reservoir heterogeneity. Next, a global flow-based upscaling is applied, where an unstructured partitioning of the original grid is created. Once the full hierarchy of levels is constructed, the simulation is started at the coarsest grid. Grid space refinement criteria is based on the local changes and can be adjusted for specific models and governing physics. The multi-level connectivity lists are redefined at each timestep and used as an input for the next. The developed Adaptive Mesh Refinement framework was implemented in Delft Advanced Research Terra Simulator which uses the Operator-Based Linearization technique. The performance of the proposed approach is illustrated for several challenging geothermal applications of practical interest.

Caves developed in carbonate units have a significant role in fluid flow, but most of these subsurface voids are below seismic resolution. We concentrated our study on four caves to determine the roles of fractures and folds in the development of karst conduits that may form flow pathways in carbonate reservoirs. We performed structural field investigations, petrographic analyses, and geometric characterization using Light Detection and Ranging (LIDAR) for caves in Neoproterozoic carbonates of the Salitre Formation, central part of the São Francisco Craton, Brazil. We found that the conduit shape, usually with an ellipsoidal cross-section, reflects the tectonic features and textural variations. Carbonate layers containing pyrite and low detritic mineral contents are generally karstified and appear to act as favorable flow pathways. Our results indicate that the development of the karst system is related to fracture corridors formed along parallel and orthogonal sets of fold hinges, which provide preferential pathways for fluid flow and contribute to the development of super-K zones. This study provides insights into the prediction of subseismic-scale voids in carbonate reservoirs, with direct application for the hydrocarbon and hydrogeology flow and storage.

The hypogenic caves developed in carbonate units have a significant structural control but most of their features are not detect by conventional methods due to their size below seismic resolution. This contribution focuses on the structural, petrographic and geometric characterization of karst conduits in Neoproterozoic carbonates of the Salitre Formation, central part of São Francisco Craton, Brazil. We address the influence of fractures and folds on the development of karst conduits through field and laboratory analysis and the application of Light Detection Ranging to characterize cave/conduit geometry. The preliminary results indicate that the process of karstification are intensified in fractures corridors developed along fold hinges, which create fluid flow corridors in carbonate units and may change petrophysical reservoir properties.

The Morro Vermelho Cave (MVC) (Brazil) developed within the Morro Vermelho karst system, which affected Neoproterozoic limestones (Salitre Formation). The MVC experienced little interactions with meteoric processes and is an example of a hypogenic cave formed during strike-slip deformation. The Salitre carbonates in the MVC experienced distributed deformation along an elongated domain overlying a buried strike-slip fault. Gently dipping, semiductile shear zones formed with decimeter-scale (3.9 in.) dolomitic veins. In our model, Mg-rich fluids flowing along the Salitre aquifer caused at the same time extensive dolomitization of the body of rock (100-m [328-ft] scale) experiencing distributed deformation. With progressive displacement, the deep strike-slip fault propagated upward causing the development of an anticline pop-up, steepening sedimentary layers, and steep 1-10-m-long (3.3-33.8-ft) fractures, which served as pathways for upward fluid flow. These steep extensional fractures made it possible for fluids flowing in lower, quartzitic aquifers to enter the carbonate aquifer causing silica deposition in rock cavities and in fractures and fault planes. Following the main stage of speleogenesis, silica deposition took over again depositing on the cave walls a continuous silica crust, rarely observed in other settings worldwide. The interplay between regional bedding-parallel flow and focused circulation of fluids along steep faults and dipping layers, and the associated rock-fluid interactions are not unique to the contractional settings presented but can also occur in association with similar faults in rifted continental margins.

In this study, an attempt is made to better understand the effect coarsening of the parameter space has on the uncertainty representation of the response. Firstly, an HF ensemble of channelized reservoir models is constructed using a Multi-Point Statistic (MPS) approach. Several levels of coarsening are generated using a flow-based upscaling algorithm. A water injection strategy is simulated for each scale of the hierarchical ensemble. Dynamic analysis is performed on a reduced representation of the response uncertainty obtained via Multidimensional Scaling (MDS). We introduce an Uncertainty Trajectory (UT), which quantifies the coarsening effect in terms of deviation from the HF ensemble response uncertainty. The UT also includes the temporal behavior of the response uncertainty of each ensemble scale. The mean integrated distance from the HF ensemble UT can be used as a measure of dissimilarity in the flow behavior of consecutive coarser ensembles scales. Reducing the number of HF flow simulations required for uncertainty quantification can be achieved via the proposed methodology and thereby greatly reducing the overall computational cost.