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.

We apply Massively Parallel Interface for MPFA-O scheme with state-of-the-art Operator-Based Linearization (OBL) approach for multiphase flow in porous media. The implementation of MPFA-O scheme enhances the modelling capabilities for non-K-orthogonal mesh. A fully implicit scheme is applied to guarantee the stability of solutions when a mass-based formulation is involved to keep the flexibility of the framework for general-purpose reservoir simulation. As the MPFA-O introduces more non-zeros elements in the Jacobian matrix than the traditional TPFA, massively parallel computations via Message Passing Interface (MPI) in this work help to guarantee competitive computational efficiency for high-fidelity geological models. Concerning the Jacobian assembly hassle, we apply the OBL approach which introduces operators combining the fluid and rock properties in the conservation equations and discretizes the operators in the physical parameter space. By computing values and derivatives of the operators via a multilinear interpolation, the assembly of Jacobian matrix and residual vector could be drastically simplified. Another benefit of the OBL is that by only evaluating operator values on the predefined nodes in the physical parameter space, the overhead related to complex phase behavior and property evaluation is significantly reduced. In the end, we present several benchmark cases to rigorously demonstrate the accuracy, convergence, and robustness of the framework and two challenging field-scale cases to further prove its computing performance and parallel scalability.

The efficient operation and management of a geothermal project can be largely affected by geological, physical, operational and economic uncertainties. Systematic uncertainty quantification (UQ) involving these parameters helps to determine the probability of the focused outputs, e.g., energy production, Net Present Value (NPV), etc. However, how to efficiently assess the specific impacts of different uncertain parameters on the outputs of a geothermal project is still not clear. In this study, we performed a comprehensive UQ to a low-enthalpy geothermal reservoir using the GPU implementation of the Delft Advanced Research Terra Simulator (DARTS) framework with stochastic Monte Carlo samplings of uncertain parameters. With processing the simulation results, large uncertainties have been found in the production temperature, pressure drop, produced energy and NPV. It is also clear from the analysis that salinity influences the producing energy and NPV via changing the amount of energy carried in the fluid. Our work shows that the uncertainty in NPV is much larger than that in produced energy, as more uncertain factors were encompassed in NPV evaluation. An attempt to substitute original 3D models with upscaled 2D models in UQ demonstrates significant differences in the stochastic response of these two approaches in representation of realistic heterogeneity. The GPU version of DARTS significantly improved the simulation performance, which guarantees the full set (10,000 times) UQ with a large model (circa 3.2 million cells) finished within a day. With this study, the importance of UQ to geothermal field development is comprehensively addressed. This work provides a framework for assessing the impacts of uncertain parameters on the concerning system output of a geothermal project and will facilitate analyses with similar procedures.

Focus is on comparing stochastic, process-based and deterministic methods for modelling heterogeneity in hydraulic properties of fluvial geothermal reservoirs. Models are considered a generalized representation of a fluvial sequence in the upper part of the Gassum Formation in northern Denmark. Two ensemble realizations of process-based and stochastic heterogeneity were simulated using the state-of-the-art numerical modelling software, Delft Advanced Research Terra Simulator (DARTS), to assess differences on a statistically relevant sample size. Simulator settings were optimized to achieve two orders of magnitude improvement in simulation time. Our general findings show that the stochastic and process-based methods produce nearly identical results in terms of predicted breakthrough time and production temperature decline for high net-to-gross ratios (N/G). Simple homogenous and layered models overestimate breakthrough and underestimate temperature decline. More complex representation of facies in process-based models show smaller variance in results and stay within the limits of ensemble runs compared to simpler facies representation. Results indicate that a multivariate Gaussian based stochastic representation of heterogeneity provides comparable thermal response to a process-based model in fluvial systems of similar quality.

