# Multiscale Simulation of Flow and Transport in Porous Media

Flow in porous media is omnipresent. It exists in blood running through our capillaries, air moving through our lungs, humidity in the clothes we wear and sanitary products we use, fertile soil that sustains our crops, subsurface aquifers and filter systems that provide clean water, and deep geological formations that yield hydrocarbons and thermal energy. Modelling to describe and understand such systems is therefore significant. Here I discuss fluid flow in porous rock, which is relevant to the geosciences, climate studies, and various other areas.

Conservation of mass and Darcy’s law describe the flow of a single, incompressible fluid as follows:

\[\nabla \cdot \vec{v} = q, \qquad \vec{v} = - \mathbf{K} \nabla p, \tag1 \]

where \(\vec{v}\) is discharge per unit area (often known as Darcy flux or velocity), \(p\) signifies fluid pressure, and \(q\) represents volumetric source terms.

Alternatively, we can write \((1)\) as a second-order, elliptic Poisson equation \(-\nabla\cdot(\textbf{K}\nabla p) =q\). The positive-definite tensor \(\mathbf{K}\)—called permeability—represents the rock’s ability to transmit fluids, which is key to determining flow patterns in aquifers and hydrocarbon reservoirs. Sedimentary rocks develop over millions of years via a layering process that compacts different sediments into a highly complex structure. As a result, permeability exhibits a strong multiscale nature—with a wide range of spatial correlations that have no clear scale separation—and up to eight or nine orders-of-magnitude variations. This makes it surprisingly difficult to discretize and solve \((1)\) with sufficient accuracy, since important changes often occur over distances that may be significantly shorter than feasible grid sizes in flow simulators.

The mathematical geosciences community has therefore devised a number of upscaling or homogenization techniques that aim to capture essential features of the permeability on a coarse scale through a mixture of analytical and numerical methods. Over the past two decades, many researchers have specifically focused on so-called *multiscale methods* that attempt to solve reduced-order pressure equations on a coarse scale — using numerically-computed basis functions that consistently incorporate subscale permeability variations. If \((1)\) is discretized with a cell-centered, finite-volume scheme as linear system \(A \mathbf{p} = \mathbf{q}\) for fine-scale pressure \(\mathbf{p}\), a multiscale method instead seeks to solve a coarser problem. This occurs via the approximation \(\mathbf{p} \approx B \mathbf{p}_c\), where \(\mathbf{p}_c\) is a pressure that one can map onto the fine scale by basis matrix \(B\). Applying linear restriction operator \(R\) finds the coarse system

\[A_c \mathbf{p}_c = \mathbf{q}_c, \qquad A_c = RAP, \qquad \mathbf{q}_c = R \mathbf{q}. \tag2\]

Naturally, the approximation’s quality and efficiency strongly depend on the choices one makes to compute \(B\)’s basis functions.

This procedure is a specialized algebraic two-level solver and does not offer much computational gain over a standard algebraic multilevel solver if we only wish to compute the pressure once. The main interest in multiscale methods arises from multiphase flow simulations, illustrated as follows by the simple prototype problem of incompressible, immiscible, two-phase flow:

\[\nabla \cdot \vec{v} = q_p, \qquad \vec{v} = - \lambda(S)\mathbf{K} \nabla p, \qquad \phi \frac{\partial S}{\partial t} + \nabla \cdot (\vec{v} f(S) ) = q_s. \tag3 \]

Here, the total mobility \(\lambda\) models permeability reduction in the presence of multiple fluid phases, and \(S\) is the volume fraction of the wetting phase. To solve \((3)\) and advance \(p\) and \(S\) in time, one must first solve the Poisson problem \(-\nabla\cdot(\lambda \textbf{K}\nabla p) =q\) multiple times with variations in \(\lambda\), source terms, and boundary conditions. Additionally, \((3)\) frequently appears as a forward problem in the highly heterogeneous, data-dense problems that pertain to uncertainty quantification, forecasting, and optimization.

