Global implicit solver for multiphase multicomponent flow in porous media with multiple gas components and general reactions

In order to study the efficiency of the various forms of trapping including mineral trapping scenarios for CO2 storage behavior in deep layers of porous media, highly nonlinear coupled diffusion-advection-reaction partial differential equations (PDEs) including kinetic and equilibrium reactions modeling the miscible multiphase multicomponent flow have to be solved. We apply the globally fully implicit PDE reduction method (PRM) developed 2007 by Kräutle and Knabner for one-phase flow, which was extended 2019 to the case of two-phase flow with a pure gas in the study of Brunner and Knabner. We extend the method to the case of an arbitrary number of gases in gaseous phase, because CO2 is not the only gas that threats the climate, and usually is accompanied by other climate killing gases. The application of the PRM leads to an equation system consisting of PDEs, ordinary differential equations, and algebraic equations. The Finite Element discretized / Finite Volume stabilized equations are separated into a local and a global system but nevertheless coupled by the resolution function and evaluated with the aid of a nested Newton solver, so our solver is fully global implicit. For the phase disappearance, we use persistent variables which lead to a semismooth formulation that is solved with a semismooth Newton method. We present scenarios of the injection of a mixture of various gases into deep layers, we investigate phase change effects in the context of various gases, and study the mineral trapping effects of the storage technique. The technical framework also applies to other fields such as nuclear waste storage or oil recovery.


Introduction
The climatic changes caused by carbon dioxide (CO 2 ) in the atmosphere is one of the presently most challenging problems in the world. In detail, not only CO 2 is dangerous This work was performed completely during the time when Markus M. Knodel was affiliated at the Chair of Applied Mathematics 1, Universität Erlangen-Nürnberg, and is not in any relation with any later institute where MMK was affiliated afterwards, except for the present affiliation of MMK, the GCSC at Frankfurt University, as the revision of this paper was performed when MMK was already affiliated there.
Markus M. Knodel markus.knodel@gcsc.uni-frankfurt.de; knodel@math.fau.de Extended author information available on the last page of the article.
for the climate, but as well other gases such as CH 4 , which are as well produced too much to date, in many cases together with CO 2 .
Therefore, several techniques are under development with the aim to capture CO 2 and other climate killing gases in an irreversible way.
A major focus lies on the injection of CO 2 and other climate killing gases into deep layers of rocks within geological formations. These formations typically are covered by a cap rock which is impermeable. Inside the deep layers, depending on the rock structure, CO 2 and the other gases experiences different chemical and physical processes. The injected gas migrates through the pore space. A substantial amount of CO 2 is chemically split into its components and carbon gets bound in different forms below the ground. Similar scenarios apply for the CO 2 usually accompanying gases such as CH 4 and H 2 S. Concerning the different classes of binding mechanisms, we refer to the literature [11,12,32,49]. The most save binding mechanism is mineral trapping, where the carbon part of CO 2 finally gets bound in mineral form as salt, so it cannot escape any more given that no further processes happen.
Our model inherits all forms trapping mechanisms, in particular mineral trapping which is our focus.
Hence, in this context, CO 2 is injected in liquid form (accompanied in most cases by other gases such as CH 4 ) at high pressure at deep layers, the consequences are highly complex fluxes of gas and liquid deep below the ground. The aim is to study when and to which part the processes in the pore space enable or disable permanent mineral trapping of CO 2 (and other climate killing gases). Numerical simulations are the only available but promising possibility to study the efficiency of trapping mechanisms (mineral trapping) for CO 2 (and CH 4 , H 2 S, . . . ) storage in deep layers. Mineral trapping mechanisms of CO 2 and related gases ask for multiphase multicomponent flow partial differential equations which incorporate diffusive and advective transport and general reactions [4,49]. General reactions imply not only kinetic reactions which are comparably easy to handle, but moreover also equilibrium reactions. The precise prediction of gas and liquid flow asks for the solution of large nonlinear coupled systems of diffusion-advection-reaction partial differential equations (PDEs), algebraic equations (AEs) and ordinary differential equations (ODEs).
Precise mathematical simulations ask for two preconditions: Advanced models for the description of the geophysical processes and advanced solution techniques in order to simulate the model in an efficient way. The choice of a suitable formulation of the equations and the choice of the solution method is important for efficient numerical solution.
Often, splitting methods are applied to compute reactive transport mechanisms. The splitting solver techniques often allow for the use of already given techniques and codes applied to the more complex problem. The computations are split in each time step into a flow problem and a reactive transport problem. Once the corresponding subproblems are solved, the physical quantities are updated. This method allows to use customized codes such as to couple codes for multiphase flow with such ones for reactive transport problems.
One approach is the so-called SNIA method, sequential non-iterative approach [8]. SNIA causes splitting errors, and therefore requires very small time step sizes to control these errors [58].
More promising is the SIA method, sequential iterative approach [37], where the subproblems are solved alternately within a time step to eliminate or reduce the splitting errors. Still, the number of iteration steps needed often is large, while the time step size has to be small to reach convergence, [29,52].
In order to overcome the principal restrictions of splitting sequential approaches, more and more the so-called global implicit approach (GIA) is applied, the case of coupled multiphase flow and reactive transport problems is getting more and more in the focus of GIA approaches [1,15,60,62].
Even though the solution of the equations at each single time step is more demanding compared to the SIA approach, the general disadvantages and limitations of the SIA approach do not apply any more. Roughly spoken, while the computation of each time step itself is more complex, the number of time steps itself decreases drastically. Anyhow, one does not have the limitations of the coupling of splitted problems. For a comparison of SIA and GIA, we refer to [3].
We apply the fully globally implicit PDE reduction method developed by Kräutle and Knabner ("PRM", published in former papers [28,35]) which enables to eliminate the nonlinear and global coupled part and in particular the unknown rates of the equilibrium reactions. These properties are based upon specific variable transformations. Separating the resulting remaining PDE/ODE/AE system into a global and a local system, we apply a nested Newton solver method [34]. The nested Newton solver enables fast and efficient computation of the dynamics of the system by means of the application of parallel solvers to the Finite Element discretized / Finite Volume stabilized [7] PDE system. The local problems can be solved perfectly in parallel. Our computations of the behavior of the concentrations of the different species of the multiphase multicomponent flow are highly resolved in space and time.
For the case of single-phase reactive transport problems, the PRM applied as most effective approach demonstrated within the GdR MoMaS reactive transport benchmark [9,10].
In a recent paper, the fully global implicit approach of the PRM was applied to the case of coupled multiphase flow and reactive transport problems [6]. Brunner and Knabner extended the PRM to the case of two-phase flow, where one gas was considered together in a fully coupled manner with an arbitrary number of species in liquid and solid phase.
In this paper, we extend the directly afore mentioned mineral trapping scenario [6] to the case of an arbitrary number of species in gaseous phase, that means we extend the PRM to the case of the injection and reactive transport of various species in gaseous phase into deep layers of porous media.
The appearance/disappearance of the gas phase and of the mineral phases lead to semismooth equations. In detail, we use an extended capillary pressure as persistent global variable along with complementarity constraints to obtain a unique and consistent description of the phase state. For the phase disappearance, we use persistent variables which lead to a semismooth formulation. These equations are solved with a semismooth Newton solver.
Also we do not need to apply a discontinuous switch of variables in the case that the gas phase appears or disappears. Moreover, we adapt the technique originating from a study of Neumann, Bastian and Ippisch [48] and integrated in our reduction technique framework as described by the paper of Brunner and Knabner [6] to the case of an arbitrary number of species in gaseous phase which are in equilibrium with their dissolved counterpart in liquid phase.
Within a major revision, extension and adaption of the code used in the course of the paper of Brunner and Knabner [6], we implemented our methods into the parallel M++ [61] based RICHYM++ framework caring for codereusabilty at a high level.
Namely, we present results for the case of the injection of various gas species, and we study phase change effects which appear within our simulations.
Note that at this stage, mechanical deformations are not included in our model, as this is a minor effect in the given scenario [14,32]. The incorporation of such effects is planned for further papers.
Our techniques allow for predictions in very complex scenarios of miscible multicomponent multiphase flow scenarios and are applicable also for similar cases such as nuclear waste storage and oil recovery. This paper is organized as follows: At first, we introduce the mathematical model, Section 2. Then we describe the application of the reduction method to the case of an arbitrary number of gases, Section 3. Next we consider the application of the semismooth nested Newton solver, the separation into local and global variables, and the resolution function, Section 4. The global variables defined are persistent and can be used as primary unknowns throughout the computation.
In order to consider simulation results achieved with our method, we study cases of examples of applications to practical scenarios, namely also the case for a reaction system related to mineral trapping in the context of the sequestration of CO 2 together with other gases, Section 5. Our examples of configurations are chosen such that they serve as proof-of-concept of our method using reliable coefficients and laws. Hence at this stage, we did not focus to select fully realistic scenarios but moreover scenarios where the complexity of the combinations in part even is higher as within geophysically realistic scenarios. In particular, our framework can be easily applied to more complex chemical reactions.
In Section 6, we discuss our GIA for miscible multiphase multicomponent flow in porous media with multiple gas components in porous media, we discuss the potential of our approach, and we discuss our techniques in the context of other approaches. In the concluding Section 7, we we give some concluding remarks concerning our GIA in the context of present and future research in the field.