We present a scalable collocated Finite Volume Method (FVM) to simulate induced seismicity as a result of pore pressure changes. A discrete system is obtained based on a fully-implicit fully-coupled description of flow, elastic deformation, and contact mechanics at fault surfaces on a flexible unstructured mesh. The cell-centered collocated scheme leads to a convenient integration of the different physical equations, as the unknowns share the same discrete locations on the mesh. Additionally, a generic multi-point flux approximation is formulated to treat heterogeneity, anisotropy, and cross-derivative terms for both flow and mechanics equations. The resulting system, though flexible and accurate, can lead to excessive computational costs for field-relevant applications. To resolve this limitation, a scalable processing algorithm is developed and presented. Several proof-of-concept numerical tests, including benchmark studies with analytical solutions, are investigated. It is found that the presented method is indeed accurate and efficient; and provides a promising framework for accurate and efficient simulation of induced seismicity in various geoscientific applications.

We develop a collocated Finite Volume Method (FVM) to study induced seismicity as a result of pore pressure fluctuations. A discrete system is obtained based on a fully-implicit coupled description of flow, elastic deformation, and contact mechanics at fault surfaces on a fully unstructured mesh. The cell-centered collocated scheme leads to convenient integration of the different physical equations, as the unknowns share the same discrete locations on the mesh. Additionally, a multi-point flux approximation is formulated in a general procedure to treat heterogeneity, anisotropy, and cross-derivative terms for both flow and mechanics equations. The resulting system, though flexible and accurate, can lead to excessive computational costs for field-relevant applications. To resolve this limitation, a scalable parallel solution algorithm is developed and presented. Several proof-of-concept numerical tests, including benchmark studies with analytical solutions, are investigated. It is found that the presented method is indeed accurate, stable and efficient; and as such promising for accurate and efficient simulation of induced seismicity.

A realistic deep low-enthalpy geothermal reservoir based on real data with high detail and complicated sedimentary structure is utilized to perform sensitivity analyses of the geological features influencing reservoir properties. We perform simulations using the Delft Advanced Research Terra Simulator (DARTS). Compelling numerical performance of DARTS makes it suitable for handling a large ensemble of models including efficient sensitivity and uncertainty analyses. The major finding is that shale facies, generally ignored in hydrocarbon reservoir simulations, can significantly extend the predictive lifetime of geothermal reservoirs exploited by deep well doublets. It is important to accurately account for the shale facies in the simulation, though with an additional computational overhead. The overburden layers can improve doublet performance, but the impact depends on reservoir heterogeneity. In addition, heterogeneity will also divert the flow path with even a minor shift in the well placement. The discharge rate, an essential parameter of geothermal operation strategy, inversely corresponds to the doublet lifetime but positively correlates with the energy production for studied parameter ranges. Low sensitivity of doublet lifetime to vertical-horizontal permeability ratio and permeability-porosity correlation is observed. All these systematic findings for a realistic geothermal field with characterization at unprecedented level of detail can help to provide a general guideline for forward simulation and farther improve the profitability of geothermal energy production in realistic deep geothermal reservoirs through computer-assisted modeling and optimization.

A physics-based data-driven model is proposed for forecasting of subsurface energy production. The model fully relies on production data and does not require any in-depth knowledge of reservoir geology or governing physics. In the proposed approach, we use the Delft Advanced Reservoir Terra Simulator (DARTS) as a workhorse for data-driven simulation. DARTS uses an operator-based linearization technique that exploits an abstract interpretation of physics benefiting computational performance. The physics-based data-driven model is trained to fit data increasing the fidelity of the model forecast and reflecting significant changes in reservoir dynamics or physics over its history. The model is examined and validated for both synthetic and real field production data. We demonstrate that the developed approach is capable of providing accurate and reliable production forecast on a daily basis, even if the exact geological information is not available.