Many research papers have proposed and analyzed various forms of multiscale methods for solving \((3)\). However, a significant gap exists between the idealized, regular, structured grids in the literature and the highly-complex grid models that geologists use to describe real rock formations. The most widespread format in industry is a semi-structured design composed of hexahedra that may degenerate, shift up and down, or collapse to zero volume to account for erosion, faults, and other deformations in the deposited rock volumes.

**Figure 1.**Illustration of some of the challenges present when grids model geological porous media. Figure courtesy of Olav Møyner.

Scientists also frequently build grids from more unstructured elements; such grids may contain local refinements that adapt to well trajectories and important geological characteristics. One can normally assume that a typical grid will consist of general polyhedral cells, non-matching interfaces, large aspect ratios, multiple neighboring cells, and several orders of magnitude of variation in cell-to-cell permeability. Figure 1 depicts a few representative examples taken from engineering models. Posing local problems to produce acceptable basis functions on these grids is a substantial challenge that is essential to implementation of multiscale methods for various applications, including engineering computations in hydrocarbon recovery, carbon dioxide sequestration, geothermal energy processes, and natural gas storage. Furthermore, typical industrial applications use multiphase, multicomponent flow models that are significantly more intricate than \((3)\) and harbor further complexity and submodels that account for wells and various forms of surface and control facilities.

As a part of an ongoing industrial collaboration, we worked to adapt existing multiscale methods to more general grids. After several partial successes, we discovered a method [3] that was sufficiently robust on general grids. The approach—illustrated in Figure 2—is based on techniques from algebraic multigrid [4] and has a surprisingly simple implementation, as requested by our industrial partners. A quick turnaround—execution within a few weeks of delivering the method manuscript—was made possible by a very flexible software environment and continuous validation of models with full industrial complexity. We easily sketched out, prototyped, and tested new models in the same environment so that seemingly promising novel ideas could “fail fast” on real data and be discarded when needed.

**Figure 2.**Generation of basis functions through restricted smoothing. Figure courtesy of Olav Møyner.

This is typical when applying mathematics to issues in industry. While mathematicians tend to excel at solving difficult problems, targeting the *right* difficult problem is a challenge. Whereas simplification of a problem is an essential part of its solution, one must possess significant domain knowledge to delineate essential aspects from the more mundane details. Reinventing all parts necessary to ride the buggy is often time-consuming and impractical for smaller groups of researchers. Thus, the distance from research to industrial application may seem unnecessarily long. The computational geosciences group at SINTEF, where I work as a researcher, is committed to producing open-source software. For example, we freely release our research platform—the Matlab Reservoir Simulation Toolbox (MRST)—as a “batteries included” package for reservoir simulation [1-2]. So far, it has contributed to 140 master’s and Ph.D. theses and 270 scientific papers by external authors. We thus believe that MRST is a significant asset for the research community and a great starting point for scientists interested in testing their problems on porous media applications.

*This article is based on the author’s SIAM Activity Group on Geosciences Early Career Prize lecture, entitled “Multiscale Simulation of Porous Media Flow: Obstacles, Opportunities, and Open-source.” He presented this talk at the 2019 SIAM Conference on Mathematical & Computational Issues in the Geosciences, which took place earlier this year in Houston, Texas.*

**References**

[1] Krogstad, S., Lie, K.-A., Møyner, O., Nilsen, H.M., Raynaud, X., & Skaflestad, B. (2015). MRST-AD – an open-source framework for rapid prototyping and evaluation of reservoir simulation problems. In *2015 Reservoir Simulation Symposium*. Houston, TX: Society of Petroleum Engineers.

[2] Lie, K.-A. (2019). *An introduction to reservoir simulation using MATLAB/GNU Octave: User guide for the MATLAB Reservoir Simulation Toolbox (MRST)*. Cambridge, England: Cambridge University Press.

[3] Møyner, O., & Lie, K.-A (2016). A multiscale restriction-smoothed basis method for high contrast porous media represented on unstructured grids. *J. Comput. Phys., 304*, 46-71.

[4] Vaněk, P., Mandel, J., & Brezina, M. (1996). Algebraic multigrid by smoothed aggregation for second and fourth order elliptic problems. *Comput., 56*(3), 179-196.