A simple, conservative, staggered grid-based pseudo-transient scheme implemented on GPUs to solve two-fluid phase flow in porous reservoirs

Understanding and quantifying the flow of multiple pore fluid phases is a key component to understanding many natural and engineered reservoir processes. A well-known example is the migration of hydrocarbon phases in sedimentary reservoirs, including natural migration, conventional production, and enhanced hydrocarbon recovery. As Carbon Capture and Storage (CCS) emerged as a potential strategy to combat climate change, understanding multiphase flow processes of supercritical CO₂ and the various pore fluids is at the forefront of scientific interest.

Our governing equations are based on the conservation of mass and momentum in two immiscible fluid phases. The fluid phases may be compressible or incompressible. We assume the inertial terms to be negligible in both phases and momentum transfer to happen in the Darcy-flow regime. Multiple fluid phase effects are introduced in the mathematical model by saturation-dependent relative permeabilities and fluid viscosities, which uniquely determine the total mobility, relative mobility, and fractional fluid flow curves. Here we consider capillary effects negligible. This results in two independent equations: an elliptic equation for fluid pressure and a hyperbolic equation for saturation.

We choose a staggered-grid-based finite-difference discretization. The elliptic fluid pressure equation is solved using pseudo-transient iterations. The hyperbolic saturation equation is solved using a conservative first-order upwind scheme. Multiple coupling and time-stepping options (e.g., explicit/implicit, first order/higher order) are tested against an analytical solution. The numerical scheme can be implemented on GPUs in a straightforward manner using the ParallelStencil.jl package in Julia. We will provide some examples of various reservoir applications, which will also be used as performance benchmarks. Finally, we will discuss the potential of the presented numerical scheme to be implemented in a broader THMC framework.