Alternative to CPU computing architectures, such as GPU, continue to evolve increasing the gap in peak memory bandwidth achievable on a conventional workstation or laptop. Such architectures are attractive for reservoir simulation, which performance is generally bounded by system memory bandwidth. However, to harvest the benefit of a new architecture, the source code has to be inevitably rewritten, sometimes almost completely. One of the biggest challenges here is to refactor the Jacobian assembly which typically involves large volumes of code and complex data processing. We demonstrate an effective and general way to simplify the linearization stage extracting complex physics-related computations from the main simulation loop and leaving only an algebraic multi-linear interpolation kernel instead. In this work, we provide the detailed description of simulation performance benefits from execution of the entire nonlinear loop on the GPU platform. We evaluate the computational performance of Delft Advanced Research Terra Simulator (DARTS) for various subsurface applications of practical interest on both CPU and GPU platforms, comparing particular workflow phases including Jacobian assembly and linear system solution with both stages of the Constraint Pressure Residual preconditioner.

Alternative to CPU computing architectures, such as GPU, continue to evolve increasing the gap in peak memory bandwidth achievable on a conventional workstation or laptop. Such architectures are attractive for reservoir simulation, which performance is generally bounded by system memory bandwidth. However, to harvest the benefit of a new architecture, the source code must be inevitably rewritten, sometimes almost completely. One of the biggest challenges here is to refactor the Jacobian assembly which typically involves large volumes of code and complex data processing. We demonstrate an effective and general way to simplify the linearization stage extracting complex physics-related computations from the main simulation loop and leaving only an algebraic multi-linear interpolation kernel instead. In this work, we provide the detailed description of simulation performance benefits from execution of the entire nonlinear loop on the GPU platform. We evaluate the computational performance of Delft Advanced Research Terra Simulator (DARTS) for various energy transition subsurface applications of practical interest on both CPU and GPU platforms, comparing particular workflow phases including Jacobian assembly and linear system solution with both stages of the Constraint Pressure Residual preconditioner.

A physics-based data-driven model is proposed in this study for the forecasting of secondary oil recovery. The model fully relies on production data and does not directly requires any in-depth knowledge of the reservoir geology or governing physics. In the proposed approach, we utilise Delft Advanced Reservoir Terra Simulator (DARTS) as a base for data-driven simulation. DARTS uses an Operator-Based Linearization technique which exploits an abstract interpretation of physics benefiting computational performance for a forward simulation. The proposed strategy was evaluated first on the two synthetic data ensembles and showed good prediction accuracy for a significantly reduced model size. Besides, the data-driven proxy methodology was compared with an advanced flow-based upscaling technique and demonstrated an improved accuracy for both ensembles. Besides, the proposed data-driven approach was examined on two realistic data sets. For the first case, the methodology demonstrates advanced predictive performance for training based on synthetic data generated from a high-fidelity simulation model with imposed random noise. To check the robustness of the proposed methodology, the control parameters for a forecast period were significantly changed in comparison to the training period. The data-driven model still manages to predict the forecast production quite close to the reference high-fidelity results. However, the training performed on another data set based on historical production from a real brownfield was not fully successful. We relate a bigger error in both training and forecast period for this model to poor data quality. The training procedure for this model led to a moderate accuracy in history matching for a long production period, where general production trends have resembled true data and water breakthrough time was restored in nearly all wells. However, there are still periods of poor accuracy, especially where shark peaks and falls are experienced.

The continuous progress of reservoir monitoring technology provides encouraging capacities to reduceuncertainties in the subsurface characterization and to mitigate risks in field development applying thereservoir simulation approach. However, it is always challenging to take full advantage of the observationdata, since an accurate representation of strong heterogeneities requires a high-resolution grid. Most ofthe discretization methods cannot handle full tensor permeability, and high nonlinearity introduced bycomplex physical process drastically reduces simulation efficiency. In this work, we develop an advancedparallel framework for reservoir simulation with the implementation of state of the art discretizationand linearization methods. We apply the multipoint flux approximation (MPFA) method to handle thefull tensor permeability in unstructured grids. To keep the fidelity of the geological model and improvecomputational efficiency, we use massively parallel computations via Message Passing Interface (MPI).Complex subsurface physics is described by mass-based formulations making the framework flexible forgeneral-purpose reservoir simulation. However, the representation of phase behavior introduces additionalworkload when compared with the phase-based formulations in the traditional approach. Here, we apply theOperator-Based Linearization (OBL) approach which not only overcomes this drawback but also turns it toan advantage. In this method, the conservation equations are described in an operator form. By constructinga library of tabulated operators, the repeated work spent on complex phase behavior and property evaluationcan be significantly reduced. We benchmark the parallel framework with analytical solutions under single-phase flow and multiphase flow. The results demonstrate that the parallel framework provides accuratesimulation results for structured and unstructured grids. We validate that MPFA implemented in our parallelframework converges to real solutions when the permeability is a full tensor. Besides, several realisticcases have been rigorously tested confirming high computational capacity, efficiency, and accuracy of theadvanced massively parallel framework for general-purpose reservoir simulation. With the implementationof MPFA and OBL approaches, the parallel framework is fully equipped for the simulation of problemswith full tensor permeability, high-heterogeneities, and complex physical processes.