The PDEs for three phases with different gas species-flow, transport and reaction
In the following, we describe a purely macroscopic model (also called on the Darcy-or RVE-scale), which indeed may be derived by upscaling with various techniques, generality and mathematical rigor. For classical results based on periodic homogenization we mention e.g. [30]. We aim at unrestricted generality concerning flow, transport and reaction, leaving out two aspects: the mechanical deformation of the aquifier (see introduction) and the alteration of the porous medium due to the ongoing surface reactions. The latter can be included at various levels of complexity (as pure macroscopic or as micro-macro model, see e.g. [22,47]) and will be the subject of a subsequent paper.
Denoting α = g, , s as the gaseous, liquid and solid phase, we consider an arbitrary number of components in gaseous, liquid and solid phases described by the molar concentrations in a first form from mass conservation per phase where sorp means nonmineral and min mineral, and the indices range for each phase α = g, , s given by The PDEs describing multicomponent multiphase flow in porous media (with porosity φ) with general reactions, i.e. the reactive transport, read for I α concentrations in a first form from mass conservation per phase Of course, the solid phase species are macroscopically immobile and thus, q s and j s vanish.
Note that we use different colors to highlight important terms.
For all phases, we have the porosity φ of the porous medium in Eq. 4. As discussed, this is a given function in space.
For gas and liquid phase α = g, , we have -the molar density ρ mol α -the saturation of the phase α, termed s α . We have gas saturation s g and liquid saturation s . note that Note that the solid concentrations a measured per liquid volume, the same holds for gas concentrations. This means that the fluid-filled part of an REV serves as reference for "solid" and "gas" concentrations. -Darcy velocity q α (noted here with red color), which describes the advection in liquid and gaseous phase α = g, , and is given by the generalized Darcy's law with g the gravitation constant (vector), p α the total pressure of the phase α, and ρ mass α the mass density of the phase α. Further we have the absolute permeability tensor K, the relative permeability k r α of phase α, and μ α is the viscosity of the phase α. The Darcy velocity appears further in nonlinear form in the -Diffusion-dispersion-tensor (noted with teal color). The diffusion is described by Millington and Quirk [46] and the dispersion is described by Bear-Scheidegger [2].
with the transversal and longitudinal dispersion coefficients A L , A T for the mechanical dispersion, and the molecular diffusion coefficient D α for each phase α, indicating that we assume the same molecular diffusion for each component of each respective phase. The latter assumption is an approximation that can be justified if dispersion dominates diffusion. The diffusiondispersion-tensor enters the -diffusive flux j α (noted also with teal color, as directly related to the diffusion-dispersion tensor) with mole fractions χ α = c α /ρ mol α for each phase α.
For all states (g, , s), the right hand side f α describes the rate of kinetically and in equilibrium described reactions (noted with magenta color, not to interchange with the red color of the Darcy velocity) which are part of the chemistry of the system:

