Øystein Strengehagen Klemetsdal

Geothermal production and storage typically involve heat transport that is strongly localized around wells, in contrast to hydrocarbon recovery where pressure and multiphase flow processes act over much larger spatial scales. This spatial separation implies that significant computational savings can be achieved by concentrating resolution in regions with high temperature gradients while coarsening the grid elsewhere. In this work, we develop and investigate a dynamic coarsening framework for geothermal reservoir simulation in which coarse cells are constructed by aggregating an underlying fine-scale geomodel into a set of nested partitions. At each timestep, a dynamically adapted grid is assembled based on indicators that measure local solution complexity, enabling refinement near wells and other active regions while retaining coarser resolution elsewhere. The framework is flexible-applicable to general unstructured grids comprising arbitrary polyhedral cells-and computationally efficient because geometric and transmissibility properties for all partitions can be precomputed. Key novel elements include a sequential splitting scheme derived directly from the discrete conservation equations, a rigorous energy-conservative mapping between scales for arbitrary polyhedral grids, and geothermal-specific coarsening indicators that can be applied either statically, based on time-of-flight, or dynamically, updated at each timestep. We evaluate the method on representative underground thermal energy storage scenarios and compare results against fully resolved fine-scale simulations using various refinement indicators and coarsening strategies. All numerical experiments are performed using the geothermal module of the open-source MATLAB Reservoir Simulation Toolbox (MRST).

Summary: Newton methods commonly used for advancing implicit schemes in time often exhibit slow convergence, especially when initial guesses are far from the solution or when there are imbalances in nonlinearities across time and space. A common approach to deal with this is to relax the Newton updates, although this does not always yield good results. Alternatively, nonlinear solvers using local-global domain-decomposition preconditioning strategies have shown greater robustness. However, the observed reductions in iteration count do not always justify the higher computational cost in practical applications. In this work, we compare four different approaches to accelerating the nonlinear convergence of Newton methods. In the first and least invasive method, the convergence is continually monitored throughout the Newton loop to determine adjustments in update relaxation. The convergence monitors are computationally inexpensive and based on quantities that are usually available in commercial reservoir simulators. The second method is a straightforward nonlinear domain-decomposition (NLDD) preconditioning method that enhances the initial guess for Newton’s method by prepending each global Newton iteration with local subdomain solves. To reduce the computational overhead, local solves are adaptively bypassed in subdomains where there are small changes in mobility. The third method adds adaptive relaxation to the global step of the NLDD method. The fourth method uses convergence monitors to adaptively adjust relaxation in Newton and, if necessary, switches on local preconditioning solves depending if convergence issues are determined to be more severe. We investigate the effectiveness of the methods in practical models relevant to oil recovery, CO2 storage, and geothermal energy storage. Our results show that runtime can be significantly reduced in both single and ensemble models by enhancing the standard Newton’s method with dynamic relaxation. However, the most effective overall approach is NLDD preconditioning, where local solves are adaptively enabled or disabled based on mobility changes.

We discuss a nonlinear domain-decomposition preconditioning method for fully implicit simulations of multicomponent porous media flow based on the additive Schwarz preconditioned exact Newton method (ASPEN). The method efficiently accelerates nonlinear convergence by resolving unbalanced nonlinearities in a local stage and long-range interactions in a global stage. ASPEN can improve robustness and significantly reduce the number of global iterations compared with standard Newton, but extra work introduced in the local steps makes each global iteration more expensive. We discuss implementation aspects for the local and global stages. We show how the global-stage Jacobian can be transformed to the same form as the fully implicit system, so that one can use standard linear preconditioners and solvers. We compare the computational performance of ASPEN to standard Newton on a series of test cases, ranging from conceptual cases with simplified geometry or flow physics to cases representative of real assets. Our overall conclusion is that ASPEN is outperformed by Newton when this method works well and converges in a few iterations. On the other hand, ASPEN avoids time-step cuts and has significantly lower runtimes in time steps where Newton struggles. A good approach to computational speedup is therefore to adaptively switch between Newton and ASPEN throughout a simulation. A few examples of switching strategies are outlined.

