Multiphase multicomponent flow in porous media with general reactions: efficient problem formulations, conservative discretizations, and convergence analysis

We develop a numerical framework for efficiently simulating partially miscible multiphase multicomponent flow in porous media with general chemical reactions, modeled by a system of coupled and strongly nonlinear partial differential equations, ordinary differential equations, and algebraic equations. The equations are transformed with the help of a reduction scheme which preserves the model and allows the unknown equilibrium reaction rates to be eliminated. In the transformed system, the algebraic equations resulting from the chemical equilibria are eliminated in terms of a nonlinear, implicitly defined resolution function. Hereby, the equilibrium conditions consisting of equations and inequality constraints are formulated as a complementarity problem and rewritten as algebraic equations using the minimum function. The existence of the resolution function is established by using the connection of the nonlinear system of algebraic equations to a constrained minimization problem based on a modified Gibbs functional, and numerical strategies to evaluate it numerically are presented. Using a global implicit approach, the nonlinear systems that remain to be solved in each time step are treated with the semismooth Newton method. With the help of the resolution function, we are ready to define a persistent set of primary variables that are valid in either phase state and for an arbitrary mineral assemblage. Different numerical tests related to hydrogen migration in deep geological repositories of radioactive waste and to CO2 sequestration show that our solver provides accurate numerical results, and that it is capable of handling the strongly nonlinear coupling of flow, transport, chemical reactions, and mass transfer across the phases. The second part of this work is concerned with the analysis of robust mixed hybrid finite element methods of lowest order for advection–diffusion problems. In the approach studied here, the classical scheme is extended with the help of the Lagrange multipliers associated with the hybrid problem formulation. More precisely, it relies on the fact that the Lagrange multipliers represent approximations of the scalar unknown on the interelement boundaries and uses them in the approximation of the advective fluxes. Using techniques from the a posteriori error analysis, we are able to establish optimal order convergence for a new class of methods including the standard method, a partial upwind method, and a full upwind method. The advantage of the specific choice of the upwind weights is the fact that the method remains local. In this case, static condensation can be employed, whereas the standard upwind-mixed method requires information from neighbor cells such that static condensation is not applicable. Concerning approximations with the BDM1 mixed finite element, we show that the use of Lagrange multipliers in the discretization of the advective term provides optimal second order convergence for the total flux variable in the L2 norm, whereas the standard mixed method is known to be of suboptimal first order accuracy only