Multiphase flow is a critical process in a wide range of applications, including carbon sequestration, contaminant remediation, and groundwater management. Typically, this pro- cess is modeled by a nonlinear system of partial differential equations derived by considering the mass conservation of each phase (e.g., oil, water), along with constitutive laws for the relationship of phase velocity to phase pressure. The problem becomes much more complex if the phases are allowed to contain multiple chemical species (also called components), as miscibility and phase transition effects need to be taken into account. The main problem with phase transition stems from the inconsistency of the primary variables such as phase pressure and phase saturation, i.e. they become ill-defined when a phase appears or dis- appears. Recently, a new approach for handling phase transition has been developed by formulating the system as a nonlinear complementarity problem (NCP). Unlike the widely used primary variable switching method (PVS), which requires a drastic reduction of the time step size when a phase appears or disappears, this approach is more robust and allows for larger time steps. One way to solve an NCP system is to reformulate the inequality constraints for the primary variables as a non-smooth equation using a complementary function

(C-function). Because of the non-smoothness of the constraint equations, a semi-smooth Newton method needs to be developed. Another feature of the NCP approach is that the set of primary variables in this approach is fixed even when there is phase transition. Not only does this improve the robustness of the nonlinear solver, it opens up the possibility to use multigrid methods to solve the resulting linear system. The disadvantage of the complementarity approach, however, is that when a phase disappears, the linear system has the structure of a saddle point problem and becomes indefinite, and current algebraic multigrid (AMG) algorithms cannot be applied directly.

In this work, we aim to address computational issues related to modeling multiphase flow in porous media. First, we develop and study efficient solution algorithms for solving the algebraic systems of equations derived from a fully coupled and time-implicit treatment of models of incompressible two-phase flow. We explore the performance of several precon- ditioners based on algebraic multigrid (AMG) for solving the linearized problem, including “black-box” AMG applied directly to the system, a new version of constrained pressure residual multigrid (CPR-AMG) preconditioning, and a new preconditioner derived using an approximate Schur complement arising from the block factorization of the Jacobian. We show that the new methods are the most robust with respect to problem character as de- termined by varying effects of capillary pressures, and we show that the block factorization preconditioner is both efficient and scales optimally with problem size. We then generalize the block factorization method and incorporate it into a multigrid framework which is based on the multigrid reduction technique to deal with linear systems resulting from the NCP approach for modeling compositional multiphase flow with phase transitions. We demon- strate the effectiveness and scalability of the method through numerical results for a case of

two-phase, two-component flow with phase appearance/disappearance. Finally, we propose a new semi-smooth Newton method which employs a smooth version of the Fischer-Burmeister function as the C-function and evaluate its performance against the semi-smooth Newton method for two C-functions: the minimum and the Fischer-Burmeister functions. We show that the new method is robust and efficient for standard benchmark problems as well as for realistic examples with highly heterogeneous media such as the SPE10 benchmark.