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.@en

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.@en

Numerical simulation is based on space and time discretization providing a trade-off between accuracy and computational performance. The operator-based Linearization (OBL) approach introduces an additional discretization domain for the physical description of fluid and rock. Delft Advanced Research Terra Simulator (DARTS) is a general-purpose reservoir simulation platform entirely built around the OBL approach. DARTS provides unique flexibility capabilities, allowing to customize physical description and control simulation process using the high-level Python programming language. Fully Implicit Method provides unconditional simulation stability, while the partial derivatives required for Jacobian assembly are computed analytically through OBL. At the same time, DARTS ensures exceptional simulation performance addressing it at three levels. Inexpensive linearization combined with a reduction in the nonlinearity is controlled by OBL discretization resolution form algorithmic level. Efficient C++ backend for critical kernels and state-of-the-art linear solver with two-stage CPR preconditioner composes the software level. Finally, the ability to perform the entire simulation loop on the GPU platform owing to the OBL approach constitutes the hardware level. The DARTS framework has already served as a platform for several academic and industrial research projects for different geo-energy applications including geothermal, CO2 sequestration and petroleum.@en

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.@en

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.@en

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.@en

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.@en

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.@en

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.@en

Simulation of geothermal processes is based on the solution of strongly nonlinear governing equations describing flow of mass and energy in the subsurface. The solution of this problem requires a linearization of governing equations. Recently, a new approach based on the operator-based multi-linear representation of conservation equations was proposed for an isothermal flow and transport. The discretized version of conservation equations was transformed to the operator form where each term of the equation was represented as a product of two operators: one dependent on physical properties of fluids (e.g. density or viscosity) and another dependent on properties, altered in space (e.g. porosity or permeability). The first type of operators was parametrized in the physical space of the simulation problem at the pre-processing stage, or adaptively during the course of simulation. Next, multi-linear interpolation was applied for physics-based operators when the linearization of the second type of operators was performed in a conventional manner. In this work, the approach has been extended for the thermal processes. The operator-based parametrization was performed using temperature, pressure and composition. The performance and robustness of the method was tested against the conventional approach on low enthalpy models of geothermal reservoirs of practical interest. @en

General purpose reservoir simulation is based on the solution of governing equations describing mass and energy transfer in the subsurface. The solution process requires the linearization of strongly nonlinear governing equations. Usually, a Newton-based method is used for the linearization. This method demands the assembly of Jacobian and residual for a fully coupled system of equations. Recently, a new linearization approach was proposed and tested for binary systems. The key idea of the Operator Based Linearization (OBL) approach is to transform the discretized mass and energy conservation equations to an operator form which separates space-dependent and state-dependent properties of governing equations. This transformation provides the opportunity to approximate the representation of exact physics (physical properties) of a problem. Specifically, each term of the conservation equations is presented as the product of two different operators. The first operator depends on the current physical state of a system and contains fluid properties, such as density, viscosity, relative permeability, etc. The second operator captures both spatially altered properties, such as permeability, and the rest of state variables, such as pressure in the discrete approximation of gradient. All state-dependent operators are uniformly parametrized within the physical space of the problem (pressure-composition intervals). During simulation process, a multi-linear interpolation is applied to approximate the first type of operators, while the second type of operators is processed based on conventional approach. In this work, we extended the approach to thermal systems with an arbitrary number of components. Besides, we significantly improved the performance of OBL employing adaptive parametrization technique. We tested the approach for truly multi-component thermal systems of practical interest. The computational performance, accuracy, and robustness of a new method were demonstrated against the conventional approach.@en

Numerical simulation of coupled multiphase multicomponent flow and transport in porous media is a crucial tool for understanding and forecasting of complex industrial applications related to the subsurface. The discretized governing equations are highly nonlinear and usually need to be solved with Newton's method, which corresponds with high computational cost and complexity. With the presence of capillary and gravity forces, the nonlinearity of the problem is amplified even further, which usually leads to a higher numerical cost. A recently proposed operator-based linearization (OBL) approach effectively improves the performance of complex physical modeling by transforming the discretized nonlinear conservation equations into a quasilinear form according to state-dependent operators. These operators are approximated by means of a discrete representation on a uniform mesh in physical parameter space. Continuous representation is achieved through the multilinear interpolation. This approach provides a unique framework for the multifidelity representation of physics in general-purpose simulation. The applicability of the OBL approach was demonstrated for various energy subsurface applications with multiphase flow of mass and heat in the presence of buoyancy and diffusive forces. In this work, the OBL approach is extended for multiphase multicomponent systems with capillarity. Through the comparisons with a legacy commercial simulator using a set of benchmark tests, we demonstrate that the extended OBL scheme significantly improves the computational efficiency with the controlled accuracy of approximation and converges to the results of the conventional continuous approach with an increased resolution of parametrization.@en

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.@en