Chemical reactions
The complete system of equations does not consist only of the PDEs in Eq. 4, but in addition of ODEs and AEs. These ODEs and AEs result from the equations of state (EOS) for the gaseous and liquid phase, the equilibrium reactions including exchange between phases, and relations between different components of the equations, such as the relations between partial pressures of the gas and the total gas pressure.

Source term with reactions
Our model includes both kinetic and equilibrium reactions. The stoichiometric matrix of the chemical reactions reads as: which has the dimension (I g + I + I s ) × J where J = J eq + J kin (11) with J eq the number of equilibrium reactions and J kin the number of kinetic reactions, and I α the number of components in phase α. Each column of S stands for one reaction, and each row of it represents a chemical component, in its phase. We may write the source/sink term of the PDEs in Eq. 4 with the rates measured in moles per liquid volume, as in Eq. 9, or in matrix form We denote all mass exchange between gas and liquid phase as "reactions", i.e., they are incorporated in the vector of rates R. Further our model includes mass exchange between liquid and liquid and liquid and solid phase, described also as reactions incorporated in the vector of rates R.

Interphase mass transfer
In case of a partially miscible multiphase multicomponent system, we have to model the interphase mass transfer between the fluid phases. This effect is of major importance for mass and volume of the gas phase. The local appearance and disappearance of the gas phase can strongly depend on on this effect. A simple realization for exchange g ↔ (gas -liquid) is Henry's law [27] (which is vector-valued): with the Henry constant H (diagonal matrix valued) and partial gas pressure p g (vector valued partial pressures). For CO 2 exchange, other laws such as the experimentally based spline description derived by Spycher-Pruess [57] should be used (however only valid for CO 2 ; we use Spycher-Pruess for CO 2 ). Note that in principle, physical gas pressures only appear once the gas phase appears and the gas saturation gets bigger then zero. In order to account for this request, in later steps, we will introduce a more generalized definition of the gas pressure, which will allow the use of an equation similar to Eq. 13 in a straightforward way.

Chemical reactions-Law of mass action (LaMA)
In most cases, the activity and the concentrations are related as a i (c) = γ i (c)c i with an activity coefficient γ i . For the case of not too large concentrations, the approximation γ i ≈ 1 is justified. Therefore, later on, in general, we will replace the activity by the concentrations itself except for the main component in liquid phase, which is water, and also not for minerals. For water and minerals, the choice a i = 1 is thermodynamically indicated, therefore we use it.
The law of mass action (LaMA) is used for equilibrium descriptions of reactions in and between liquid and solid phase, with vector of equilibrium constants K and stoichiometric matrix S.

Kinetic reactions
In case of kinetic reactions and assuming that only homogeneous reactions between components of the liquid phase and heterogeneous reactions between components of the liquid phase and nonmineral components of the solid phase participate in kinetic reactions, a kinetic version of the LaMA reads for each reaction j , where the positive parameters k f j and k b j describe forward and backward rate constants.

Equilibrium reactions
We may reformulate the LaMA as with the equilibrium constant K j > 0 for the j th component of the equilibrium constant vector K.
Note that the reaction rates R eq j (c) are additional unknowns, to be eliminated.
For the exchange in and between liquid and solid phase, using ln(c ) = (ln(c 1 ), . . . , ln(c I )), (16) this equation specifies to: -Aquatic (mobility) reactions between different aqueous species ↔ : -sorption reactions between aqueous species and nonmineral solids: ↔ s sorp : (18) and -mineral reactions between aqueous species and mineral solids ↔ s min : Φ min , where the mineral reactions are particular. (We assume that one and only one mineral participates in each mineral reaction.) To account for the mineral equilibrium reactions, we define the function To cover that both the saturated case and the unsaturated case we use the semismooth equation, or, alternatively, a complementary condition, or rewritten as vector equation Concluding, we have the equilibrium conditions: Having in mind that the solid-liquid reactions, here described as volume reactions, are actually upscaled microscopic surface reactions, we expect specific surface factors in the reaction rates f for a kinetic description. In the multiphase case here, these factors rely on microscopic information, not available here (see e.g. [47]). To render them constant and incorporable in e.g the stoichiometric constants, requires structural assumptions, as e.g. only the liquid phase is in contact microscopically with the solid phase, which we will assume in the following.

Summary of the reactions
In detail, we have the following details of the different reaction types corresponding to the stoichiometric matrices and reaction coefficients: This indicates several restrictions, which can be omitted easily: There are no homogeneous reactions in the gas phase, there is no exchange gas-solid (see above). There are no kinetic reactions between gas and fluid phase, there are no homogeneous reactions in the solid phase, and mineral reactions are all equilibrium.
The reaction rates structure read as follows: Ane we have the total number of equilibrium reactions J eq = J ex + J mob + J sorp + J min (27) and the total number of reactions Finally, the source/sink terms may be rewritten as

Equations of state-EOS
The mass or molar densities of each phase can depend on various factors. For the gas, a usual approach is to use the equation of state (EOS) (30) where f g (p g ) can be considered in general as a generic compressibility law which needs to be further specified for simulations. The most simple version is the ideal gas law. Experimentally based splines such as the EOS of Duan-Moeller [16] (strictly speaking: only pure CO 2 injection) are the most realistic ones. We use them for CO 2 injection.
For the EOS of the liquid, the molar density of the liquid phase usually can be described as a function of the concentrations of its components, using a generic compressibility law f χ to be further specified. (For example, also constant molar mass density of the liquid phase in principle could be used instead.) However, more involved experimentally based spline descriptions exist, namely Garcias law if we consider only CO 2 exchange. The function f is given in case of Garcias law [20] by (32) and is given by a spline.

Capillary pressure
The capillary pressure defined as the difference between liquid and gaseous pressure Physically, p g and p c are only meaningful if the gas phase exists (s g > 0). An essential part of macroscopic multiphase fluid modeling is the assumption that p c is a function of (one of the) saturations We introduce extended (persistent) variablesp c ,p g .
both for s g > 0 and s g = 0.
As long asp c is below the entry pressure threshold, there is no gas phase, then we havẽ where the (nota bene) function p c (s g ) is given by the used model of soil water retention. Within our implementations, we use the Brooks-Corey [5] and the van Genuchten model [59]. As long as Eq. 36 holds, the liquid is undersaturated. Once the entry pressure has been exceeded, we have a gas phase and the relation between capillary pressure and saturation reads where we are now in the saturated state.
Combining these equations, we can write ). (38) This complementary problem can be rewritten as As Eqs. 23, 39 is not differentiable everywhere, but still semismooth [31].

