In this chapter, we discuss how two-way coupled fluid flow and geomechanics can be modeled in the MATLAB Reservoir Simulation Toolbox (MRST) using the ad-mechanics module. A brief introduction to linear poroelasticity is provided, which is a common framework for studying geomechanics in the context of reservoir management or groundwater applications. We review commonly used poroelastic coefficients and moduli and present a handy tool that removes the need to manually navigate the large number of poroelastic relationships to compute values of needed parameters. The chapter further provides three examples where well-known model cases in linear elasticity and poroelasticity are modeled in MRST and compared with results from analytical estimates. These examples include the compression of a dry sample (a linear elastic problem) as well as the compression of a wet sample (Terzaghi’s problem) and Mandel’s problem.

In this chapter, we present fv-unsat, a multipoint finite-volume–based solver for unsaturated flow in deformable and nondeformable porous media. The latter is described using the mixed form of Richards’ equation, whereas the former by the equations of unsaturated poroelasticity. The module aims at flexibility, relying heavily on discrete operators and equations, exploiting the automatic differentiation framework provided by the MATLAB Reservoir Simulation Toolbox (MRST). Our examples cover two numerical convergence tests and two three-dimensional practical applications, including the water infiltration process in a nondeformable soil column and a realistic desiccation process of a deformable clay sample using atmospheric boundary conditions. The resulting convergence rates are in agreement with previously reported rates for single-phase models, and the practical applications capture the physical processes accurately.

The compositional module in the MATLAB Reservoir Simulation Toolbox (MRST) implements two different formulations of a three-phase compositional system that consists of a pair of multicomponent phases and an optional immisicible phase. In petroleum engineering, the aqueous phase is taken to be immiscible and the hydrocarbon liquid and vapor phases are governed by an equation of state (EoS). The overall composition formulation uses pressure and overall mole fractions as primary variables, whereas the natural variable formulation relies on solving for phase mole fractions and phase saturations simultaneously. Thermodynamic behavior is modeled using $K$-values or a (standard) cubic EoS. In the chapter, you will learn about the model equations, choice of primary variables, and numerical strategies for solving the thermodynamic problem, alone or coupled to the flow equations. We discuss details of the implementation, which builds upon the object-oriented, automatic differentiation (AD-OO) framework and utilizes state functions and generic model classes for increased modularity. We also present a few relatively simple simulation examples to illustrate typical behavior and teach you how to set up simulation cases yourself.

Simulation of multiphase flow in fractured reservoirs still poses a challenge due to the different timescales of fluid flow in fractures and matrix. Common approaches to modeling fractures in reservoir simulators include the discrete fracture and matrix (DFM) method, where the fractures are explicitly represented as lower-dimensional elements in the computational mesh, and multicontinuum approaches (e.g., dual-porosity and dual-permeability models) where the behavior of the fractures and matrix are integrated and treated as distinct continua. The latter requires models (bespoke “transfer functions”) that upscale the multiphase transfer between fracture and matrix. There are several formulations for transfer functions available in the literature, and they are often application dependent. Here, we propose a unified framework for simulation of flow in fractured media. The framework makes no distinction between dual-continuum and DFM methods, treating fractures and one or more matrix domains as flowing domains and virtual domains. Transfer functions are reinterpreted as fluxes between cells of different domains. This enables us to create an abstraction that encompasses both methods and makes it easy to build hybridized models including different regions with different matrix/fracture interaction concepts. We present a series of cases to illustrate the main differences between both modeling approaches and the benefit of a flexible implementation that enables the development of a fit-for-purpose simulator for fractured reservoirs.

This chapter explains how one can formulate nonlinear finite-volume (NFV) methods, as advanced discretization schemes, to solve the flow equation in porous media. These schemes are of particular interest because apart from being consistent, they are monotone by design. We explain the basic ideas of the NFV methods: how to construct one-sided fluxes, interpolate using harmonic averaging points, and obtain unique discrete fluxes through grid faces with convex combinations of one-sided fluxes. We outline key functions in the accompanied nfvm module in the MATLAB Reservoir Simulation Toolbox (MRST) and show some examples of how the method is applied.

The multiscale restricted-smoothed basis (MsRSB) method is the current state-of-the-art within multiscale methods. MsRSB is very robust and versatile can be used either as an approximate coarse-scale solver having mass-conservative subscale resolution or as an iterative fine-scale solver that will provide mass-conservative solutions for any given tolerance. The performance of the method has been demonstrated on incompressible two-phase flow, compressible two- and three-phase black-oil models, as well as compositional models. It has also been demonstrated that the method can utilize combinations of multiple prolongation operators; e.g., corresponding to coarse grids with different resolutions, adapting to geological features, adapting to wells, or moving displacement fronts. This chapter explains the basic ideas of the MsRSB method, including methods to construct coarse partitions, prolongation, and restriction operators; reduction of the fine-scale flow equations to a coarse-scale system; and formulation as part of a two-level iterative solver. We outline the key functions in the module and show various examples of how the method can be used as an iterative solver for incompressible and compressible flow on 2D rectilinear grids, unstructured grids, and 3D stratigraphic grids.