The efficient development of a geothermal field can be largely affected by the inherent geological and physical uncertainties. Besides, the uncertain operational and economic parameters can also impact the profit of a project. Systematic uncertainty quantification involving these parameters helps to determine the probability of concerning outputs. In this study, a low-enthalpy geothermal reservoir with strong heterogeneity, located in the West Netherlands Basin, is selected as the research area. Detailed geological model is constructed based on various static data including seismic and log interpretation. However, significant uncertainties still exist in definition of the model parameters, mainly reservoir permeability and porosity. Besides, the fluid properties have not been sampled in this field and can vary in the range between brackish to highly saline water. Also, the heat price and operational investment fluctuate with time and add up to uncertainty. Taking all interested parameters into consideration, the Monte Carlo method is utilized to select specific input data set. The forward simulations are powered by the GPU version of Delft Advance Research Terra Simulator (DARTS), which provides efficient simulation capabilities for geothermal applications. Through this investigation, a wide range of production temperature has been observed due to the uncertainty of the input parameters.

Numerical modelling of multiphase multicomponent flow coupled with mass and energy transport in porous media is crucially important for many applications including oil recovery, carbon storage and geothermal. To deliver robust simulation results, a fully or adaptive implicit method is usually employed, creating a highly nonlinear system of equations. It is then solved with the Newton-Raphson method, which requires a linearization procedure to assemble a Jacobian matrix. Operator Based Linearization (OBL) approach allows detaching property computations from the linearization stage by using piece-wise multilinear approximations of state-dependent operators related to complex physics. The values of operators used for interpolation are computed adaptively in the parameter-space domain, which is uniformly discretized with the desired accuracy. As the result, the simulation performance does not depend on the cost of property computations, making it possible to use expensive equation-of-state formulations (e.g., fugacity-activity thermodynamic models) or even black-box chemical packages (e.g., PHREEQC) for an accurate representation of governing physics without penalizing runtime. On the other hand, the implementation of the simulation framework is significantly simplified, which allows improving the simulation performance further by executing the complete simulation loop on GPU architecture. The integrated open-source framework Delft Advanced Research Terra Simulator (DARTS) is built around the OBL concept and provides a flexible, modular and computationally efficient modelling package. In this work, we evaluate the computational performance of DARTS for various subsurface applications of practical interests on both CPU and GPU platforms. We provide a detailed performance comparison of particular workflow pieces composing Jacobian assembly and linear system solution, including both stages of Constrained Pressure Residual solver.

Accurate prediction of temperature and pressure distribution is essential for geothermal reservoir exploitation with cold water re-injection. Depending on our knowledge about the heterogeneous structure of the subsurface, the reservoir development scheme can be optimized and the overall lifetime of the geothermal field can be extended. In this study, we present Delft Advanced Research Terra Simulator (DARTS), which provides fast and accurate energy production evaluation for geothermal applications. This simulation framework is suitable for uncertainty analysis with a large ensemble of models. In DARTS, we select the molar formulation with pressure and enthalpy as primary variables. Besides, the fully-coupled fully-implicit two-point flux approximation on unstructured grids is implemented to solve the mass and energy conservation equations. For the nonlinear solution, we employ the recently developed Operator-Based Linearization (OBL) approach. In our work, DARTS is compared with the state-of-the-art simulation frameworks using a set of benchmark tests. We demonstrate that DARTS achieves a good match for both low- and high-enthalpy conditions in comparison to other simulators. At the same time, DARTS provides high performance and flexibility of the code due to the OBL approach, which makes it particularly useful for uncertainty quantification in processes involving complex physics.