Final basic equation system
Defining the transport operator L α describing transport due to advection, dispersion and diffusion for gaseous and liquid phase, we arrive at the final equation system Note that for the sake of clearness, we have inserted the entire arguments of the exchange reaction condition Φ ex (p g , χ ) to account for the consistent variable choicep g rather then p g , to extend (13).
For the sake of clearness, we further add the superscripts partial and total for the gas pressures.
In order to highlight the still present and problematic equilibrium reaction coefficients, we mark them with magenta color, the same is valid for the equilibrium conditions (which itself are not problematic). The new introduced transport operator is marked with blue color.

The orthogonal matrices
The aim of the PRM [28,35] is to reduce the number of nonlinearly coupled PDEs and in particular to remove the rates of the equilibrium reactions to enhance efficiency. A major property of the PRM compared to other methods such as the direct substitution approach (DSA) [45,52] is that the number of nonzero entries in the Jacobian is much smaller within our approach, because the transformation we apply does not create additional nonlinearties [28,29,35]. in opposite to approaches such as the DSA. This reduces strongly the computational costs for computing and inverting the matrix.
Starting with S α = S α A α to extract the linear independent part S α of the stoichiometric matrices S α of each phase α, the linear dependent parts are expressed with the A α matrices [29]. By assumption, we allow linear dependencies only among kinetic reactions, thus we have where consist of a subset of the columns of S kin respectively S kin s . As described above, there exist or can be created matrices A α , α = g, , s with and they have the block structure where Let S ⊥ α consist of columns that form a basis of the space that is orthogonal to the column space of S α . Hence, the orthogonal property holds. Additionally, matrices B α , α = g, , s are constructed with the same size as S * α with the property such that the columns {B α , S ⊥ α } form a basis of the complete space. A standard choice is in the following B α = S * α .

Multiplication of the equations, variable transformations, reduction procedure
We will first explain in brief the principle of the reduction procedure for a small subset of the PDEs, and afterwards apply this principle to the complete set of PDEs in order to get rid of the equilibrium reaction rates.

Principle of the reduction procedure
In order to explain the idea of the transformation, let us focus on Eqs. (41) and (42) at the moment. Multiplying Eqs. (41) and (42) with and afterwards transforming the concentrations to reaction invariant η and reaction participating ξ components transfers the PDEs (79) The afore indicated subtractions in Eq. 79 are the core of the PRM. The transformation and subtraction procedure leads to a strongly reduced number of PDEs, as the unknown equilibrium reaction rate appears in exactly one equation, which can be dropped to the expense of the equilibrium reactions. The sub-system Eq. 79 reduces to:

Transformations of the complete system
Now let us consider all the PDEs in Eqs. 41-43 and perform the above transformations (see [35]), details can be found in [33]. Now the new variables may be introduced as: We split the ξ 's into subvectors Using the transformed reaction variant respectively invariant variables, the PDE system in Eqs. 41-43 can be rewritten in terms of the new variables in Eq. 82 At first, the unknown R mob and Eq. 88 can be dropped because the reaction rates of R mob appear only in this equation. (as our system is overdetermined, we may drop some equations).
Then R ex , R sorp and R min drop out in the differences of these equations and the Eqs. 85, 93, and 94 can be dropped again for the same reasons as in the case of Eq. 88.
For reasons due to numerical robustness explained in [6,28], we apply the additional transformation using additional variables and, in order to be along with the equations implemented by Brunner [6] in the special case of only one gas in gaseous phase, we further substitutê Setting the system reads Finally, we recall the back transformation to get the entire physical concentrations where c

Global and local system and variables, resolution function, nested Newton
Many of the variables are spatially local. We can take advantage by eliminating these variables from the spatially connected global system by a nonlinear Schur-complement approach in terms of a so-called resolution function for the local system. We "split" our huge equation system Eqs. (99)-(121) into a "global" and a "local" problem, with "global" (primary) and "local" (secondary) variables. Both systems are solved by means of applying Newton solvers. Hence, we apply a nested Newton solver which facilitates the solution of the equation system strongly and increases efficiency substantially.
Our strategy ensures the global implicit strategy and is the basis of numerical stability.

Nested Newton and resolution function
We split the equation system Eqs.
One can consider the local variables as function of the global ones with the help of the spatially local equations which define the resolution function by means of the (approximate) solution of L(Ξ glob , Ξ loc (Ξ glob )) = 0.
It is a general and implicitly defined version of the constitutive laws, its (approximative) solution, generalizes the "flash calculations".
Thus the system reduces to This coupling causes that the Jacobian of the global system is computed as Therefore, the global Jacobian depends on the local system via the resolution function.
which is used to compute J glob Eq. 130. Note that the system of equations to be solved for its evaluation are small compared to the global system as the coupling (over the species) at most stems from the single reactions. Note that we write the global equations and variables with magenta color, the local ones with blue color, and the coupling term via the resolution function with red color.
The application of two Newton procedures, one to solve the global equations, and one to solve the local ones, coupled by means of the resolution function, is called a nested Newton solver.