The current situation with green gas emission requires the development of low carbon energy solutions. However, a significant part of the modern energy industry still relies on fossil fuels. To combine these two contradictory targets, we investigate a strategy based on a combination of CO2 sequestration with Enhanced Oil Recovery (EOR) in the hydrocarbon reservoirs. In such technology, the development of miscibility is the most attractive strategy from both technological and economic aspects. Modeling of this process involves solving complex nonlinear problem describing compositional flow and transport in highly heterogeneous porous media. An accurate capture of the miscibility development usually requires an extensive number of components to be present in the compositional problem which makes simulation run-time prohibitive for optimization. Here, we apply a multi-scale reconstructing of compositional transport to the optimization of CO2 injection. In this approach, a prolongation operator, based on the parametrization of injection and production tie-lines, is constructed following the fractional flow theory. This operator is tabulated as a function of pressure and pseudocomposition which then is used in the Operator-Based Linearization (OBL) framework for simulation. As a result, a pseudo two-component solution of the multidimensional problem will match the position of trailing and leading shocks of the original problem which helps to accurately predict phase distribution. The reconstructed multicomponent solution can be used then as an effective proxy-model mimicking the behavior of the original multicomponent system. Next, we use this proxy-model in the optimization procedure which helps to improve the performance of the process in several folds. An additional benefit of the proposed methodology is based on the fact that important technological features of CO2 injection process can be captured with lower degrees of freedom which makes the optimization solution more feasible.@en

As oil production continues worldwide, more oil fields require complex EOR methods to achieve outlined recovery factors. Reservoir engineers are dealing more often with problems involving thermal multiphase multi-component flow models tightly coupled with complex phase behavior. Such modeling implies the solution of governing laws describing mass and energy transfer in the subsurface, which in turn requires the linearization of strongly nonlinear systems of equations. The recently proposed Operator-Based Linearization (OBL) framework suggests an unconventional strategy using the discrete representation of physics. The terms of governing PDEs, discretized in time and space, which depend only on state variables, are approximated by piece-wise multilinear operators. Since the current physical state fully defines operators for a given problem, each operator can be parametrized over the multidimensional space of nonlinear unknowns for a given distribution of supporting points. Onwards, the values of operators, along with their derivatives with respect to nonlinear unknowns, are obtained from the parametrization using multilinear interpolation and are used for Jacobian assembly in the course of a simulation. Previously, the distribution of supporting points was always uniform, requiring higher parametrization resolution to provide accurate and consistent interpolation of an operator around its most nonlinear regions in parameter space. In addition, when the resolution is low, the system can lose hyperbolicity causing convergence issues. In this work, we apply the prior knowledge of underlying physics to distribute the supporting points according to the tie-simplex behavior of the multiphase mixture in parameter space. The approach allows to decrease the parametrization resolution keeping the same accuracy. In addition, the OBL framework is extended to describe multisegment wells working under different controls. We test the accuracy of the developed framework for truly multi-component systems of practical interest.@en

Non-isothermal multiphase compositional simulation is based on the solution of governing equations describing mass and energy transfer in the subsurface. The solution strategy requires a linearization of strongly nonlinear governing equations describing the process. Usually, a Newton-based method is used for the linearization that demands an assembly of a Jacobian matrix and residuals for a fully coupled system of equations. Recently, a new linearization approach was proposed for compositional problems and tested for simulation of binary compositional and low-enthalpy geothermal flow. The key idea of the approach is the transformation of discretised mass conservation equations to an operator form with separate space-dependent and state-dependent components. This transformation provides an opportunity for an approximate representation of exact physics (physical properties) of the problem. Specifically, each term of conservation equations is represented as a product of two different operators. The first operator depends on a current physical state of a system and contains different properties such as density, viscosity, relative permeability, etc. The second operator captures both spatially altered properties such as permeability and the rest of state variables such as pressure in the discrete approximation of the gradient. At the pre-processing stage, all state-dependent operators are uniformly parametrized within the physical space of the problem (pressure-composition intervals). During the simulation process, a multi-linear interpolation is applied to approximate the first type of operators, while the second type of operators is processed based on the conventional approach. In this work, we have extended this approach to general purpose simulation. We introduced the operator-based parametrization of mass and energy conservations equation based on the pressure, composition, temperature, and porosity. In addition, the approach has been extended and tested on truly multi-component systems of practical interest. The accuracy and robustness of the new method have been tested against the results of simulations based on the conventional approach.@en

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.@en

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++.@en

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.@en