Olav Møyner

Summary: The 11th SPE Comparative Solution Project (SPE11) was established to benchmark the ability of reservoir simulators to capture the key physical processes of geological carbon dioxide (CO2) storage. A central question in this project was whether grid refinement could yield converged solutions for quantities of practical interest. In this work, we revisit Case B, perceived to be the easiest benchmark, using the open-source simulator OPM Flow (The Open Porous Media Initiative). Building on the most highly resolved submission to SPE11, we conduct systematic refinement studies on facies-aligned corner-point grids, extending the resolution down to the submeter scale and beyond in a localized subdomain. Convergence is assessed not for dense reporting fields, but rather for the 13 sparse measures defined in SPE11, which act as proxies for operational risks such as overburden disturbance, leakage, dissolution dynamics, and migration into low-permeability seals. We find that pressures at observation points converge robustly even on coarse grids, whereas measures related to phase partitioning and dissolution exhibit strong grid sensitivity, and in several cases, no clear convergence, even at the finest resolutions considered. In particular, the onset and development of convective dissolution fingers remain highly resolution-dependent, underscoring the challenges of resolving small-scale processes in field-scale models. These results highlight that while some operationally critical quantities can be simulated with confidence, others remain elusive under straightforward grid refinement. The findings emphasize the need for complementary approaches—such as upscaled models or proxy measures—to ensure reliable predictions for large-scale CO2 storage operations.

Summary: We present a systematic investigation into how the choice of grid type and discretization scheme affects predictions in the 11th SPE Comparative Solution Project (SPE11). Focusing on Case B, we compare grid types with varying conformity with internal aquifer geometry: standard Cartesian grids, two types of cut-cell grids, perpendicular bisector (PEBI) grids, hybrid quadrilateral-triangle grids, and Delaunay-triangulated grids. These grids are coupled with state-of-the-art flux approximation schemes: The standard two-point flux approximation (TPFA) scheme, which is inconsistent on non-K-orthogonal grids, and three consistent schemes—multipoint flux approximation (MPFA), average multipoint flux approximation (AvgMPFA), and nonlinear two-point flux approximation (NTPFA). For transport terms, we evaluate a second-order weighted essentially nonoscillatory (WENO) scheme as an alternative to the standard single-point upstream-mobility weighting (SPU) scheme as a means to reduce numerical smearing and directional bias. Simulations are conducted using isothermal black oil and multicomponent K-value thermal flow models implemented in the open-source MRST and JutulDarcy simulators. (The thermal model assumes a constant thermal gradient for simplicity.) We also include results from Case C, emphasizing computational efficiency. We observe that variability in the prediction of overall plume migration and measurables used as proxies for assessing risk, such as pointwise pressure buildup, is modest across grid types and discretization schemes. However, the choice of grid type significantly affects the formation of the self-enhancing dissolution fingers that drive convective mixing beneath the CO2 plume—a key focus of the comparative solution project. These fingers form at locations where the discrete grid representation of the CO2-brine interface deviates from the true interface. The greater the number of deviation points, the more fingers will be triggered. Likewise, the more irregular the local representation of the interface, the stronger the fingers develop. In comparing different grid types and discretizations, we conclude that Cartesian grids with the standard TPFA/SPU discretization are likely the best choice for high-resolution simulations. Most simulators are optimized for this combination, and deviations from K-orthogonality are relatively minor for Case C and entirely absent in Case B. For simulations at more modest resolutions, we recommend using a grid type that better conforms to faults and facies boundaries, preferably a cut-cell grid. While such grids are generally not K-orthogonal, using a consistent discretization like AvgMPFA can help mitigate inconsistency errors that might otherwise compromise simulation accuracy.