Geothermal energy is renewable, always on, and available anywhere (at least in principle). Hot underground aquifers are therefore appealing as a source of green energy but also for large-scale energy storage, which is important to buffer the seasonal energetic imbalance associated with the use of renewable energies. The viability of a geothermal exploitation project is determined by a number of factors such as energy efficiency, storage capacity, economical aspects (e.g., drilling and operational costs), and compliance with legal regulations. Such assessments require a detailed characterization of the geology and physical properties of both the aquifer and aquiclude, groundwater chemistry, and flow properties. Proper understanding of these processes depends on accurate and flexible numerical simulation tools. In this chapter, we present geothermal, a module for geothermal simulations of low- to moderate-enthalpy geothermal systems. The module implements the equations for conservation of energy and conservation of mass for water and salt (NaCl), along with pressure-, temperature- and NaCl-dependent viscosity and density and other functionalities specific to geothermal problems. We demonstrate the accuracy of the module by benchmarking it with TOUGH2, a widely used groundwater flow simulator. We also show how geothermal can be used to simulate different geothermal applications.

The ad-core module in MRST offers an object-oriented framework for rapid prototyping of new reservoir simulators based on automatic differentiation (AD-OO). The framework simplifies the task of changing and extending existing simulation models in MRST or implementing brand new ones. The MRST textbook (Lie, Cambridge University Press, 2019, Chapter 12) presents a model hierarchy for the black-oil equations, discretized by a standard fully implicit method, and describes how to (automatically) select time steps and configure linear and nonlinear solvers. Herein, we present a further modularization of the AD-OO framework that aims to simplify the implementation of more complex flow models and other types of discretizations and solution strategies. To this end, we view the reservoir simulator as a graph of functional relationships and their dependencies and introduce the new concept of so-called state functions to define these functional relationships and compute discrete quantities required for the linearized governing equations. Using the graph perspective, it is relatively simple to not only visualize and understand the data flow of highly complex reservoir simulators, but also replace components of the graph and/or extend the graph with new branches as needed. The result is a versatile family of reservoir simulators that can easily be configured to run different types of multiphase, multicomponent models and at the same time support a number of different spatial and temporal discretizations. The state-function concept also has a built-in compute cache that helps you to systematically eliminate redundant function evaluations. The chapter explains the new concept in detail and exemplifies its use by showcasing implicit, explicit, and adaptive-implicit versions of the same physical processes. We also demonstrate the use of consistent and high-resolution schemes to improve simulation accuracy. Applications to complex flow physics (EOR models, compositional flow, fractured reservoirs) are discussed in other chapters.

The upr module in the MATLAB Reservoir Simulation Toolbox (MRST) can construct unstructured Voronoi grids that conform to polygonal boundaries and geometric constraints in arbitrarydimensions prescribed inside the reservoir volume. The resulting volumetric tessellations are usually realized as locally orthogonal, perpendicular bisector (PEBI) grids, in which cell faces can be aligned to accurately preserve objects of codimension one (curves in 2D and surfaces in 3D) and/or cell centroids can be set to follow curves in 2D or 3D. This enables you to accurately model faults, let grid cells follow horizontal and multilateral well paths, or create lower-dimensional or volumetric representations of fracture networks. The module offers methods for improving grid quality, like configurable policies for treating intersecting geometric object and handling conflicts among constraints, methods for locating and removing conflicting generating points, as well as force-based and energy-minimization approaches for optimizing the grid cells. You can use \mcode{upr} to create a consistent hierarchy of grids that represent the reservoir volume, the constraining geometric objects (surfaces and curves), as well as their intersections. The hierarchy is built such that the cell faces of a given (sub)grid conform to the cells of all bounding subgrids of one dimension lower.

Many different factors contribute to elapsed runtime of reservoir simulators. Once the cases become larger and more complex, the required wait time for results can become prohibitive. This chapter discusses three features recently introduced into the MRST AD-OO framework to make simulation of large cases more efficient. By analysing the sparsity pattern of the Jacobians for some of the most common operations involved in computing residual flow equations, we have developed different implementations of automatic differentiation that offer better memory usage and requires fewer floating point operations. Using these so-called AD backends ensures (much) faster assembly of linearized systems. Likewise, these systems can be solved much faster by utilizing external packages for linear algebra; herein, primarily represented by the AMGCL header-only C++ library for solving large sparse linear systems with algebraic multigrid (AMG) methods. Last, but not least, the new ``packed problem'' format simplifies the management of multiple simulation cases and enables automatic restart of simulations and an ability for early inspection of results from large batch simulations. Altogether, these features are essential if you are working with bigger simulation models and want timely results that persist across MATLAB sessions.

We explain how you can use discontinuous Galerkin methods to formulate implicit higher-order discretizations of transport equations on stratigraphic and polytopal grids and outline how this is implemented in the dg module of the MATLAB Reservoir Simulation Toolbox (MRST).