Various novel computing architectures, like massively parallel and multi-core, as well as computing accelerators, like GPUs or TPUs, keep regularly expanding. In order to exploit the benefits of these architectures to the full extent and speed up reservoir simulation, the source code has to be inevitably rewritten, sometimes almost completely. We demonstrate how to extract complex physics-related computations from the main simulation loop, leaving only an algebraic multilinear interpolation kernel instead. In combination with linear solvers, which usually have made available soon once the new architecture is introduced, the approach accommodates execution of the entire nonlinear loop on the latest hardware and computational architectures. We describe the integrated simulation framework built on top of this technique and show the applicability of the approach to various challenging physical and chemical problems. All simulation engines along with linear solvers, well controls, interpolation engines, and state operator evaluators are implemented in C++11 and exposed into Python coupling the flexibility of the script language with the performance of C++.

Numerical simulation is one of the most important tools required for financial and operational management of geothermal reservoirs. The modern geothermal industry is challenged to run large ensembles of numerical models for uncertainty analysis, causing simulation performance to become a critical issue. Geothermal reservoir modeling requires the solution of governing equations describing the conservation of mass and energy. The robust, accurate and computationally efficient implementation of this solution suggests an implicit time-approximation scheme, which introduces nonlinearity into the system of equations to be solved. The most commonly used approach to solving the system of nonlinear equations is based on Newton's method and involves linearization with respect to nonlinear unknowns. This stage is the most complicated for implementation and usually becomes the source of various errors. A new linearization approach – operator-based linearization – was recently proposed for non-isothermal flow and transport. The governing equations, discretized in space and time, were transformed to the operator form where each term of the equation was specified as the product of two operators. The first operator comprises physical properties of rock and fluids, such as density or viscosity, which depend only on the current state of a grid block, fully defined by the values of nonlinear unknowns. The second operator includes all terms that were not included in the first operators, and depends on both the state and spatial position of a control volume. Next, the first type of operators was parametrized over the physical space of a simulation problem. The representation of highly nonlinear physics was achieved by using multi-linear interpolation, which replaces the continuous representation of parametrized operators. The linearization of the second type of operators was applied in the conventional manner. In this work, we investigated the applicability of this approach to the geothermal processes, specifically for low-enthalpy and high-enthalpy geothermal doublet models with hydrocarbon co-production. The performance and robustness of the new method were tested against the conventional approach on a geothermal reservoir of practical interest. This approach shows significant improvement of geothermal simulation performance, while errors, introduced by coarsening in physics, remain under control. The simplicity of implementation on emerging computational architectures and nonlinearity reduction provide advanced opportunities for uncertainty quantification and risk analysis of geothermal projects.

The nonlinear nature of flow and transport in porous media requires a linearization of the governing numerical-model equations. We propose a new linearization approach and apply it to complex thermal/compositional problems. The key idea of the approach is the transformation of discretized mass- and energy-conservation equations to an operator form with separate space-dependent and statedependent components. The state-dependent operators are parameterized using a uniformly distributed mesh in parameter space. Multilinear interpolation is used during simulation for a continuous reconstruction of state-dependent operators that are used in the assembly of the Jacobian and residual of the nonlinear problem. This approach approximates exact physics of a simulation problem, which is similar to an approximate representation of space and time discretization performed in conventional simulation. Maintaining control of the error in approximate physics, we perform an adaptive parameterization to improve the performance and flexibility of the method. In addition, we extend the method to compositional problems with buoyancy. We demonstrate the robustness and convergence of the approach using problems of practical interest.