Local and global system and variables
In order to apply the afore introduced nested Newton solver, we split the remaining equation system and the variables into local and global parts.
An intuitive natural splitting would be to choose the global system consisting out of the PDEs remaining after the reduction procedure, whereas the local system consisting out of the ODEs and the AEs. However, due to investigations concerning the existence of the resolution function which was elaborated in the predecessor publication of some of ours [6], a few AEs have to be shifted to the global system. Still, we had to enlarge the system elaborated by our predecessor work, in that we have more equations due to the introduction of various gas species.
We display the separation we use, we show information on the resolution function, and we discuss the changes indicated by the challenges in case of multiple gas phases in gas phase.
In the forthcoming, those equations which are new compared to the Brunner/Knabner case [6] are indicated in red.
At first, we show the new form of splitted equations and variables.
The global system reads: and the global variables are: The global system consists of the PDEs, the liquid EOS, and the gas exchange equilibrium reactions (Henry's law and/or Spycher-Pruess). New (not part of [6]) equations are indicated in red.
The local system contains the law φ cap for the capillary pressurep c , various relations between the variables, and the LaMA based equilibrium reaction conditions for liquidsolid interactions.
The local system reads: and the local variables are: For this choice, the local variables are uniquely determined by the global variables in terms of the resolution function. As this is possible regardless of the phase state and mineral assemblage, the global variables are persistent and no switching procedure is necessary if the gas phase appears or disappears or if a mineral precipitates or dissolves.
We apply a semismooth Newton solver [34] to solve both systems of local and global equations and variables (cf. former papers [6,28,35]). The Newton iteration for the local problem (at each grid point) is nested in the Newton iteration for the global problem.
As the local equations do not incorporate spatial couplings, they can be solved in parallel with perfect weak and strong scalability (if one ignores the global system, nevertheless in total strong speedups are possible due to this property in parallel).
The dependency of the local to the global variables ∂Ξ loc ∂Ξ glob based on the resolution function in Eq. 127 is indispensable to keep a Newton's method in total, otherwise the procedure would reduce to a nonlinear block Gauss-Seidel method.

Algorithmic differences compared to Brunner/Knabner study [6]
Compared to Brunner and Knabner [6], there were important challenges and extensions necessary to extend the case of one gas in gaseous phase to the case of an arbitrary number of gases in gaseous phase which are in equilibrium interchange with their liquid counterpart.
In particular, note that the equation type for the gases in gaseous phase changes. Whereas in case of only one gas, the gas PDE is hyperbolic, as there is no diffusion, in case of several gases, diffusion appears also for the gases in gaseous phase, so the gas PDEs now are of parabolic type.

New variables, new equations
Because in [6], the gas molar density could be identified with the gas concentration, there was no need to apply the reduction procedure for the gaseous phase PDE, and to define reaction variant and reaction invariant unknownsξ g , η g . Partial gas pressure and total gas pressure were identical, the mass density of the gas just was the product of the molar density with the mass.
In the case of various gases, the reduction procedure has to be extended also to the gaseous phase PDEs. Several new variables need to be introduced, we need to introduce new equations for the new relations between different new variables. And several equations which before were scalar valued now are vector valued.
The new or modified variables arise from the distinction between partial and total gas pressure, new definitions of molar and mass densities of the gas and the reaction variant and invariant unknownsξ g , η g arising from the reduction procedure applied to to the gas phase PDEs to allow for removal of the unknown equilibrium reaction rates.
Consequently, these new equations need to be distributed in a consistent way to the local and global equation system.
The new variables are ⎛ After the reduction procedure, the new equations originating from the PDEs and from the new variables read: To explain their origin, we repeated these new equations derived above.
The distribution of the new equations to local-global system is motivated as follows: A part of the new equations before was present as scalar equation already, but are vector valued equation now. This holds true for the gas liquid exchange law equations (Henry's law, Spycher-Pruess, . . . ) and for the relation of the reaction variable ξ andξ g , which equals the molar gas density in case of only one gas.
While the equation for the reaction invariants η g drop in case of vanishing equilibrium reactions between different gases in gaseous phase, the input values for the gas EOS effectively changes, depending on the gas EOS used. Further, now we need to relate partial gas pressures and total gas pressure, we get new equations relating molar and mass density of the gases in gaseous phase, and (what is not evident at the first glance) we need to introduce a relation equation between the partial pressures, the molar density of the gas, the concentrations of the different gases in gaseous phase, and the total gas pressure.
Concerning the distribution of these equations between local and global system, it was possible to put the major part to the local system.
The equations for the exchange between gas and liquid describing the relation between partial gas pressures and liquid molar densities now are vector valued, and we put these equations to the global system. The PDEs for the relation between gas and liquid exchange reaction variants necessarily have to be part of the global system, and change from scalar valued to vector valued.
The rest of the new equations are all local and it is possible to put them to the local system.
Note that due to the new equations, also the resolution function changes, and the corresponding contribution to the global Jacobian originating from the resolution function had to be adapted. Also the local system solution procedure need to be adapted.
We have performed comparison computations of our new algorithm with the results of [6] for the case of only one gas, i.e. CO 2 , to ensure consistence, and the results matched (data not shown).

Extended Modified Newton method
The relation between the generalized capillary pressurep c and the gas saturation s g , as given in Eqs. 33-39, is known for causing convergence problems for Newton's method. The reason is that the graph of the s g -p c -relation, cf. Eqs. 36-37, shows a kink at the point (s g ,p c ) = (0, p c (0)), i.e., between the branches s g = 0,p c ≤ p c (0) and s g ≥ 0, p c = p c (s g ). When at a certain computational point the gas saturation changes from zero in the previous timestep to s g > 0 in the current timestep, this kink causes a strong overestimation ofp c which disturbs the values of other variables and leads to nonconvergence of Newton's method. To overcome the convergence problems, a modification of Newton's method was proposed in [48] (cf. Section 5.5.2 p. 49-50), by resetting the overshooting value forp c to a certain threshold when necessary. This modification was adapted in [6] (cf. p. 135) to the situation of nested Newton solvers (local and global equations/variable), but up to now it was only established for the case of one gas species in the gaseous phase in interchange with its liquid counterpart. Indeed, [6] also had to adapt the correction for the gas pressure to ensure robustness due to the context of the nested Newton's method solver.
Since this issue is a crucial bottleneck of the complete solution procedure, it was important part of our work to find out how to generalize this technique to our scenario of several gaseous species, i.e., to scenarios where the partial pressures of various gases play a role. Indeed, the application of the resetting technique top c applies as before. However, the adaption of the partial pressures of the gases is a nontrivial question. We solved it by resetting all partial pressures with the same reduction factor to establish consistent variables corrections. It turned out that this settles the convergence issues, but we are aware of the fact that this solution still has some heuristic elements, and we intend to develop more advanced techniques for the gas phase appearance in the context of various gases, to allow for even better robustness while avoiding variable switches.

Notes on possible extensions
The method can be made more flexible, by approximating the Jacobian by elimination of certain off-diagonal entries which may correspond to fix-point handling of certain non-linearities. One may expect a deviation in order of convergence or even loss of convergence, but as long as convergence is preserved, the limit is the correct solution. This might be applied e.g. for more generalized activity coefficients or the alteration influenced porosity.

Simulation results
In order to demonstrate the efficiency of our extension of the framework presented by Brunner and Knabner [6] for one gas species, we present the results of two example computations with several gas species that get injected into the porous medium.
Simulation examples-overview The first example is comparably simple, as we consider the case of two injected generic gases dissolving in water which obey Henry's law. The only liquid components present are those in equilibrium with their two gaseous counterparts, and water. Solid components are not part of this example, but it demonstrates that our reduction method works also for more than one gaseous species.
The second example displays the case of three gases in equilibrium which are injected into the porous medium, obeying Henry's law. In this example, the injection time is not equal to the simulation time in contrast to the other examples of this study. The aim is to demonstrate the effects which appear after the stop of the injection, and to study the effects of the phase transitions for the case of various gases. The third example is highly complex, as we use the example Brunner and Knabner published for the case of one gas, CO 2 , but we enrich it with two additional generic gases (which might be, for example, CH 4 and H 2 S), and chemistry corresponding to these additional gases which are in equilibrium with their aqueous counterpart, but participate also within liquid to liquid aquatic reactions. The dissolved components further are in exchange with their counterpart in solid mineral form. This example enables to study the efficiency of our method for the case of a mineral trapping scenario in case that several gases are injected. As our study has the aim of a proof-of-concept study, we use Spycher-Pruess as exchange law for CO 2 , while the other two gases obey Henry's law. In all cases, we use the law of Brooks-Corey for the relation between the saturations and the partial pressure.
The parameters common for all three examples are listed in Table 1. Note that whereas in general, the absolute permeability is a tensor, we assume this tensor to be a multiple of the unit matrix, i.e.
reflecting the assumed macroscopic isotropy of the medium. This allows us better isolate the effects of the gas capture mechanisms from other flow effects.

Discretization methods and solvers
The PRM in principle can be discretized with any kind of discretization method. In our study, the discretization is performed by means of an adaptive implicit Euler method in time and Q1 Finite Elements where we use a Finite Volume upwind stabilization for the advective parts [7], which is equivalent to a standard cell-centered finite volume method [6] in case of a linear transport equation with a cellwise constant diffusion coefficient.
To solve the equation system arising from the global Newton, we use BiCGStab as linear solver with a SuperLU preconditioner. The global Jacobian is computed by means of analytic differentiation, also concerning that part referring to the resolution function contribution.
Caring for code-reusabilty, we implemented the afore explained algorithm into our parallel M++ [61] based RICHYM++ framework.

Physical time simulated and time step size
Our total simulation time of the physical process is 6 · 10 6 s, i.e. around 70 days, for all three examples we present.
The adaptive time step size choice strategy is the same as Brunner and Knabner used [6].

Domain and boundary conditions
The simulations presented in this study were performed at a 600 × 100 m sized 2D computational domain (cf. Fig. 1). The coarse grid of the discretized domain consists of rectangles of size 20 × 20 m (squares). (For the computations, we perform several refinements.) The upper boundary is 800 m below ground, the lower one 900 m.
The initial value of the liquid pressure is given in both examples as the hydrostatic profile p 0 (x, z) = 10 5 + (900 − z)ρ std mass, · 9.81 Pa.
Further, in all three examples, the initial concentrations are constant allover the computational domain.
Denoting by The boundary values of all other variables follow directly from the transformations performed along the reduction procedure.
Concluding, this means that we have For the nonzero influx values of the gases in liquid phase, we refer to the values given for the simulation examples presented below.
The Dirichlet boundary values at D are the same as the initial values for all components.

Injection of two gases into a deep layer
In this simulation example, we consider the injection of two gases in liquid form which are in equilibrium with their gaseous counterpart. The porous media domain is filled with water, but no minerals or any other chemically interacting components are present in the geological configuration.
Hence the chemistry under consideration reads: where the non-interacting water in liquid phase is noted to recall for its presence. The parameters we use are displayed in Tables 1 and 2, and the relation between liquid molar fractions and partial pressures as indicated by Henry's law and the parameters we choose for this simple example are displayed in Fig. 2.
In this simple case, we use a two-fold refined grid. 2 We display the results at two time snap shots: for one day and 70 days of physical time, the latter at the end of the simulation. Figures 3 and 4 display these simulation results.
We observe that soon after the injection of the gases, a gas phase appears. This gas phase moves upwards and accumulates below the impermeable upper boundary. Of course, as we do not have any solid components in this scenario respectively components choice, no mineral trapping scenario can appear.
We also consider evaluations of concentrations along trajectories to study the simulation results in more detail, cf. Fig. 5. An interesting observation is that that gas/liquid with highest Henry constant (slope for the relation between partial pressure and liquid mole fraction) shows increase at interface. Even though the peak width decreases slightly with increasing grid refinement level, the effect itself is stable also for higher grid levels. As well, Fig. 5 gives interesting insights into the phase change effects at the gas front between that region where the gas has already arrived, and that region where the gas was not transported so far by means of the advection dominated processes of our examples. (Note that we use the notions "gas front" and "interface between gas phase presence and gas phase absence" as some sort of synonym, because the simulation screenshots show a relatively sharp "jump" of gas concentrations, i.e. there is not a smooth transition from blue (low concentration) to red (high concentration), but  this happens within a small spatial zone, the "interface"-in opposite to diffusion dominated cases.) Anyhow, as this example is more for qualitative demonstration purposes, we leave exact studies for higher grid levels to the third, more complicated example, in order to study the complex questions at a complex example. The results shown here are taken from level 2 of grid refinement.

Injection of three gases with limited injection time
In the second example, we consider the case of the injection of three gases into the porous media domain. The three gases obey Henry's law and are in equilibrium with their liquid counterpart, where the gases are dissolved in water, but water is chemically inert. We call the gases CO 2 , CH 4 , and H 2 S, but mention that in part, the quantitative values of the properties of these gases do not coincide with their physical ones.
The injection time is 3·10 6 s, i.e. about 5 weeks of physical time, then injection stops, and the total physical simulation time is 6·10 6 s, i.e. about 10 weeks. This scenario allows for specific observations concerning the phase transitions.
The reaction chemistry we assume is as follows: The parameters specific to the given example are listed in Table 3.
The assumed, heuristic relations between the partial pressures of the gases in gaseous phase and the molar fraction of their corresponding counterpart in liquid phase are displayed in Fig. 6, i.e. the Henry's law relations.  It is clearly observed that after switching off of injection, still the capillary pressure increases in a broad region of the domain, as the gases further distribute. Figure 9 displays the evaluation of the three gas concentrations in gaseous phase and the gas saturation with higher resolved graphics for various time points to allow for better insight into the effects of phase change.
We observe that as long as gas is injected, the gas front and the interface between gas phase presence and absence effectively coincides. After when the gas injection is switched off after half of the simulated time, the total amount of gas reduces again due to the outflow at the zero Dirichlet conditions on the right boundary. In particular, the gas migrates to the northern, i.e. upper, boundary, where there are no flux conditions due to the cap rock and accumulates there. The gas saturation vanishes step by step at lower regions, the interface between presence of the gas phase and absence of a gas phase still remains relatively sharp, and migrates from bottom to top. Whereas at injection phase, the gas saturation has its maximum at the injection site, and more and more at the upper cap rock boundary, the maximum completely locates directly below the upper boundary, i.e. the cap rock. The relatively sharp interface region between the presence of the gas phase and the absence of a gas phase is due to the advection dominated transport processes. The migration and accumulation below the impermeable cap rock is due the gravitation, which causes that the gases in gaseous phase (which are in equilibrium with their liquid counterpart) migrate into the upper direction. (Note that the gas in gaseous phase does not "disappear" into the liquid phase due to the equilibrium conditions.) Concluding, we observe important phase change effects.

Injection of various gases and mineral trapping
The example we present at last demonstrates that our method is capable to simulate very efficiently one of the most promising technique of reduction of climate killing gas in the atmosphere, the injection of gas into deep layers causing mineral trapping of climate killing gases. As we can This is the realistic case of gas injection into deep layers, as usually CO 2 is not injected alone, but in combination with additional gases, among those CH 4 and SO 2 usually are the prominent ones besides CO 2 [56]. Therefore, as a proof of concept, the case we consider extend the case of Brunner and Knabner [6] in that we add two additional gases which are in equilibrium with their counterpart and which participate further, within the liquid phase, with other components. Their behavior is similar to CH 4 and H 2 S [56], therefore we note them with these names.
For these new gases, we assume Henry's law for simplicity, and for CO 2 we use the spline description given by Spycher-Pruess, as in the case of Brunner and Knabner for one gas [6]. The relations between mole fractions in  6 Henry's law: Relation between molar fractions in liquid phase and partial pressures in gas phase for the example with three dissolved gases in equilibrium and water without interaction, and restricted injection time liquid phase and partial pressures, i.e. the basis for the exchange laws, are displayed in Fig. 10.
As our study is directed to be a proof-of-concept study, and as in reality, CO 2 is that component with the highest injection rate among the gases, for the EOS we keep close to [6]: For the gas EOS, we assume the validity of the Duan Moeller law [16] for the relation of gas molar density and total gas pressure rather than partial gas pressure of CO 2 , and the validity of Garcia's law [20] for the relation of the molar density of the liquid phase and the molar fraction of dissolved CO 2 in liquid phase, i.e. for the liquid EOS. On this course, we make the mathematical complexity higher because we have complex splines, and thus establish a proof-of-concept.
Finally, the viscosity of the gas phase is described by the approach of Fenghour et al. [19], cf. also the paper of Knabner and Brunner [6].
Our simulations are already close to the real case of multiple gas injection into deep layers combined with reactive multiphase flow and therefore, our numerical techniques likely will be an advanced candidate for upcoming two phase reactive flow computations within realistic scenarios. Hence, our case extends elements of Brunner and Knabner [6] with heuristic elements. We take the chemistry of this study, but we add two additional gases. Most values which were used in the study of Brunner and Knabner [6] are still used. For the new parameters, we use in major part heuristic ones.
The parameters used in this simulation are listed in Tables 1 and 4.
The chemistry we use in this example is given by MinA + 3H + ←→ Me 3+ + SiO 2 (185) The blue lines in Eqs. 179-186 arise from the [6] case, the red ones are inspired by [56]. All 3 gases are injected into the porous domain (note that the benchmark idea of Sin et al. [56] does not contain gas injection).
Note: MinA is a Silicate, MinB a Carbonate, and Me could be e.g. Aluminium.
The reaction system consists of three minerals (Min A, Min B, and Calcite) and, including water, nine aqueous reacting species. We have three gases in interchange with their liquid counterpart. The first, second and third reaction represent the interphase mass exchange of CO 2 , H 2 S and CH 4 between their respective gas and liquid phase. The fourth and the sixth reaction allow a transition of dissolved CO 2 into HCO − 3 and calcite. They both affect and are affected by the pH, i. e., the concentration of H + .
This holds true as well for the case of the fifth reaction, the reaction of the new gas H 2 S with hydrogen and HS − . If a Silicate (Min A) is present at the initial state, it dissolves at high H + concentrations. This causes the release of metal ions by the seventh reaction. Finally, the metal ions can In Table 5, we compare the properties of this study to other ones popular in the field. Indeed, the study of [6] was able to simulate the injection of one gas in combination with mineral trapping, but not several gases. The study of [56] considered the case of various gases, but no mineral trapping and no gas injection. Our study combines all these effects in an efficient manner. Figures 11 and 12 display screenshots of the evaluation, we display the results at two time snap shots for the four- Like in the case of Brunner and Knabner [6] for one gas species, soon after the beginning of the injection of now various gases, the gas phase appears. The gase phase, similar to the case of Brunner, moves to the top where the impermeable boundary is located, and accumulates below this upper boundary. The mineral trapping of CO 2 happens in a similar way as in Brunners case, in that the dissolved CO 2 causes the appearance of H + , which in our case is strengthened by the dissolved H 2 S, which also causes the appearance of additional H + . So the pH decreases even more than in the case with only CO 2 . This induces the dissolution of Calcite and MinA. Metal ions and bicarbonate are released, and this causes the precipitation of MinB. All expected major effects are well reproduced by our simulations, such as the appearance of the gas phase, the precipitation and the dissolution of minerals.
Note that also in case of complete equilibrium, material can be transferred from one phase to another one, as for example precipitation of mineral reduces the amount of concentration in liquid phase, which again causes the immediate change from gaseous to liquid phase due to the equilibrium conditions. As well, the advective and diffusive transport of the gaseous and liquid phase does not happen with the exact same "speed", leading to concentrations which are not in equilibrium between gaseous and liquid phase in case of kinetic reactions. However, due to the AEs of the equilibrium conditions (Henry, Spycher-Pruess,. . . ), the equilibrium is imposed instantaneously again. This indicates phase change effects, which however are not directly "visible". Concluding, there are several processes which cause phase change effects, however it is difficult to "visualize" these effects directly.
In Fig. 13, we study the behavior of the computations along a trajectory and the respective grid convergence. Walking along a pathway defined by the diagonal from the lower left to the upper right, we observe that the results are quite stable with increasing grid level. We display the behavior of the gas saturation over grid level 1 to 4. We display a generalized relative difference R ij (w) along the path length of the trajectory w = x 2 + y 2 defined by (188) Figure 13 shows clearly that with increasing grid level, the results get more and more stable. Having such small differences after such a high number of physical time passed within the simulations demonstrates the excellent robustness of our algorithm.

Discussion of the results
Mineral storage of climate killing gases is a challenging topic in science presently. Mathematical formulations and simulations play a key tool to figure out the most efficient way to reduce the climate killing gases in the atmosphere. We demonstrated that our globally implicit solver, the PDE reduction method developed by Kräutle and Knabner [35] (PRM) applies very well to the very complex case of the injection of various climate killing gases into deep layers of rocks.
The geophysics of the models we simulate describe partially miscible multiphase multicomponent flow in porous media with general geochemical reactions and possibly appearing/disappearing gas phase, and general geochemical reactions, such as mineral dissolution and precipitation, i.e. nearly any kind of equilibrium and kinetic reactions. To evaluate the highly complex equations, we apply the globally implicit solver for multiphase multicomponent flow based on the refined reduction method developed by Kräutle and Knabner [35]. We presented the extension of the technique published in Brunner and Knabner [6] from one to various gas species which get injected into deep layers leading to mineral trapping.
In particular, we extended the concept of the incorporation of one gas within the framework of the PRM as published by Brunner and Knabner [6] to the case of an arbitrary number of gases. To do so, the numerical methods published in [6] needed to be extended substantially. Further, our numerical experiments showed substantial phase change effects in the case of various gases.
As from the very beginning of the application of the PRM [35], the transformation of the variables into reaction variant and invariant variables and the following subtraction procedure strongly reduce the size of the nonlinearly coupled global system. The division of the system into a local and global system efficiently treated but nevertheless coupled (not splitted!) by means of the resolution function further enables a fully global implicit solver with a huge number of advantages compared to classical approaches: At first, the global system contains all PDEs and some algebraic equations. However, due to the procedure we applied before, it contains only few nonlinearities, much less than other approaches such as the direct substition appreach DSA, and therefore, the structure of the Jacobian matrix is much simpler and in particular, much more sparse. This allows a comparably fast solution of the arising systems of linear equations which arise in the course of the application of the semismooth nested Newton solver. In our approach, we use a Finite Element discretization in space

Fig. 13
Grid convergence of gas saturation in Eq. 188 under spatial refinement along the diagonal from the lower left to the upper right compared at the final time step result at t = 6 · 10 6 s. The results are getting more and more stable with increasing grid level. The small peak in the relative difference between level 3 and level 4 is caused by the small shift of the interface, which is not surprising after such a long physical time simulated, namely as the small peak at the interface nearly already disappears at the highest grid level. Simulation result and trajectory displayed for grid at level 4. Note that the difference at the right boundary is related to the Dirichlet values there and the different size of the grid cells, so should not surprise at all where the advection parts are stabilized by means of Finite Volumes. For the discretization in time, we use the implicit Euler method. The local system iteratively decouples in very small subsystems at each node or basic element of the discretization, depending on the discretization type. This allows for the perfectly parallel computation of a substantial part of the system. Along the lines of Brunner and Knabner [6] (based on [48]), the case of the appearance and disappearance of phases can be handled without any discontinuous switch of variables.

Conclusions and outlook
We elaborated and presented the extension of the PRM [35] (which was extended to two-phase reactive transport for one gas in [6] already) to the case of an arbitrary number of gases in gaseous phase. We presented simulations where we studied phase change effects for the case of the injection of various gases into deep layers of porous media, and we presented results for the injection of various climate killing gases into deep layers in the context of mineral trapping scenarios. These gases are in exchange with their dissolved counterpart in liquid phase. We considered the case of a substantial number of liquid and solid components in the rock, and studied the efficient application of our technique to the case of the presented mineral trapping scenarios which happen in the course of the injection of these gases. Our computations are stable and robust. Such computations are at the high end of present technology.
As our study serves as a proof of concept, in part, we used heuristic parameters. Anyhow our approach likely is one of the very first successful application of globally implicit solvers for the case of the injection and reactive transport of multiple gas components into deep layers of porous media inducing mineral trapping and we are approaching highly realistic scenarios.
Indeed, up to now, to our knowledge, still there exists no benchmark for two-phase reactive flow scenarios in porous media which explore the effects of the injection of various gases into porous media. Our numerical techniques might be used for future computations in this direction. M.M.K. thanks Erlangen University for the kind funding, and the RRZE of Erlangen University for the supplied computing time at the Emmy cluster, and Thomas Zeiser and Markus Wittmann from the RRZE for their kind help.
Funding Open Access funding enabled and organized by Projekt DEAL.

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.