Shallow geothermal reservoirs are excellent candidates for low-enthalpy energy storage, and can serve as heat batteries providing constant discharge of base heat, as well as rapid discharge of heat in periods of high demand. Recharging can be done by pumping down hot water, heated using excess heat from e.g. waste incinerators. In addition to having a very low carbon footprint, such systems also require limited surface infrastructure, and can easily be placed near or under the end-user, such as below residential buildings. The geological setting is typically complex, with horizons, faults, and intertwined patterns of natural fractures, and the nearwell region is often hydraulically fractured to enhance inter-well communication. Therefore, in order to fully utilize the potential of shallow geothermal heat storage, numerical simulations are imperative. In this work, we show how to practically model such systems, including generation of computational grids with a large number of wells and fractures, numerical discretizations with discrete fractures, and complex storage strategies with multiple wells working together under common group targets. We also discuss how adjoint-based methods can be used to tune model parameters (e.g. well injectivities, rock properties, and hydraulic fracture conductivities) so that the model fits observed data, and to find well controls (e.g. rates and temperatures) that optimize storage operations. The methodology is demonstrated using one artificial and two real underground thermal energy storage projects currently under development, and we highlight important challenges and our suggested solutions related to each of them. Thematic collection: This article is part of the Earth as a thermal battery: future directions in subsurface thermal energy storage systems collection available at: https://www.lyellcollection.org/topic/collections/thermal-energy

Domain decomposition methods are widely used as preconditioners for Krylov subspace linear solvers. In the simulation of porous media flow there has recently been a growing interest in nonlinear preconditioning methods for Newton’s method. In this work, we perform a numerical study of a spatial additive Schwarz preconditioned exact Newton (ASPEN) method as a nonlinear preconditioner for Newton’s method applied to both fully implicit or sequential implicit schemes for simulating immiscible and compositional multiphase flow. We first review the ASPEN method and discuss how the resulting linearized global equations can be recast so that one can use standard preconditioners developed for the underlying model equations. We observe that the local fully implicit or sequential implicit updates efficiently handle the local nonlinearities, whereas long-range interactions are resolved by the global ASPEN update. The combination of the two updates leads to a very competitive algorithm. We illustrate the behavior of the algorithm for conceptual one and two-dimensional cases, as well as realistic three dimensional models. A complexity analysis demonstrates that Newton’s method with a fully implicit scheme preconditioned by ASPEN is a very robust and scalable alternative to the well-established Newton’s method for fully implicit schemes.

Numerical smearing is oftentimes a challenge in reservoir simulation, particularly for complex tertiary recovery strategies. We present a new high-resolution method that uses dynamic coarsening of a fine underlying grid in combination with local timestepping to provide resolution in time and space. The method can be applied to stratigraphic and general unstructured grids, is efficient, introduces minimal computational overhead, and is applicable to flow models seen in practical reservoir engineering. Technically, the method is based on three concepts: We assess the applicability of the method through a set of representative cases, ranging from conceptual to realistic, with complex fluid physics and reservoir geology, and demonstrate that the method can be used to reduce computational time and still retain high resolution in spatial/temporal zones and quantities of interest.

The fully implicit method is the most commonly used approach to solve black-oil problems in reservoir simulation. The method requires repeated linearization of large nonlinear systems and produces ill-condi\-tioned linear systems. We present a strategy to reduce computational time that relies on two key ideas: (\textit{i}) a sequential formulation that decouples flow and transport into separate subproblems, and (\textit{ii}) a highly efficient Gauss--Seidel solver for the transport problems. This solver uses intercell fluxes to reorder the grid cells according to their upstream neighbors, and groups cells that are mutually dependent because of counter-current flow into local clusters. The cells and local clusters can then be solved in sequence, starting from the inflow and moving gradually downstream, since each new cell or local cluster will only depend on upstream neighbors that have already been computed. Altogether, this gives optimal localization and control of the nonlinear solution process. This method has been successfully applied to real-field problems using the standard first-order finite volume discretization. Here, we extend the idea to first-order dG methods on fully unstructured grids. We also demonstrate proof of concept for the reordering idea by applying it to the full simulation model of the Norne oil field, using a prototype variant of the open-source OPM Flow simulator.

A number of different multiscale methods have been developed as a robust alternative to upscaling and as a means for accelerated reservoir simulation of high-resolution geomodels. In their basic setup, multiscale methods use a restriction operator to construct a reduced system of flow equations on a coarser grid, and a prolongation operator to map pressure unknowns from the coarse grid back to the original simulation grid. The prolongation operator consists of basis functions computed numerically by solving localized flow problems. One can use the resulting multiscale solver both as a CPR-preconditioner in fully implicit simulators or as an efficient approximate iterative linear solver in a sequential setting. The latter approach has been successful implemented in a commercial simulator. Recently, we have shown that you can obtain significantly faster convergence if you instead of using a single pair of prolongation-restriction operators apply a sequence of such operators, where some of the operators adapt to faults, fractures, facies, or other geobodies. Herein, we present how you can accelerate the convergence even further, if you also include additional basis functions that capture local changes in the pressure.