Reservoir simulation is traditionally based on the assumption that water is an inert phase, while hydrocarbon components split into oil and gas phases. This approach is usually reasonable when modeling conventional hydrocarbon recovery, but specific applications may require accounting for mass exchange between the water and hydrocarbon phases.

We here present the extension of our Graphics Processing Units (GPUs) compositional reservoir simulator (Esler et al., 2021) to support gas-water equilibrium. Specifically, the Søreide and Whitson equation of state (EOS) (Søreide and Whitson, 1992) was implemented to compute mutual solubilities of hydrocarbon/brine mixtures. The impact of salinity on phase equilibrium is accounted for, with salt being treated as an active tracer. The simulator uses a mass-variables formulation, meaning that little modifications to the construction of transport equations and Jacobian assembly was required; most of the required code changes are localized in the EOS module for the computation of component fugacities, and phase properties such as partial molar fractions and partial molar volumes.

Treating salt as an active tracer instead of defining a further pseudo-component has an important advantage with the Søreide and Whitson EOS. If salinity changes as in water vaporization processes, our choice ensures that flash iterations can still be cast as a Gibbs Minimization problem with salt being a constant parameter. On the contrary, salinity would change as flash iterations progress, jeopardizing the thermodynamic consistency of the phase equilibria. The overall reservoir simulation system of equations is still accurate to first order in time, at the cost of possibly slight volume imbalances at the end of converged timesteps.

The accuracy of the implementation with respect to conventional CPU ones is proved using a wide range of problems where hydrocarbon-water mass exchange play an important role in the physics of the recovery/storage process. In particular, we focused on CO2 sequestration in saline aquifers, where solubility trapping is a key mechanism.

A key conclusion of this work is that the extreme performance of GPU-based reservoir simulation naturally transfers to new fields of study, which is critical when modeling saline aquifers whose extent is an order of magnitude larger than that of typical oil and gas fields.

During injection of supercritical CO2 into saline aquifers or depleted gas reservoirs, the complex interaction of CO2 and impurities with reservoir fluids plays a very important role and can significantly alternate the injectivity. Brine evaporation into the CO2-rich phase can lead to salt precipitation which will reduce the effective permeability of the porous rock. A tangible cooling of the near-wellbore region due to the Joule-Thomson effect can lead to hydrate formation which will reduce injectivity even more. Complex phase behavior of supercritical CO2 with brine and hydrocarbon components in highly heterogeneous porous media accompanied by all these phenomena will strongly affect pressure distribution which is in turn related to mechanical risks.

In this work, we present a unified simulation framework for modelling near-wellbore effects induced by supercritical CO2 injection developed in the Delft Advanced Research Terra Simulator (DARTS) platform. This framework uses the Operator-Based Linearization (OBL) technique for incorporating all complex physical phenomena in a fully coupled fully implicit manner. A general multicomponent multiphase flash based on a combination of classic cubic equations of state (e.g., Peng-Robinson) for hydrocarbon/CO2-rich phases and an activity model for the aqueous phase is implemented. Hydrate phase behavior is modelled using a modified Van der Waals-Platteeuw hydrate equation of state. Formation dry-out and salt precipitation are incorporated by using the Element Balance approach coupled with thermodynamics. Thermophysical property correlations relevant to the thermodynamic conditions of interest are implemented and validated against lab experiments.

We demonstrate that all important physical phenomena, such as the Joule-Thomson effect, hydrate formation and salt precipitation can be effectively captured by the OBL approach. We use several existing numerical benchmarks to validate the accuracy of the developed framework in the dynamic representation of all these effects. The interplay between these complex phenomena and reservoir heterogeneity is demonstrated in an unstructured heterogeneous near-wellbore reservoir model.

We construct a representative single matrix block model surrounded by fully penetrating vertical fractures that analytically captures both fracture and block depletion with fracture-matrix mass transfer. This is done with and without internal fracture structure in the matrix block. We solve the 1D Green’s function for a fracture system enveloping a matrix block in terms of the time evolution of average fracture pressure. We, likewise, solve the Green’s function for the average pressure in a matrix block surrounded by a constant pressure boundary. Average pressure is then easily transferred by material balance into cumulative production or instantaneous flowrate. Primary variables in assembling the interacting systems model are the volume ratio, Vf/Vm, permeability ratio, kf/kx, and geometry, (a/b)(ky/kx), with the last term accounting for both block shape and permeability anisotropy. The single block model is extended to include the influence on depletion rate for matrix blocks with internal, sub-scale fractures by solving for the Green’s function average pressure with a net-zero flux discrete fracture model. Internal fractures are infinite conductivity features that allow for increased mass transport to the matrix boundary and more rapid depletion. The single matrix block model is migrated to one for heterogeneous systems using superposition and matrix block distributions. We illustrate the signatures in pressure and Bourdet derivative for homogeneous and heterogeneous models. The influence of subscale fractures can be fully captured in influence of kx, ky/kx, or effective block size and shape, provided fracture characteristics are upscaled using this or other approaches.