Standard solution strategies in reservoir simulators all rely on a global linearization and the solution of large (ill-conditioned) linear systems, particularly for fully implicit formulations. Such solves are expensive, and typically require significant interprocess communication when implemented on a distributed memory architecture. Modern nonlinear domain decomposition (NLDD) is based on adaptive alternation between local and global solves and has proved to be an effective way to reduce the number of global linear solves in serial mode. Herein, we discuss the extension of NLDD to a distributed setting, exemplified by the operational reservoir simulator OPM Flow. We discuss which parts are reused from the standard Newton solver and which modifications are required to deploy the methodology in a standard simulator. We apply our new methods to a series of open benchmark and real-life asset models. In serial, NLDD benefits significantly from using a Gauss–Seidel approach in which subdomains are solved in order of decreasing pressure. This strategy can only be applied effectively within each process in a distributed setting, but not between processes, which potentially poses a challenge to scalability. However, we observe excellent improvement in efficiency of the nonlinear solution phase also in parallel, in particular when slight adjustments are made in the choice of linear solver. Our implementations are freely available as open source, enabling readers to reproduce many of the examples on their own hardware and test the methods on their own models.

Reservoir simulators are highly complex computer programs and are often treated as “black boxes” that act on standardized input formats. In part, this is due to the closed-source nature of commercial offerings commonly used in industry, but even when the source-code is available, modification of a fully featured simulator can be a daunting task: A reservoir simulation problem, when fully specified, can have a great number of parameters and functional relationships that are tightly coupled together in most implementations. This is an obstacle to workflows that go beyond forecasting and into optimization, parameter fitting and sensitivity analysis where the solution quantity of interest should be paired with gradients taken with respect to forces and parameters. The monolithic design of typical reservoir simulation codes cannot be separated from the technical and mathematical choices made when discretizing the governing equations: The implicitness required to handle pressure changes and strong coupling between different solution variables means that the Jacobian of the governing equations must be implemented with respect to the primary variables: An error-prone process when done manually, and highly time-consuming when considering coupled systems. Automatic differentiation (AD) is a common cure to this problem that can generate Jacobians, but often comes at a severe cost in runtime performance. A more practical issue is that high computational performance for simulation problems relevant to modern geoenergy operations necessitates implementation in a compiled language while optimization and history matching workflows that integrate many data sources are most naturally expressed in high-level scripting languages. We describe the design of JutulDarcy, an open-source porous media simulator written in the Julia language. JutulDarcy is designed from the ground up with extensibility and computation of sensitivities in mind, starting from the assumption that high-performance automatic differentiation enables us to radically rethink the way simulators are designed. Highlights of this design include arbitrary execution graphs for the equation terms, in which dependence relationships can be re-configured on the fly, easy coupling of models, computation of sensitivities with respect to all model parameters through an adjoint method with a parameter declaration system, and an assembly process with performance that can exceed that of compiled simulators with hand-coded Jacobians. We validate the simulator on a number of standard test cases (immiscible, black-oil, CO2 and compositional flow) and demonstrate MPI parallel performance on multi-million cell models. In addition, we offer a frank assessment of the relative strengths and weaknesses of the Julia programming language for scientific computing.

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.

Integrating microbial activity into underground hydrogen storage models is crucial for simulating long-term reservoir behavior. In this work, we present a coupled framework that incorporates bio-geochemical reactions and compositional flow models within the Matlab Reservoir Simulation Toolbox (MRST). Microbial growth and decay are modeled using a double Monod formulation, with populations influenced by hydrogen and carbon dioxide availability. First, a refined Equation of State (EoS) is employed to accurately capture hydrogen dissolution, thereby improving phase behavior and modeling of microbial activity. The model is then discretized using a cell-centered finite-volume method with implicit Euler time discretization. A fully coupled fully implicit strategy is considered. Our implementation builds upon MRST's compositional module by incorporating the Søreide-Whitson EoS, microbial reaction kinetics, and specific effects such as bio-clogging and molecular diffusion. Through a series of 1D, 2D and 3D simulations, we analyze the effects of microbial-induced bio-geochemical transformations on underground hydrogen storage in porous media.These results highlight that accounting for bio-geochemical effects can substantially impact hydrogen loss, purity, and overall storage performance.