The IMMIS mode applies to multiple completely immiscible phases. The code PIMS, parallel immiscible multiphase flow simulator, is a simplified version of the MPHASE mode in which the dependency on thermodynamic relations have been removed, since for immiscible systems the solubility is identically zero for each component. In this case the number of components is equal to the number of phases, or degrees of freedom associated with each node for an isothermal system. The immiscible property removes the variable switching strategy used in MPHASE, which may be the most numerically difficult part of PFLOTRAN, and may cause problems for multi-level solvers. The governing equations solved by PIMS are given by

(1)\[\frac{{{\partial}}}{{{\partial}}t}\big(\varphi\rho_{{\alpha}}^{} s_{{\alpha}}^{}\big) + {\boldsymbol{\nabla}}\cdot \big(\rho_{{\alpha}}^{} {\boldsymbol{q}}_{{\alpha}}\big) = Q_{{\alpha}},\]

where the subscript \({{\alpha}}\) denotes an immiscible phase.

In this equation \(\varphi\) is porosity, \(s_{{\alpha}}\), \(\rho_{{\alpha}}\) refer to the \({{\alpha}}\)th phase saturation and density, respectively, \({\boldsymbol{q}}_{{\alpha}}\) is the Darcy velocity of the \({{\alpha}}\)th phase given by

(2)\[{\boldsymbol{q}}_{{\alpha}}= -\frac{kk_{{\alpha}}}{\mu_{{\alpha}}} \big({\boldsymbol{\nabla}}p-\rho_{{\alpha}}g \hat{\boldsymbol{z}}\big),\]

with permeability \(k\), relative permeability \(k_{{\alpha}}\), fluid viscosity \(\mu_{{\alpha}}\), and \(Q_{{\alpha}}\) is the source/sink term. The selection of primary variables are pressure \(p\) and \(n-1\) independent phase saturation variables \(s_{{\alpha}}, {{\alpha}}=1,...,n-1\) with

(3)\[\sum_{{{\alpha}}=1}^n s_{{\alpha}}= 1.\]

The mass conservation equations are coupled to the energy balance equation given by

(4)\[\frac{{{\partial}}}{{{\partial}}t} \Big(\varphi\sum_{{\alpha}}s_{{\alpha}}\rho_{{\alpha}}U_{{\alpha}}+ (1-\varphi) \rho_r C_r T\Big) + {\boldsymbol{\nabla}}\cdot\Big(\sum_{{\alpha}}\rho_{{\alpha}}{\boldsymbol{q}}_{{\alpha}}H_{{\alpha}}- \kappa{\boldsymbol{\nabla}}T\Big) = Q_e,\]

where \(U_{{\alpha}}\), \(H_{{\alpha}}\) denote the internal energy and enthalpy of the \({{\alpha}}\)th fluid phase, \(\kappa\) denotes the thermal conductivity of the bulk porous medium, \(\rho_r\), \(C_r\) denote the rock density and heat capacity, and \(T\) refers to the temperature. Thus the number of equations is equal to number of phases plus one, which is equal to the number of unknowns: (\(p\), \(T\), \(s_1\), …, \(s_{n-1}\)).