Upon injection into a reservoir, CO2 migrates and interacts with the host rock and pore water (Benson et al. 2005, 2012). Supercritical or gaseous CO2 dissolves in the pore aqueous solution as a function of the pressure, temperature, and salinity conditions. This process changes the local chemistry; in that it significantly decreases the pH and promotes geochemical reactions such as the dissolution of primary minerals and precipitation of secondary phases. These reactions then alter the porous network and can have a significant feedback on the migration of fluids (water and gas) in the reservoir. In addition, the dissolution of CO2 in the aqueous phase changes the overall composition of both fluids, thus impacting the gravitational and viscous forces, as well as the motion of both phases. The importance of multiphase flow and reactive transport is thus critical to carbon capture and storage (CCS), and moreover, it is necessary for CCS operators to consider these induced chemical and physical processes in predictive simulations, thus ensuring that CO2 does not leak into surrounding aquifers or the reservoir surface. Hence, accurate numerical multiphase flow and reactive transport modeling approaches are of paramount importance. Reactive transport simulations provide a unique means of quantitative and proper consideration of the coupling between the gas migration and its interactions with the pore water, in conjunction with the reactivity of the porous host rock (Garcia 2003; Ennis-King and Paterson 2005; Lu and Lichtner 2005; Nordbotten et al. 2005; Le Gallo et al. 2006; Trenty et al. 2006; Hassanzadeh et al. 2007; Pruess and Spycher 2007; Gaus et al. 2008; Nordbotten et al. 2008; Goerke et al. 2011; Sin 2015; Sin et al. 2017b). This chapter presents an overview of the current mathematical and numerical deterministic approaches to the modeling of such multiphase multicomponent reactive systems. Applications within the context of CO2 storage, in addition to several other examples, are used to illustrate the concepts employed, their limitations, and the evolution of the simulation methods developed in the recent decades.

The presence of a gaseous phase, in addition to an aqueous and solid phase, adds complexity to the physical and chemical processes that should be considered in a reactive transport approach. Reactive transport models were developed initially to model saturated systems, e.g., only under consideration of the liquid-solid interactions. Several methods were developed accordingly for the modeling of single-phase reactive transport. For problems with large numbers of components, complex geochemical reaction paths, and a high degree of coupling between the transport and chemistry, the set of equations results in large nonlinear systems that require efficient mathematical and numerical techniques (Steefel 2019, this volume). Significant research efforts have been directed toward the reduction of this computational complexity, especially with respect to the dimensions of the resulting systems. In this context, methods based on operator splitting between transport and reactions have typically been more advantageous for implementation due to their simplicity when compared with methods that collectively address the full set of equations (Yeh and Tripathi 1989, 1991; Xu and Pruess 1998; van der Lee et al. 2003; Molins et al. 2004; Kräutle and Knabner 2007; Molins and Mayer 2007; Steefel et al. 2014). The consideration of multiphase systems yields additional nonlinearities in flow- and geochemistry-related equations. The solution of the multiphase multicomponent flow and reactive transport (MMF&RT) increases in complexity in accordance with the evolution of the reactive system. The vapor–liquid equilibrium of multicomponent systems requires sophisticated models to handle vanishing phases and the related degeneration of equations that create numerical instabilities. Accurate modeling of the thermodynamic state for multicomponent fluids is therefore necessary for both the geochemistry and the flow, as it contributes to the phase equilibrium (transfers among phases and associated mass balance), in addition to the fluid and transport properties. Furthermore, it can have phenomenological effects, e.g., the changing of the flow regime, fluid dynamics, and matrix structure. Coupling methods for thermodynamic, hydrodynamic, and chemical modeling are therefore required, and they should ensure numerical stability and robustness.

In the modeling of multiphase systems, numerical issues can come both from the chemical heterogeneity and from the resolution of the set of nonlinear flow equations. The complexity of MMF&RT modeling can thus be attributed to both the complex physics and the high-dimensional nature of these problems. Moreover, several critical phenomena can be reproduced only by high-resolution simulations. Over the last four decades, numerous approaches have been developed and tested with global implicit or sequential iterative numerical methods (Peszynska and Sun 2002; Nghiem et al. 2004; Hao et al. 2012; Fan et al. 2012; Wei 2012; Wheeler et al. 2012; Xu et al. 2012; Nardi et al. 2014; Steefel et al. 2014; Hron et al. 2015; Parkhurst and Wissmeier 2015; Ahusborde and El Ossmani 2017; Sin et al. 2017a; Brunner and Knabner 2019). A range of multiphase multicomponent flow and transport solvers (e.g., White et al. 2012; Xu et al. 2014; Sin et al. 2017a; Voskov 2017) and geochemical solvers (e.g., Appelo et al. 2014; Corvisier et al. 2014; Leal et al. 2014) were recently developed to model CO2 geological sequestration systems, among other applications. In this paper, a detailed state-of-the-art review focused on coupling issues is presented, and several recommendations are provided for future developments.

This chapter is divided into three main parts: (1) governing processes in multiphase multicomponent flow and reactive transport (MMF&RT), (2) mathematical and numerical approaches of MMF&RT modeling, and (3) applications with respect to CO2 storage. Thereafter, discussion of other applications that illustrate the ability of existing methods to represent general non-ideal multicomponent gaseous systems within a wide range of temperature and pressure conditions, in addition to their possible interactions with water/brine/rock, are presented.


The governing processes and present issues of multiphase flow and geochemical modeling are introduced in this section. The definition of wettability is first presented, which is difficult for several systems; followed by the definitions of capillary pressure, residual and capillary trapping, and heterogeneity. In addition, the relative permeability is introduced, which is a characteristic property of multiphase flow. Recently, novel experimental and numerical technologies have provided further insight into different concepts such as residual non-wetting saturation, maximum relative permeability, governing forces, hysteresis, the impact of heterogeneity, and upscaling from pore- and core-scales. The different regimes of gas diffusion are then discussed in conjunction with corresponding modeling methods. Thereafter, basic formulations of multiphase flow are presented, from immiscible to compositional flows. After a brief presentation of geochemical modeling, e.g., mass balance and phase equilibria, a discussion of issues with respect to the modeling of multicomponent systems, such as the prediction of fugacity and activity coefficients, is given. Furthermore, thermodynamic properties such as density and viscosity are dependent on the composition, pressure, and temperature; and the challenge of describing appropriate equations of state (EOS) for multicomponent systems is significant. Hence, several methods used to predict these properties are discussed.

Capillary pressure and wettability

Porosity ϕ quantifies the space occupied by fluids or a void within a representative elementary volume (REV). When two fluids are separated by a distinct interface, these fluids are defined as immiscible, and the interfacial tension (IFT) results in a curved fluid–fluid surface within each pore. The cohesive and adhesive forces form a contact angle θ between the fluid–fluid interface and the solid face. If the contact angle θ is acute, the fluid is defined as wetting; otherwise, it is defined as non-wetting. For example, a gas and a liquid both in contact with the solid surfaces of the porous matrix create a gas–liquid interface with a curvature that results from balanced interfacial tensions.

In a system in which the number of fluids (Nf) is greater than one, each fluid represents a phase, and the system is referred to as multiphase. The saturation Sa is the ratio of the volume occupied by the phase a as a gas or liquid in a REV relative to the volume of pore space in this REV. Under the assumption that fluids occupy the entire pore space, the following relationship is satisfied:


The gas and liquid fluids in contact with the solid phase apply different pressures on the curved interface that separates them. This variation in pressure from one side to the other is defined as the capillary pressure pc:


where pg and pl are the gas and liquid pressures, respectively. The capillary pressure at a point in the interface can be evaluated with the mean radius of curvature and a local interfacial tension g using the Laplace equation, which takes the form of the Young–Laplace equation for a capillary tube of radius r:


When a gas enters into a water-saturated porous medium, a certain pressure is required to displace the wetting fluid that occupies the porous space. This desaturation process is referred to as drainage. Displacement of the wetting fluid by the non-wetting fluid is described as the imbibition process (Fig. 1a). The relationship that links the capillary pressure and the liquid saturation is dependent on the direction of the displacement (e.g., drainage or imbibition) and thus hysteresis effects are commonly exhibited.

During the passage of a wetting phase (imbibition), a portion of non-wetting fluid can be trapped in the pore space. This capillary trapping process is highly dependent on the contact angle. Capillary pressure and wettability have been extensively investigated for oil–water systems (Leverett 1939; Mungan 1966; Batychy and McCaffery 1978). The contact angle is dependent on the mineralogy and surface roughness, and it is difficult to estimate. However, Armstrong et al. (2012) demonstrated that the contact angles of oil–water systems can be directly measured using X-ray computed tomography (CT). For the supercritical CO2–water system (scCO2–water), most of the recent core-flooding experiments confirm that scCO2 remains a non-wetting phase during water displacement (Krevor et al. 2012; Pentland et al. 2011; Pini et al. 2012; Andrew et al. 2014; Øren et al. 2019). Pini et al. (2012) observed contact angles of 20–45° at 25 °C and 0–5° at 50 °C for a CO2–water system by fitting the contact angle for the IFT values obtained from Berea sandstone. Espinoza and Santamarina (2010) reported a contact angle of 30° for CO2–brine–calcite at 0–10 MPa and 23.35 °C. Broseta et al. (2012) observed contact angles of 60–75° during imbibition at 0.5–14 MPa and 35 °C. Andrew et al. (2014) confirmed that scCO2–brine–carbonate is a water-wet system at 10 MPa and 50 °C by measuring contact angles of 35–55° using micro-CT. Finally, Øren et al. (2019) observed angle ranges 20–30° and 30–40° for CO2–brine analogue octane:brine within reservoir samples from the Paaratte sandstone formation in the Otway Basin, Australia, using a quasi-static pore-scale model based on in-situ three-dimensional (3D) pore-scale images.

Brooks and Corey (1964) measured the capillary pressure head (pl / ρlg, where ρl is the liquid density and g is the gravity constant) for different types of soil. Similarities between the retention curves were observed, and the capillary pressure was found to be dependent on the effective saturation Sl, which can be expressed as follows:


where pe is the entry pressure (Fig. 1a); which is the pressure required to initiate a displacement of the wetting phase (e.g., water) by the non-wetting phase (e.g., air, oil, and CO2). The effective saturation is defined by


where Slr is the residual or irreducible saturation of the wetting phase. The Brooks–Corey model for capillary pressure was designed for drainage processes, where λ is a fitting parameter that varies with respect to the soil texture, e.g., the pore size distribution, within the range 0.5–7.5 (Brooks and Corey 1964). Experimental studies conducted on Berea sandstone yielded λ in a range of 0.64–2.4 (Krevor et al. 2012; Pini and Benson 2013). Pini and Benson (2013) carried out core flooding experiments using different fluids: gN2, gCO2, and scCO2 at 50 °C and 2.4 MPa, 2.4 MPa, and 9 MPa, respectively (where “g” represents gaseous, and “sc” represents supercritical). The entry pressure increased in accordance with an increase in the interfacial tension: 4.94 kPa (41 mN/m), 6.8 kPa (57 mN/m), and 8.13 kPa (65 mN/m) for scCO2, gCO2, and gN2, respectively.

The trapping process was formulated by Land (1968) relating the initial and residual saturations as follows:


where Sgr is the residual (trapped) saturation of the non-wetting phase; Sg,init is the initial or turning point saturation; and C is the Land coefficient, which ranges from 0.9–2.5 for Paaratte sandstone (Øren et al. 2019). The Land coefficient has been found to range from 2–7 for Berea sandstone, and from 2–5 for other sandstone and carbonate rocks (Krevor et al. 2015). The collected data revealed that the trapped CO2 saturation minimum is roughly 0.1, and most rocks can immobilize 30%–40% of CO2. The Land model (Eqn. 6) is the most widely used, although there are other trapping models (Carlson 1981; Spiteri et al. 2008).

The Brooks–Corey model (Eqn. 4) commonly describes drainage processes. In contrast, multiple methods were developed for the imbibition process to capture the hysteresis behavior and capillary trapping depending on the initial saturation (Killough 1976; Ruprecht et al. 2014; Pini and Benson 2017). Pini and Benson (2017) proposed a novel imbibition model based on effective mobile non-wetting saturation Sl,mob, such that Sg,mob + Sl,mob = 1:

pc,imb(Sl)=pc,init(S¯l,mob11λ)for S¯l,initS¯l,mob1

where Sl,mob is the effective mobile wetting phase saturation, and Sl,init is the turning-point saturation. Moreover, pc,init can be expressed based on the equality of imbibition and drainage capillary pressures at the turning point: pc,imb = pc of a drainage curve (Eqn. 4). This Brooks–Corey based formulation (Eqns. 4 and 7) was used to analyze the results of capillary pressure heterogeneity and the hysteresis effect for a scCO2–water system on Berea sandstone at the core and sub-core scales. A method of translation of physical properties and their heterogeneity for further numerical simulators was suggested by Pini and Benson (2017) using a capillary scaling method. Such studies (Pini et al. 2012; Pini and Benson 2013, 2017; Kuo and Benson 2015; Hejazi et al. 2019) describe a methodology that combines experiments and modeling for the analysis and mapping of capillary heterogeneity at the core and sub-core scales on sandstone and carbonate rocks.

The heterogeneity of rock samples is related to the grain shape and tortuosity, e.g., size and distribution of pore throats and bodies. This is characterized mostly by heterogeneous permeability K at a macroscopic scale. For multiphase flow, capillary pressure and relative permeability are major indicators of heterogeneity. The experimental data of capillary pressure are typically scaled by the Leverett J-function (Leverett 1941) J(S)=K/ϕpc(S)/γ; thus, they can be compared for similar textural properties (Pentland et al. 2011; Pini and Benson 2017). Using micro-CT imaging and numerical simulations, Andrew et al. (2014), Pini and Benson (2017), and Øren et al. (2019) demonstrated that the spatial distribution of heterogeneity at the micro-scale has a significant influence on the capillary pressure and relative permeability at the core scale. Recent experimental results for nine sandstone samples (Ni et al. 2019) confirmed that CO2 trapping decreases in accordance with an increase in porosity, and increases in accordance with an increase in the degree of heterogeneity. Jackson and Krevor (2019) demonstrated that experimental observations of capillary heterogeneity, and therefore relative permeability, which includes the capillary heterogeneity at the core scale, can be translated to the meter scale. Thus core-averaged methods can be used to estimate the capillary pressure, relative permeability during drainage and imbibition, and residual trapping.

Several analytical and numerical studies have emphasized the effects of capillary trapping on CO2 storage, the plume extent, and fate (Doughty 2007; Hesse et al. 2008; Juanes et al. 2010; Nordbotten and Dahle 2011; Golding et al. 2013). The heterogeneity of porous media has a significant influence on the CO2 dispersal and trapping due to local impermeable barriers that cannot be captured using simple upscaling techniques of vertical permeability reduction (Green and Ennis-King 2010; Hesse and Woods 2010). Trevisan et al. (2015) carried out reservoir scale experiments to investigate the impact of capillary heterogeneity on flow dynamics during drainage and imbibition. Moreover, heterogeneity, which is a variation of effective permeability, may yield a smaller onset time of convection than that in an equivalent homogeneous system (Green and Ennis-King 2010). Mouche et al. (2010) proposed a method for upscaling the vertical migration of CO2 in a layered porous medium at capillary-dominant conditions. Rabinovich et al. (2015) developed a method for upscaling capillary pressure and rock relative permeability curves from the core to coarse scale, and from the sub-core to core scale. This approach allows for the modeling of drainage processes that consider high degrees of capillary heterogeneity at the fine scale, and do not consider buoyancy/gravity and the hysteresis effect. The modeling of heterogeneous media at the reservoir scale, in addition to obtaining and integrating information on the characteristic functions of different rocks at the core, sub-core, and pore scales, is thus necessary; although challenging.

Relative permeability

Darcy's law extended to multiphase flow is formulated as uα = − krαKα(∇pα − ραg) where the effective permeability krαK is a permeability tensor of phase α, and μα and ρα are the viscosity and density of phase a, respectively. The relative permeability krα is typically assumed to be a non-directional fluid–matrix property, although it can be a tensor. By construction, the relative permeability is dependent on the pore structure and saturation. Typical examples of non-wetting and wetting relative permeability curves are presented in Figure 1b. The relative permeability of the wetting phase decreases rapidly with the progress of drainage. The wetting fluid is immobile at significantly high saturations Slr (Fig. 1b), whereas large pores are occupied by the non-wetting phase. This behavior is a characteristic of strongly water-wet systems. The relative permeability krg increases in accordance with an increase in saturation Sg up to its maximum value krg,max at the maximum non-wetting phase saturation, e.g., the irreducible wetting phase saturation Sl,irr. These endpoint values are critical for rock characterization, e.g., the estimation of the storage capacity, residual trapping, and plume evolution in the CCS context.

Extensive measurements of relative permeability were carried out for CO2–brine and H2S–brine systems on sandstone, carbonate, shale, and anhydrite rocks from Alberta (Bachu and Bennion 2008; Bachu 2013), which yielded very low krg,maxvalues of drainage curves (approximately 0.1–0.6). Different values were obtained for scCO2–brine/Berea by Krevor et al. (2012) and Pini and Benson (2013): Sl,irr of approximately 0.45 and 0.35, and krg,max of 0.38 and 0.9, respectively. Reported discrepancies between endpoint values can be attributed to the different measurement methods employed, in addition to the different theoretical assumptions of averaging. Krevor et al. (2012) suggested that endpoint values of gas saturation and permeability can be measured when sufficiently high capillary pressures are reached. This is attributed to the low viscosity of CO2 in comparison with the viscosity of water (viscosity ratio μwater / μCO2 is approximately 24–30 under high temperature and pressure conditions). Based on Darcy's law, changes in the capillary pressure are governed by changes in the viscous pressure during the experiments, as demonstrated by Pini and Benson (2013) using (steady-state) scCO2, gCO2, and gN2 flooding experiments within a wide range of pressure, temperature, and salinity conditions. High capillary pressure conditions were obtained with large flooding rates. High values of relative permeabilities krg,max and capillary pressures were simultaneously measured at the inlet face of the core. Moreover, the measurements of the saturation obtained using the X-ray CT data processing were directly related to the measured capillary pressure and relative permeability. The required high capillary pressure and high flow rate conditions may be difficult to achieve due to limited experimental capabilities (Krevor et al. 2012). This steady-state technique combined with the X-ray CT scanning described by Ramakrishnan and Cappiello (1991) and Pini and Benson (2013) was applied for measurements of the characteristic properties and allowed for the quantification of the endpoint behavior.

Sensitivity analysis experiments conducted on scCO2, gCO2, and gN2–brine systems (Pini and Benson 2013) indicated that the IFT has a negligible influence on the relative permeability over the range of 40–65 mN/m, which corresponds to a wide range of reservoir pressure (P) and temperature (T) conditions. Consequently, analogue fluids can be used for the characterization of relative permeability. Al-Menhali et al. (2015) conducted a set of experiments to characterize the capillary pressure and relative permeability in CO2–brine and N2–water/Berea systems within a wide range of reservoir conditions (P = 5–20 MPa, T = 25–50 °C, and 0–5 mol/kg of NaCl). A slight difference between the N2–brine and CO2–brine relative permeability curves was observed, which may be due to small sample heterogeneities and other measurement uncertainties. However, in the case of the relative permeabilities measured for different gas viscosities at a fixed IFT, significant differences between the permeability curves under different viscosity conditions were demonstrated (Reynolds and Krevor 2015). Most importantly, the observed changes in the relative permeability were due to a shift in the balance between the capillary-dominated and viscous-dominated conditions, and not variations in the fluid properties.

The capillary number Ncapillary, which is a ratio between the capillary to viscous forces, can be used to characterize the behavior of a system. There are different definitions of Ncapillary, depending on the targeted application. For example, Virnovsky et al. (2004) proposed Ncapillary = HΔp / (L|Δpc(fw)|), where H and L are longitudinal and transverse length scales, respectively; the pressure drop Δp represents the viscous forces; and |Δpc(fw)| is the magnitude of characteristic difference in the capillary pressure. Reynolds and Krevor (2015) demonstrated that if the system approaches a viscous limit with increasing capillary number Ncapillary (by increasing viscous pressure drop Δp), the measured relative permeability is actually the intrinsic relative permeability of the rock, and the saturation of the non-wetting phase is then insensitive to capillary heterogeneity and governed by the absolute permeability K. On the contrary, when the system is under capillary-limited conditions, e.g., Ncapillary ≪ 1, the effective relative permeability is measured, which is highly dependent on capillary heterogeneity, and therefore on the flow rate and fluid properties. Most relative permeability curves (Krevor et al. 2012; Pini and Benson 2013; Reynolds and Krevor 2015) are measured at the capillary limit and are then sensitive to capillary heterogeneity and reservoir conditions. Reynolds and Krevor (2015) proposed the use of scaled relative permeability for the characterization of heterogeneous cores, e.g., the measurement of the intrinsic and effective permeabilities of representative rock samples at the required reservoir conditions to better describe flow at the reservoir scale. The importance of sample heterogeneity was also highlighted by Zhang et al. (2013, 2017), who measured the relative permeabilities of Berea samples with bedding and lamination and of a heterogeneous core sample with a low permeability. Rabinovich et al. (2015) presented an upscaling technique of capillary pressure curves for heterogeneous samples at the capillary limit, not considering gravity and the hysteresis effect. Kuo and Benson (2015) proposed a two-dimensional (2D) semi-analytical method to predict the average saturation for 3D two-phase flow models of heterogeneous cores. A recent upscaling approach (Jackson et al. 2018; Jackson and Krevor 2019) was designed for the characterization of heterogeneous cores using combined experimental data and numerical modeling. The heterogeneity of a Bentheimer sandstone sample and that of a Bunter sandstone sample were defined with the viscous limit relative permeability and scaled capillary pressure curves. The characteristic properties were estimated within a wide range of Ncapillary values that can be used for flow simulations. In addition, the hysteretic curves were also obtained for a Bunter sandstone core.

Most of the experimental relative permeability data, especially the drainage curves, are generally fitted to the Brooks–Corey model (Brooks and Corey 1964; Krevor et al. 2012; Pini and Benson 2013; Reynolds and Krevor 2015; Jackson et al. 2018; Jackson and Krevor 2019):


The relationships presented above correlate the permeabilities to the capillary pressure (Eqn. 4). Although data can easily be fitted to such a simple power model, relationships such as Equations (8) and (9) demonstrate the corresponding behavior of a capillary-controlled system. On the contrary when the capillary forces are absent, viscous forces govern the system, which leads to a uniform displacement of fluid. With respect to relative permeability curves, this is indicated by two diagonal lines. The models and frameworks commonly used for the construction of relative permeability curves and also of capillary pressure as well (Brooks and Corey 1964; Burdine 1953; van Genuchten 1980; Mualem 1976) can be generalized accordingly (Dury et al. 1999) as follows:


where Ki, i = 1, 2, 3 are fitting parameters. The terms on the right-hand side raised to the power K3 in Equations (10) and (11) are correlated with the pore connectivity and tortuosity; the terms on the right-hand side raised to the power K2 correspond to the capillarity effect. The same forms (Eqns. 10 and 11) can be employed for hysteretic relative permeability curves.

Several experimental and numerical studies have highlighted the hysteretic behavior of CO2–brine systems, the shape of these curves, and residual trapping. An important difference between the drainage and imbibition processes was observed (Bachu 2013; Bachu and Bennion 2008; Ruprecht et al. 2014; Krevor et al. 2015; Jackson and Krevor 2019; Øren et al. 2019). Ruprecht et al. (2014) analyzed three core flooding cycles of alternated drainage and imbibition processes on a Berea sandstone and found that residual trapping and the CO2 relative permeability significantly vary during the cycles. However, the water relative permeability did not exhibit hysteretic behavior. Øren et al. (2019) confirmed the hysteresis of the non-wetting relative permeability of a Paaratte sandstone using two different pore-scale modeling approaches. Interestingly, the observed non-wetting relative permeability curve of the imbibition process was higher than that of the drainage process, whereas it is typically opposite. This was attributed to the low aspect ratio of the Paaratte sample. The Killough model was employed to construct scanning curves for the capillary pressure and relative permeability, and the residual saturation was calculated by the Land equation (Eqn. 6). Moreover, Paterson et al. (2013) successfully carried out a test on residual trapping at the field scale, at the CO2CRC Otway site, Australia, during which 150 Mg of CO2 was injected; followed by the injection of 454 Mg of formation water, which forced residual and dissolution trapping. Residual saturation values of about 0.15–0.23 were estimated using different observation techniques, which confirmed the importance of residual and dissolution trapping mechanisms. Despite these recent studies, our current knowledge of the behavior of hysteretic curves, scanning curves, and hysteresis of residual trapping is still limiting, and a comprehensive database for cyclic core flooding on representative rock samples is thus required.


Fick's law, which is expressed as Jα,k = −ραDαXα,k, is generally used to describe molecular diffusion in liquid and gas phases; where Xa,k denotes the mass fraction of species k in fluid phase a, and Dα is the diffusion tensor. Fick's law can also be expressed in terms of the concentration or mole fraction. The molecular diffusion Dα can be formulated for each phase by an effective diffusion coefficient:


where τ is the tortuosity, which is generally represented by the Millington and Quirk model (Millington and Quirk 1961); and Dα is the molecular diffusion coefficient of fluid phase α. In this form, the diffusion coefficients in the aqueous phase are an average species-independent parameter required to maintain the electroneutrality of the solution (Tournassat and Steefel 2019, this volume). The molecular diffusion coefficient for the gas can be species-dependent, given that the gas species are neutral. Diffusion in the gas phase can be more complex than the molecular diffusion described by Fick's law. In particular, it may include Knudsen and nonequimolar diffusion. The Knudsen regime or free molecular flow occurs when the gas mean free path λ is significantly greater than the pore radius λp: Kn = λ / λp > 10. The Knudsen diffusion is dependent on the molecular weight, gas temperature, and pore radius. The kinetic theory of gases is then applied, and the dusty gas model of multicomponent systems can be employed (Cunningham and Williams 1980; Mason and Malinauskas 1983). Nonequimolar diffusion stems from a difference in molecular weight: lighter molecules are characterized by more rapid diffusion than heavier molecules. However, these diffusive slip fluxes do not separate gases.

For small Knudsen numbers, 10−3 < Kn < 10, when the mean free path of gas molecules is approximately the same as the pore radius or less, transition or viscous slip flow occurs. Due to the slip effect, which is referred to as the Klinkenberg effect, Darcy's law typically underestimates gas flow. Klinkenberg (1941) proposed several models of gas slip permeability, which correct the gas permeability with respect to the mean pressure and a fitting parameter of slippage. The Klinkenberg first-order correction of permeability evolves into a modified gas velocity, as follows:


where bK is the Klinkenberg first-order parameter. Gas flow in tight porous media can differ significantly from liquid flow, given that the effective permeability is pressure dependent and varies with respect to the gas composition. There are several forms of the Klinkenberg parameter. In particular, it can be presented with respect to the kinetic viscosity of the fluid, temperature, and molecular weight of the gas (Ertekin et al. 1986). Simplified forms with empirical parameters have been used to study the apparent permeability of gases such as hydrogen, helium, nitrogen, air, argon, and carbon dioxide in low permeable sands, sedimentary rocks, and shale formations (Heid et al. 1950; Florence et al. 2007; Ghanizadeh et al. 2014). Civan (2010) derived a correlation between the apparent gas permeability, Klinkenberg slippage parameter, and tortuosity. The apparent gas permeability can be greater than the initial intrinsic permeability by a factor of 103. Therefore, the Klinkenberg slippage effect is typically considered for the modeling of CO2 sequestration in unconventional reservoirs such as coal bed methane (Liu et al. 2015; Fan et al. 2019) and shale gas reservoirs (Sun et al. 2013).

General formulations of multiphase flow

The formulations of multiphase flow briefly presented below are widely used for the modeling of hydrogeochemical and petroleum engineering problems. Mathematical and numerical methods of immiscible flow have been extensively investigated by the petroleum industry (Craig 1971; Barenblatt et al. 1972; Bear 1972; Peaceman 1977; Aziz and Settari 1979; Hassanizadeh and Gray 1979; Chavent and Jaffré 1986; de Marsily 1986; Whitaker 1986; Lake 1989; Chen et al. 2006). Several major formulations of immiscible flow are discussed in this section, as they are employed for compositional flow. The optional unsaturated flow is also derived, as it is applicable to a wide range of other environmental problems.

Immiscible multiphase flow.

Mass conservation equations are required for each phase in a system of Nf fluid phases:


where mα = ϕSαρα is the mass of a fluid phase α, α∈ {0,..., Nf}, Fα = ραuαis the mass flux, and qa is the source or sink term. The phase pressures pα are related by the capillary pressure relationships:


There are various formulations of the set of equations, and their structures modify the coupling between equations. Several examples are presented below, starting from a general formulation with respect to pressure and saturation.

The pressure-saturation formulation is derived from the substitution of the Darcy velocities in the mass conservation Equation (14):


For each phase α, a natural variable is selected, namely, Sα or pα resulting in Nf nonlinear equations with Nf unknowns. The equations are strongly coupled through characteristic functions pc and krα. Primary variables should be carefully selected, as this leads to specific limitations. The modeling of a two-phase system allows for pairs (pl, Sg) and (pg, Sl). If (pl, Sg) are primary variables, the conservation equation of the non-wetting phase contains a term krgpc with singularities at the endpoints. The application of the Brooks–Corey (BC) or the van Genuchten (vG) model yields the following:


When selecting (pg, Sl) as primary variables, the term krlpc has singularities at Sl = 1.


After the linearization of the resulting system of equations written in discrete residual form, numerical problems at the saturation endpoints then appear in the Jacobian matrix (Sin 2015). The pressure–pressure formulation can be used to prevent restrictions on saturation range, as discussed below. Alternatively, this issue can be solved by employing numerical methods, e.g., the use of spline curves.

The decoupled pressure and saturation formulations separate the pressure and saturation equations. Summarizing the conservation equations, the temporal derivative of saturation can be eliminated:


where the total velocity u=Σαuα. The Darcy velocities for two phases are defined as follows:


where λα = krαα is the mobility, and fα=λα/Σβλβ is the fractional flow.

Equation (19) is solved by considering a phase pressure or global pressure as a primary variable. If a phase pressure is selected, the saturation equation is then hyperbolic, and coupling of the system is strong. If the global pressure is a primary variable, the saturation equation is weakly coupled to the pressure equation.

When a phase pressure ρα is selected as a primary variable, Equation (16) is formulated for the Darcy velocity uα using Equations (20) or (21). This yields a pure hyperbolic equation for saturation. When using a phase pressure, the pressure equation contains a capillary term fαpc. The result is that the phase pressure and saturation are decoupled.

Decoupled global pressure and saturation is also possible when the total velocity is given by u=λtotK(pΣαfαραg), where the total pressure p=Σαfαpα. The phase velocities can then be expressed with respect to the total velocity and substituted into one of the mass balance equations. Hence, the pressure and saturation equations are weakly coupled through the mobilities and fractional flow. However, boundary and initial conditions have no physical meaning, and the pre- and post-treatment of variables is then required. These techniques for the immiscible multiphase flow allow for strong or weak mathematical decoupling, which can be employed for compositional flow.

Unsaturated flow.

The modeling of unsaturated flow is critical in soil science, environmental engineering, and groundwater hydrology, e.g., to predict fluid motion in the vadose zone, changes in groundwater systems due to root water uptake, and the evaporation of water from soil. The vapor movement is significantly less important in soil science. The air pressure above the water table corresponds to the atmospheric pressure. Water is under tension and attracted upward by capillary forces. The total hydraulic head htotal consists of the pressure head h and elevation head z. When the pressure head is negative, water molecules are under tension, which indicates the capillary fringe. Air motion is neglected, which is an important hypothesis and differs from the multiphase flow formulation above: pair = 0, pw = −pc. This yields the following:


where ψ = pc / ρw g = pc / γw is the capillary pressure head.

The Darcy velocity can be expressed as follows:


where Kcond(h) = Kγw / μw is the hydraulic conductivity, which is related to the intrinsic permeability by the gravitational term γw and the water viscosity μw. Obviously, Kcond(h) can be modeled using the Brooks–Corey model and other empirical models (Eqns. 811). Consequently, the Darcy velocity can be formulated with respect to the moisture θ = ϕSw:


where Dc(θ) = Kcond(θ)/C(θ) is the capillary diffusivity, and C(h) = dθ/dh is the moisture capacity. The water balance in the vadose zone is described by Richards equation, as follows:


The equation above can be expressed in three different forms using Equation (24): mixed, h-based, and θ-based. The h-based formulation is not conservative, and the discretization of the accumulation term C∂h/∂t can lead to numerical errors. Although the θ-based equation is conservative, it fails in saturated conditions. Therefore, the mixed formulation is generally advised. Given that Equation (25) can be coupled with multicomponent solute transport, it covers numerous applications in the vadose zone (Suarez and Simunek 1996). However, assumptions such as inert gasses and the incompressibility of fluids are limiting for a wide range of problems that involve gas dynamics, the heterogeneous composition of fluids, and important changes in their composition due to physical and chemical processes.

Compositional flow.

When considering miscible fluids (i.e., no sharp interface between fluids and phase changes), the compositional flow should be solved. For a system of Nc components, the mass conservation equations and closure relationship for mass fractions can be expressed as follows:


The conservation equations can also be expressed with respect to mole fractions or concentrations.

The compositional flow problem consists of Nc equations and 2 Nf constitutive relationships (Eqns. 1, 15, 26, and 27), and it describes a general multiphase multicomponent flow (MMF) with 2 Nf + NcNf unknowns: pα, Sα, and Xα,k. The missing Nc (Nf − 1) equations correspond to the phase equilibrium relationship detailed below. Numerical and mathematical methods for MMF are presented in the section on multiphase multicomponent flow (MMF) approaches, followed by the coupling methods with reactive transport or geochemical reactions.

Geochemical reactions

In most geochemical or reactive transport codes, reservoirs are considered infinite with respect to gases, which is equivalent to the fixing of the corresponding dissolved concentrations in the aqueous solution. This assumption holds when modeling open or semi-open reactors, into which gases are continuously pumped to maintain the global pressure constant, or shallow aquifers in contact with the atmosphere. In particular, the mass balances that involve gaseous species are not considered, and the dissolved aqueous species just obey the corresponding mass action law. The system of equations required to be solved is then significantly simpler. Nevertheless, for closed reactors or for the general two-phase reactive-transport problem, this approach is not satisfactory, and codes intended for these purposes have been subsequently improved. The focus of this chapter is on the CO2 geological storage application; hence, several aspects are therefore highlighted in this context. Many numerical studies were carried out using dissolved impurities and an artificial gas phase (Xu et al. 2007; Xiao et al. 2009; André et al. 2012). When CO2 is dissolved, CO2–water–rock interactions can be handled with a classical geochemical solver. However, to process changes in the pressure, temperature, and salinity conditions, accurate models are required to predict the gas–water equilibrium; especially when gasses other than CO2 are considered. In addition, recent developments of gas–water–rock interactions allow for several other contexts (e.g., acidic–basic/oxidative–reducing conditions, high–low pressure/temperature/salinity, and kinetically-controlled biotic reactions involving gas); thus leading to various applications (natural gas reservoirs, hydrogen storage in salt caverns, enhance coal-bed methane, and soil gas migration).

Chemical equilibrium and mass balances.

Geochemical solvers or reactive transport codes can manage chemical equilibrium calculations in two manners: based on a stoichiometric approach (PHREEQC, Parkhurst and Appelo 2013; CHESS, van der Lee 2009) or on a nonstoichiometric formulation (GEM-Selektor, Kulik et al. 2004).

On one hand, the stoichiometric approach requires the expression of all the chemical reactions and associated mass actions laws with respect to activities. On the other hand, the nonstoichiometric approach is based on the minimization of the total Gibbs energy of the considered systems, which is calculated from the chemical potential of all the species. In both cases, solvers also include mass or mole balance equations that are expressed with respect to the elements or total components, which are not presented here. They also generally comprise mineral precipitation/dissolution, solid solutions, kinetically controlled reactions, sorption, and isotopes.

Water–gas equilibrium.

Many models have been developed since the early 1990s to predict the water–gas equilibrium; and in particular, mutual solubilities using various strategies. Studies are largely focused on hydrocarbons and CO2, due to applications for the oil and gas industries (Michelsen and Mollerup 2007); or more recently, for the greenhouse gas geological storage and due to the quantity of experimental data available to fit or test models. Søreide and Whitson (1992) adapted the Peng–Robinson EOS (Peng and Robinson 1976); Duan and Sun (2003) derived a virial EOS; and Spycher et al. (2003) employed a cubic EOS (Redlich and Kwong 1949). Various geochemical solvers were recently developed with specific modules to incorporate the cubic EOS to address non-ideal gas phases at high pressures and temperatures (Appelo et al. 2014; Corvisier et al. 2014; Leal et al. 2014). The phase equilibrium is presented below.

The local equilibrium of a system requires a minimum of the Gibbs free energy at a constant temperature T, pressure p, and composition. However, the system could be monophasic or composed of multiple phases. In geochemical modeling, it is common to use the heterogeneous approach: an activity model describes the activity coefficient in the liquid phase, whereas the fugacity coefficient in the vapor is modeled by an equation of state. Given that not only pure water is considered but also various ionic species, the asymmetric convention (φ−γ) is used.

When considering non-stoichiometric geochemical solvers, the activities of the aqueous species (dissolved gases) and fugacities of the gases appear in the chemical potential functions:


where Ci is a concentration quantity (e.g., molality), γi is an activity coefficient, yi is the gaseous molar fraction, and φig is a fugacity coefficient. As stated by Gibbs, at local equilibrium, the chemical potentials of each component are equal in the multiphase system at constant temperature, pressure, and composition. This is a criterion, and the solution of the equilibrium conditions must then be coupled with a stability criterion. Lewis and Randall (1923) defined the fugacity of a component as follows:

fi(T,p,x)=pexp(gi(T,p,x)gi,purePerfect Gas(T,p,x)RT).

where gi is the Gibb's energy and x is the composition. The fugacity therefore has the dimension of a pressure, and the fugacities of each component are equal at equilibrium:


and have to be determined for each component in each fluid phase. For stoichiometric solvers, considering a specific reference state via the Henry's law constant and the equality between fugacities produces a mass action law:


where xi and yi are the mole fractions of species i in the aqueous and the gaseous phases respectively, and Kih is the Henry's law constant of species i. The Henry's constants could be expressed as a function of temperature, pressure, and even salinity. In addition, it can be derived from the constants expressed at the vapor saturation pressure psat and the molar volume at infinite dilution υi:


Detailed derivations of vapor–liquid equilibria, concepts, and techniques were reported by Firoozabadi (1999) and Michelsen and Mollerup (2007).

Equation of state.

The EOS for the gas phase is solved with respect to the gas volume at fixed pressure or the pressure at constant volume. The fugacity coefficients can then be calculated. A perfect gas EOS is not sufficient when addressing relatively significant pressures and temperatures, and it is necessary to consider non-ideal gases. The cubic EOS is useful for reactive transport, given that the cubic equation can be directly solved, e.g., the use of an analytical Cardan method (Nickalls 1993, Spycher et al. 2003). When considering a mixture, a mixing rule can be employed, in addition to calibrated binary interaction parameters, thus rendering the cubic EOS general and well-suited for multicomponent systems.

The classic cubic equations are the Soave–Redlich–Kwong (SRK) (Soave 1972) and Peng–Robinson (PR) (Peng and Robinson 1976; Robinson and Peng 1978) equations. Various improvements of the EOS were proposed using different alpha functions adapted for precise compounds and more complex mixing rules. The predictive PR equation (Qian et al. 2013) improves the mixing rules proposing the temperature dependent binary parameters calculated via the global contribution method (GCM). Hence, the GCM is advantageous, given that the vapor–liquid equilibrium (VLE) can be calculated for various systems without introducing new empirical fitting parameters. There are accurate virial EOS models such as the Benedict–Webb–Rubin EOS (Benedict et al. 1940), which have been widely used for CO2 (Duan and Sun 2003). However, they are typically different for each gaseous compound, and they require iterative solutions.

The classic cubic EOS formulation and mixing rule are general, i.e., applicable to every compound. Their ability to accurately represent the CO2–H2O and CO2–H2O–NaCl systems was demonstrated (Spycher et al. 2003; Spycher and Pruess 2005; Corvisier et al. 2014, 2017; Hajiw et al. 2018), in addition to other systems that include various compounds and gaseous mixtures. A comparison between experimental data and several other methods yielded good results (Hajiw et al. 2018).

Aqueous activity correction.

Activity coefficients in solution also require models. The simplest approaches (Debye-Hückel 1923; Davies and Shedlovsky 1964; Colston et al. 1990; B-Dot, Helgeson 1969) allow for the modeling of diluted solutions, as they are based on long-distance electrostatic interactions, which are predominant in this case. For highly saline solutions, wherein the ionic strength is higher that 3 M, short-distance non-electrostatic interactions should be considered. The Pitzer (Pitzer 1973, 1991) and specific ion interaction (SIT) (Brønsted 1922; Grenthe and Puigdomenech 1997) models yield relatively good results for these solutions.

For neutral species such as dissolved gases, the Pitzer model yields the following expression:


where λij, μijk, and ζijk are the interaction parameters between i, j, and k species.

The SIT model limits to binary interactions and appears similar to first-order Pitzer models. Activity coefficients for such neutral species are expressed as follows:


where εij are the binary interaction parameters between the i and j species.

Contrary to several activity models (Setchenow, Davies, or B-Dot), the two abovementioned models take into consideration and allow for the calibration of the influence of ions on the dissolved gas activity, and therefore, the gas solubility. In many CO2 solubility models, the activity of dissolved CO2 is calculated using specific models (Drummond 1982; Rumpf et al. 1994) based on available experimental data. Nevertheless, the SIT or Pitzer models allow for the use of the same model for various electrolytes and dissolved gases; thus simplifying the databases.

In dilute solutions, the activity of water is typically set as 1.0. However, its value deviates from 1 in saline solutions. The water activity can be expressed directly as a function of ionic strength or the osmotic coefficient (Helgeson 1969; Pitzer 1973, 1991), which gives reasonable results.


The quality of the representation of the gas–water equilibrium is significantly dependent on the parameters, and therefore on the database combined with the geochemical solver. Simulations that involve a poorly defined Henry's constant do not result in accurate predictions, and this is an important limit (mainly related to the current lack of sufficient and consistent experimental data) that modelers should consider.

With reference to the literature, the parameters associated with gases used in the EOS (critical temperature and pressure, and acentric factor) do not vary significantly among various sources, and they can be retrieved from the Yaws handbook (Yaws 1999).

The interaction parameters of gases, which are dependent on the selected EOS and mixing rule; and those of electrolytes, which are dependent on the selected aqueous activity model, come from experiments (binary mixtures PTy data and activity/solubility data, respectively), and are therefore significantly dependent on the quality of the selected data.

Henry constants or reference chemical potentials, and their dependence with respect to pressure are also drawn from solubility experimental data (Hajiw et al. 2018).

First 0D simulations.

Geochemical simulators, as described here, allow for the simulation of many gas–water–brine systems, including CO2–H2O and CO2–H2O–NaCl. Significant research attention has been directed to these systems, and consistent data are increasingly available to help us select the appropriate EOS and activity models as demonstrated here below.

As an example, simulations were run using Chess on CO2 solubility and co-solubility data at 50 °C, 100 °C, and 150 °C (Figs. 2 and 3), and the following modeling options were tested: the perfect gas EOS, PR EOS, and the pressure correction of thermodynamic equilibrium constants. These illustrate that it is necessary to use an accurate EOS, especially when the temperature and pressure increase; and highlight the importance of appropriately correcting the Henry's constants with respect to pressure. As an illustration, the absolute average deviation (AAD) between the solubility experimental data and the model decreased from 85% using the perfect gas EOS and no pressure correction to 4% using the PR EOS and pressure correction.

Other simulations were run with Chess on CO2 solubility data in NaCl brines at 50 °C, 100 °C, and 150 °C, and 150 bar (Fig. 4) using the B-Dot aqueous activity correction and the SIT model. These simulations highlight the importance of appropriately representing the electrolyte and its effects on the gas solubility (the salting-out effect). In this case, the AAD decreased from 7% with the B-Dot model to 3% with the SIT model. For NaCl brines, B-Dot simulations appropriately represent the data at significantly high salinities, but it may be noted that this is not the case for other brines (e.g., CaCl2, MgCl2).


Variations in the composition of phases have a direct influence on the saturation, solubilities and co-solubilities (e.g., the dew and bubble-point pressures and the PT phase envelope). Other thermodynamic properties such as the density and viscosity are also dependent on the pressure, temperature, and composition; and they feature in the fluid dynamics calculations. The degree of influence varies with respect to the type of fluid. For example, the density of a CO2 and SO2 mixture can significantly increase with an increase in the SO2 content. On the contrary, the density of a highly non-ideal mixture of CO2 and H2 can decrease with an increase in the H2 content. The densities of the CO2-rich mixture with 2–10 mol% of H2 increases in accordance with an increase in the pressure, and decreases in accordance with an increase in the temperature and the H2 concentration: at 313.15 K and 9 MPa, adding 2% of H2 to the pure CO2 decreases the molar density by 25% (Sanchez-Vicente et al. 2013). It is therefore necessary to consider these modifications and to accurately select a thermodynamic model.

The modeling of the density of pure bodies can be achieved by employing an EOS or empirical correlation; although there is no EOS that can universally predict the thermodynamic state and properties of any multicomponent system. For CO2, CO2–H2O, H2O–NaCl, and CO2–H2O–NaCl, several specific models were developed at certain temperature and pressure conditions. For example, the density and viscosity of formation waters in sedimentary basins (Adams and Bachu 2002; McCain Jr. 1991), the density of CO2–H2O (Garcia 2003) and the density of CO2, CO2–H2O, and CO2–H2O–NaCl (Duan et al. 2008). However, these models are not applicable to other multicomponent systems.

The GCM is typically an advantageous method for modeling of multicomponent mixtures, given that the form of the EOS and mixing rule is general for any number of components. The SRK (Soave 1972) and PR (Robinson and Peng 1978) yield relatively accurate estimations of gas densities for pure bodies and gas mixtures (Sin 2015). Nazeri et al. (2017) obtained good results using the PR and SRK models for the gas, supercritical, and liquid densities of multicomponent systems that contained CO2, CO, N2, O2, Ar, H2, CH4, C2H6, C3H8, SO2, and H2S. However, in the presence of water, the estimations would be significantly different. The predictions of the liquid densities for water systems made the PR or SRK EOS approaches inaccurate, as they are not adapted for highly polar systems.

Another GCM-based model is the volume translated PR (Ahlers and Gmehling 2001, 2002a,b; Ahlers et al. 2004), which is based on the pseudo-volume (Péneloux et al. 1982). The model therefore improves the prediction of volumetric properties. In addition, the model can be employed for the asymmetric VLE approach (Ahlers and Gmehling 2002a,b). The model demonstrates a good prediction of the liquid density of water systems with a rather wide range of components (Ahlers and Gmehling 2001; Sin 2015). Chapoy et al. (2013) developed a state viscosity model and density model with volume correction for the CO2, O2, Ar, and water system. The Cubic-Plus-Association EOS can also consider molecular structures with hydrogen bonding. Moreover, in combination with the PR model, a better modeling of both the VLE and liquid densities of water systems can be achieved. This was demonstrated for hydrocarbons and other compounds (Hajiw et al. 2015, 2018). The drawback is that it is computationally complex. Further development of GCM-based models is required, and more experimental data of multicomponent mixtures are necessary for validation.

Feedback on fluid and matrix properties

When CO2 dissolves into reservoirs, it acidifies the water and initiates dissolution/precipitation reactions with rocks. Furthermore, if the injected gas is initially dry (as is typical for captured CO2), water vapor is formed to satisfy the gas–water equilibrium conditions. Some of potential co-injected impurities modify and promote geochemical reactions. This indicates that gas and its geochemical reactions have a significant influence on the entire porous medium, e.g., host water/brine and rocks (Jung et al. 2013; Sterpenich et al. 2013; Lu et al. 2014; Pearce et al. 2015, 2016a,b; Lu et al. 2016) and engineered materials such as cement-based wells (Jacquemet et al. 2005). Models should be able to follow such modifications of the fluid and matrix and accurately quantify their effect on associated properties using a general or specific EOS.

Dissolution/precipitation reactions (not only induced by gas) modify hydrological properties such as the porosity, reactive surface area, and permeability at the pore- and core-scales (Gouze and Luquot 2011; Xu et al. 2017). Moreover, CO2 injection into a heterogeneous medium may create wormholes, depending on the relationship between geochemical reactivity and fluid motion (Snippe et al. 2017). At the Darcy scale, the change in porosity due to the dissolution/precipitation of mineral m is typically given by ϕ=ϕ0Σm=1Nm(cmcm0)/ρm, where ϕ0 is the initial porosity; cm0 and cm are the initial and current concentrations of mineral m (unit: mol per volume); and ρm is the molar density of mineral m. In addition, the porosity varies with respect to pressure: ϕ = ϕ*ecr(p−-p*), where cr is the rock compressibility; ϕ* is the reference porosity; and p* is the reference pressure. The intrinsic permeability–porosity relationship is generally represented by the Kozeny–Carman model. These relationships are actually more complex than this relatively simple description, especially at small scales. For example, Xu et al. (2017) demonstrated that mineral dissolution can enlarge pore throats, increase porosity and permeability, decrease capillary pressure, and influence relative permeability due to the injection of dissolved CO2 in a Chaunoy sandstone core from the Paris Basin. A depressurization test revealed that gas exsolution and mineral precipitation may occur. However, this is completely different from injection of sc/gCO2 (Jacquemet 2006; Jacquemet et al. 2008; Snippe et al. 2017), as experimental and numerical studies confirmed that the influence of scCO2 injection on porosity and permeability is negligible. The modeling of such geochemically induced changes requires further experimental studies, more data, and an in-depth investigation to better understand these couplings and the evolving matrix structure, e.g., moving boundaries, reactive surfaces, tortuosity, diffusivity and the mass and volume balances/constraints. Further discussion of the modeling of evolving matrix properties and implications with respect to CCS are discussed in Seigneur et al. (2019, this volume) and Cama et al. (2019, this volume), respectively.

Changes in the gas composition have an influence on the gas volume and partial pressures, and therefore the compressibility and viscosity. This impacts the fluid motion: mass fluxes, gravity effects, and viscous forces. Molecular fluxes depend on a composition of the system, Hence, the interplay between these forces evolves into a heterogeneous composition of fluids, e.g., solubility-induced partitioning (see Applications section for further discussion). A notable example of the composition-dependent mass fluxes is density driven flow, which plays an important role in geological CO2 storage. The formation water becomes denser with dissolved CO2. Therefore, the accumulation of dissolved CO2 at the top of an aquifer leads to instabilities. Upon reaching an onset time, denser fluid tends to sink. The convective motion created as a consequence enhances the dissolution rate. The density of fluids should then be accurately modeled to predict the onset time and convective fluxes, and to estimate convective dissolution trapping. Farajzadeh et al. (2007), Kneafsey and Pruess (2010), Neufeld et al. (2010), MacMinn et al. (2012), and Wang et al. (2017) experimentally confirmed the importance of convective fluxes. Several studies were carried out, which involved mathematical and numerical stability analyses of convective mixing for homogeneous and heterogeneous media, with and without capillary transition zones (Lindeberg and Wessel-Berg 1997; Xu et al. 2006; Hesse et al. 2008; Pau et al. 2010; Elenius et al. 2012, 2015; Li and Jiang 2014); for an anisotropic medium (Ennis-King and Paterson 2005); in open and closed systems (Riaz et al. 2006; Wen et al. 2018) and 3D heterogeneous domains (Green and Ennis-King 2018). At present, there is insufficient information as to the influence of gases other than CO2 on convective mixing. Co-injected impurities can significantly change gas and liquid densities, affect the convection flux, and increase or decrease the onset time (see Applications section for further discussion). Coupled geochemical reactions and multiphase flow also influence gas exsolution, which can be induced during depressurization. Experimental studies have shown that the relative permeability for CO2 exsolved from brine curves is significantly different from those of drainage process (Falta et al. 2013; Zuo et al. 2013, 2017; Xu et al. 2017).


The previous section indicates that MMF&RT is a complex problem including mathematical and numerical issues with respect to the composition of flow and reactive transport. There are several formulations of MMF that attempt to accurately solve the evolving flow system and phase equilibrium. For MMF&RT, two coupling families exist: the global implicit and operator splitting approaches. Both families have advantages and drawbacks, and, the development of new coupling methods, formulations, and numerical techniques continues.

Multiphase multicomponent flow (MMF) approaches

The balance equations for multiphase multicomponent flow (Eqn. 26) can be formulated on mole fractions xα,k of species k in phase α, as follows:


where ρmol,α is the molar density of fluid α.


The balance system is equivalent to the set of mass conservation equations (26) and can be expressed with respect to concentrations, given the following:


where Mk is the molar mass of species k.

For the system of Nc conservation equations (Eqn. 35), there are 2Nf + NfNc unknowns, which are referred to as natural variables (Coats 1980): Sα, pα, and xα,k. Considering the closing relationships for pressure (Eqn. 15), saturation (Eqn. 1), and mole fractions (Eqn. 36) (2Nf equations), the balance system (Eqn. 35) should be complemented by Nc(Nf − 1) phase equilibrium relationships. Using the Gibbs phase rule, an isothermal system with a locally known number of phases Nf can be described using a minimum of NcNf + 2 variables. For an isothermal two-phase system, Nc variables are required to define a thermodynamic state. There is a large selection of variables. The most evident is the pressure and Nc − 1 variables of composition. The discretized nonlinear conservation equations should therefore be solved on the selected primary variables by employing an iterative method such as Newton's method. The residual is linearized, and the Jacobian is assembled for the iterative procedure. Approaches of MMF modeling are discussed further in the Primary variables section.

Phase equilibrium.

The phase equilibrium of fluid phases can be solved using an EOS and/or activity models depending on the selected approach, (Water–gas equilibrium section). Another method, which is referred to as flash calculations, is also widely used. As it is based on K-values, which are defined as Ki = yi / xi, the data are generally available and tabulated for hydrocarbon mixtures; where xi and yi typically denote the liquid and gas mole fractions, respectively. The K-values can be derived from the (φ−φ) and (φ−γ) approaches:


It can be noted that the K-values are therefore composition-independent under the hypothesis of an ideal mixture. The phase mole fractions V and L are defined such that V + L = 1. This allows for the introduction of a new variable zi, which is the overall mole fraction of species i:


The mole fractions can then be expressed with respect to zi, Ki, and V:


The Rachford–Rice equation yields a constraint on the mole fractions of the system:


which is solved with respect to V at a known pressure, temperature, and composition z using Newton's method. The mole fractions are then updated, and the phase equilibrium relationships are calculated to provide new K-values. The process is repeated until convergence occurs. The flash problem and stability analysis were extensively described by Firoozabadi (1999) and Michelsen and Mollerup (2007).

The phase equilibrium problem can be solved directly with the system of compositional flow (e.g., it is added to the Jacobian structure), or separately for each cell using any phase equilibrium method. An important advantage of using EOS models is that most of them can be solved analytically, e.g., for the cubic EOS. The flash problem requires an iterative method when a real mixture is modeled; otherwise, the K-values are constant, and the solution is significantly simplified.


The balance equations (35) are nonlinear and strongly coupled due to the relative permeability, capillary pressure, and fluid properties. The equations are (near-)parabolic, which is similar to the system of immiscible two-phase flow (Eqn. 16). The fully implicit method (FIM) is generally employed for the time discretization, e.g., all the equations are solved simultaneously. The fluxes are treated using the two-point or multi-point flux approximation. The transmissibility approximation affects the numerical stability and accuracy. Implicit upstream ensures the stability and increases the risk of truncation errors at the same time (Blair and Weinaug 1969; Allen 1984). The relative permeabilities, viscosites, and densities at the interface between two cells follow the flow direction uαnij > 0, where nij is the outer normal to the interface. The phase densitites in the gravity term ραg can be weighted with respect to the effective phase volume (Coats 1980). When modeling a heterogeneous medium, the discontinuity can be handled with harmonic averaging applied to Kkrα / μα or just using the intrinsic permeability. The system of nonlinear equations obtained from the time and space discretization should then be linearized with respect to the selected unknowns for Newton's method. The obtained linear system can be solved using a generalized minimal residual method (GMRES) algorithm (Saad and Schultz 1986) with a pre-conditioner such as ILU0 incomplete factorization (Saad 2003).

When the coupling between equations is weak, they can be solved separately. Similar to a decoupled pressure equation of immiscible flow (Eqn. 19), the equations in (35) are summed to derive a pressure equation. The total and phase velocities are defined using the fractional flow (Eqns. 20 and 21). This evolves into one pressure equation and a set of saturation/composition equations that are weakly coupled, and can therefore be solved semi-implicitly using implicit pressure explicit composition/saturation (IMPEC/IMPES). The pressure equation is first solved, followed by the composition equations with fixed fluxes. Under the assumption of incompressible flow, the pressure equation is elliptic (Eqn. 19). If the capillary pressure is neglected, the saturation/composition equations are purely hyperbolic and nonlinear in comparison with the linear transport equations. Given that this decreases the dependencies of the equations, they can be separated and efficiently solved using different schemes and preconditioners that are appropriate for each type of equation (Chen et al. 2006). However, if the equations are highly coupled, the solution of sequentially decoupled equations leads to the retardation of property calculations.

Despite the computational complexity of the implicit methods, they allow for a large time step, which is advantageous. The scheme is unconditionally stable; however, truncation errors may be introduced at large time steps. To reduce the numerical complexity, the flow system can be reformulated for the IMPEC. The explicitness indicates a restriction on the time step controlled by a Courant–Friedrichs–Lewy (CFL) criterion. Formulations of the CFL condition were detailed by Settari and Aziz (1975), Peaceman (1977), and Coats (2003). The adaptive implicit method combines both methods using different levels of implicitness for each cell (Thomas and Thurnau 1983; Cao 2002).

In general, codes employ the finite volume method (FVM), as it is conservative and easy to implement. The main drawback of FVM schemes is that they are not adapted for modeling complex geometries and anisotropic porous medium. The FVMs were trevised to produce a new category of cell-centered finite volume methods (Hermeline 2007; Aavatsmark et al. 2008; Eymard et al. 2009). The methods are based on multi-point flux approximations, they ensure continuous pressure across the control volume interfaces, and are consistent on triangular meshes (Friis and Edwards 2011; Pal and Edwards 2011; Eymard et al. 2012). Finite element methods can result in violation of local conservation properties, which is important when modeling flow and reactive transport. Mixed finite element and discontinuous Galerkin methods can ensure local mass conservation (Raviart and Thomas 1977) and (Wheeler and Yotov 2006); the latter was extended to hexahedral meshes by Ingram et al. (2010).

Primary variables.

Natural variable formulation (NVF) is widely used due to its direct implementation: Sα, pα, and xα,k; although it requires variable switching when a phase appears or disappears (Coats 1980). When the primary variables are changed, the Jacobian should be reassembled to avoid a set that is numerically inefficient and may lead to inconsistencies. The NVF is not unique. Various formulations were developed by combining/replacing variables/equations to solve the problem with respect to the mass or volume variables. Several alternatives to the NVF are reviewed below.

The pressure pα and component masses mα formulation (Acs et al. 1985) was designed for the compositional flow using Nc + 1 primary variables: a pressure and the component masses in moles, and mi=ϕVΣαρmol,αSαxα,i. Replacing the mole fractions by the component masses implies an additional equation that constraints the volume. This results in Nc + 1 equations.

The pressure pα and overall composition zi (Eqn. 39) formulation exploits the expression of the thermodynamic state and composition flow with respect to pα and zi. This can be achieved using the flash calculations and by re-writing the flow with respect to zi. Other molar formulations with respect to pressure, phase fractions V and L, zi, and the composition of Nf − 1 phases or ln Ki are also possible.

The overall composition formulation with the flash calculations was compared to the NVF. The performance was tested on a problem related to CO2 injection (Lu et al. 2010). A higher convergence rate of the overall composition formulation than the NVF with switching variables was observed. Voskov and Tchelepi (2012) compared the NVF, molar, and overall compositions using several problems that prioritize the NVF.

A tie-line based method was employed to parameterize the compositional space of the thermodynamic equilibrium problem (Entov et al. 2002). The technique was extended to an EOS-based two-phase multicomponent flow system, where the tie-lines are adaptively constructed, thus replacing or/and accelerating the phase stability calculations (Voskov and Tchelepi 2009).

A formulation with respect to the pressure and total concentrations using a concept of extended saturation was presented by Abadpour and Panfilov (2009). Voskov (2012) employed this approach for modeling of EOS-based compositional flow coupled with the negative flash method developed by Whitson and Michelsen (1989). The phase fraction and saturation can be negative or larger than one.

Extended NVFs were devised using Henry's law/solubility and the inverse capillary pressure function. First proposed by Ippisch (2003), the use of the solubility as a function of pressure was extended by Bourgeat et al. (2009): ρl,i = ρlXl,i = ρl,i(pl + pc(Sg)), after which the saturation Sg = Sgl,i, pl). Jaffré and Sboui (2010) introduced three variables with complementarity conditions for two-phase two-component flow: pl, Sl, Xl,i. Angelini (2010) used extended pressures pl, pg (PPF). Neumann et al. (2013) considered pc, pg as primary variables using the inverse capillary pressure and nonlinear solubility functions.

The pressure pα, saturation Sα, and fugacities fα,i (PSF) formulation given by Lauser et al. (2011) proposes a fixed set of primary variables. The phase transition is controlled by complementarity conditions, and the flash calculations are excluded. The problem is solved using a semi-smooth Newton's method. Fugacity coefficients ϕα,i are calculated using EOS models.

Masson et al. (2014) presented a comparison between the NVF, PSF, and PPF for the modeling of different physical problems, namely, a drying process by suction, gas injection, and gas migration in a heterogeneous medium. The convergence test results revealed that the NVF and PSF demonstrated better performances than that of the PPF. Ben Gharbia et al. (2015) studied the NVF coupled with the negative flash and the PSF coupled with the cubic EOS. There were numerical difficulties with respect to the PSF during the phase state transition, and it was slightly less efficient than the NVF, which is still a reliable option. The practical usage and intended applications can serve a guide for the selection of a convenient set of primary variables.

Coupling MMF&RT

Coupling methods of MMF&RT should fully utilize the accumulated knowledge of reservoir engineering and reactive transport modeling. Their mathematical and numerical approaches were designed for the solution of different problems at different scales. Reservoir simulations are typically large-scale oriented and based on the global implicit approach (GIA) using flash calculations and EOS. Most reactive transport codes decouple the transport and the chemistry by the operator-splitting approach (OSA), thus rendering it flexible for extension. Coupling MMF&RT inherits all the complexity of both problems and the experience of both communities

Global implicit approach.

The global implicit approach consists of the simultaneous solution of the nonlinear and linear partial differential equations (PDEs) and ordinary differential equations (ODEs) of transport/flow and the chemical system. Another method for solving the coupled system, which is referred to as the direct substitution approach (DSA), is derived from the RT modeling. The chemistry variables and the basis species concentrations C replace the total concentrations T in the transport equations. The total concentration of basis species is the mathematical sum of the basis species concentration over the entire system that can be constructed using the tableaux method described by Morel and Hering (1993). The chemistry is partly integrated into the transport equations, which are solved together with the heterogeneous reactions with respect to the basis species. Lichtner (1996) presented a general DSA for the solution of the RT problem, which can be applied to multiphase multicomponent systems. The concept can be expressed as follows:


where F is the linear transport or MMF operator; Seq and Skin typically denote the stoichiometric matrices; and Req and Rkin(C) are the reaction rates for the equilibrium and kinetic reactions. The formula was modified after Kräutle and Knabner (2007). The MMF&RT problem can then be solved with respect to the basis species C using Newton's method.

Global methods attempt to eliminate the equilibrium/kinetic rates. Sevougian et al. (1993), Lichtner (1996), Molins et al. (2004), Steefel and MacQuarrie (1996), and Saaltink et al. (1998) demonstrated that this can be achieved by the construction of an orthogonal matrix or linear combination of equations. For example, ESeq = 0, where E is the equilibrium rate annihilation matrix. These matrix transformations simplify and can reduce the number of differential equations; however, the equations are still coupled if there are heterogeneous reactions between the aqueous solution and solids. For example, a mineral can precipitate, dissolve, and re-precipitate depending on the aqueous composition and saturation. Kräutle and Knabner (2007) found a method to decouple a subsystem of mobile species equations from the global problem. Two sets of variables were introduced: a linear combination of only mobile species, and a linear combination of only immobile species. Hence, there was no combination of mobile and immobile species. This substitution was computationally advantageous, and the PDE system was well structured. The method was then revised to model the disappearance or appearance of minerals, and it was recently extended to two-phase modeling (Kräutle 2011; Brunner and Knabner 2019). Issues due to the disappearance of minerals and gases were investigated for the geochemical solvers and reactive transport codes, and mathematical and numerical solutions were proposed (Kräutle 2011; Leal et al. 2013)

Reduction methods of the DSA have been used to couple the MMF and geochemistry. Nghiem et al. (2004) applied a classic DSA to the multiphase multicomponent problem under the consideration of modeling both equilibrium and kinetic reactions, dissolution and precipitation of minerals, and porosity changes due to changes in the geochemistry and permeability using the Kozeny–Carman model. The importance of efficient linear solvers was highlighted, as they are required for the solution of the resulting sparse matrix. The mass balance equations with respect to the basis species are completed with a volume constraint equation, phase equilibrium, and kinetic rate relationships. The DSA-based integration of RT in the code relies on the existing numerical methods, such as the adaptive implicit method (AIM) developed for a compositional simulator generalized equation-of-state model for greenhouse gas (GEM-GHG). Another example of a matrix reduction method was presented by Fan et al. (2012). The balance equations are transformed by the elimination of equilibrium rates. In addition, the resulting system is solved with respect to the elements rather than the primary species, as the composition of all the species can be expressed with respect to the elements/atoms and corresponding number of atoms (Sevougian et al. 1993). The element-based method was developed in the automatic-differentiation-based general purpose research simulator (ADGPRS) framework using a two-stage pre-conditioner (Cao et al. 2005) and the AIM.

Despite the linear property of the transport operator, several reactive transport codes opted for the global implicit approach (GIA). MIN3P (Mayer et al. 2002); and PFLOTRAN (Lichtner et al. 2015) use this approach as well as CrunchFlow (Steefel 2009). Moreover, CrunchFlow and PFLOTRAN have an OSA option. Only PFLOTRAN has an MMF&RT coupling; and the GIA and OSA exist. Nevertheless, the reservoir modeling of large nonlinear systems using advanced solution techniques allows for the extension of fully implicit MMFs to chemical reactions. Although the computational complexity of the GIA is transparent, the increase in computational facilities allow for exploitation of the accuracy and robustness of the GIA. In addition to the abovementioned couplings, GEM-GHG (Nghiem et al. 2004) and ADGPRS (Fan et al. 2012) are coupled with GFLASH, which provides the equilibrium calculations based on the Gibbs free energy minimization (Voskov 2017). A fully implicit EOS-based compositional flow simulator, which is referred to as GPAS, was also extended for chemical flooding (Pope et al. 2005). However, the chemical reactions are modeled in the aqueous phase only. The coupling was achieved using a FIM and hybrid method. The hybrid option can be used to implicitly solve the subset of balance equations for hydrocarbons, and then explicitly solve the remainder. A concept of dominant and minor components is generally used to decouple flow and transport equations, as discussed below.

Operator splitting approach.

Similar to the classic OSA used for RT coupling, the resulting system of MMF&RT can be decoupled by alternating the MMF and RT operators. The flow system is first solved, thus yielding the intermediate pressures, saturations, composition, velocities, fluid properties, and other secondary variables. The compositional formulation can be used, thus making it applicable to existing reservoir simulators and geochemical codes. For the extension of geochemical codes, it is sometimes assumed that there is a set of dominant components responsible for MMF, while other components are minor, and their impact on the flow can be neglected. This reduces the number of nonlinear equations, and only linear transport is solved for minor components. For example, two components, namely, CO2 and H2O, can be selected as dominant for the two-component two-phase flow. However, this assumption can limit the modeling of the complex dynamics of multicomponent systems (see Applications section for further discussion).

The reactive transport part ℜ may be sequentially linked and can be solved using a GIA or OSA (Yeh and Tripathi 1989). This allows for existing RT codes to be used as input/output modules. A general OSA has the following form:


where n is the time step number, and k corresponds to an iteration number. In the non-iterative coupling (SNIA), k = 0, and the system proceeds to the next step after the RT. Alternatively, the objective of the sequential iterative approach (SIA) is to iteratively solve the flow for an updated chemical state, and then for the RT until convergence. Using the SNIA leads to persistent errors for reactive problems (Valocchi and Malmstead 1992). Another drawback is the time-step restriction, which is necessary to ensure accuracy of information transfer between the operators. The reliability of the OSA, and that of SIA in particular, was demonstrated for reactive transport and MMF&RT problems (Carrayrou et al. 2010; Saaltink et al. 2013). The chemistry has a significant influence on the fluid displacement, and vice versa. The OSA-based coupling depends on the applied internal methods, their connectivity, and on the modeled physical process.

Flexibility is a critical feature of the OSA, where new functionalities are easy to implement, and the operators/modules can be verified independently. The method is extensively applied in RT codes (Steefel et al. 2014). Several geochemical codes employed the method for modeling unsaturated flow, such as the bubble model developed in MIN3P (Mayer et al. 2012; Molins and Mayer 2007).

There are several examples of OSA-based couplings between MMF and RT:

These examples illustrate that most MMF&RT couplings are not iterative between flow and RT. This may limit chemical processes, and therefore their influence on thermodynamics and hydrodynamics, similar to SNIA-based coupling for RT. The composition-dependent properties should be modeled implicitly to prevent a discrepancy between the MMF and RT operators. Consequently, the FIM is a favorable option. Implicitness provides codes with the unconditional stability property; thus, large time steps are allowed. This is of practical usage for modeling the rapid mechanical displacement of fluids. However, the reactive time step can be a restriction. When the range of physical processes occur over different characteristic scales, time-step splitting between operators can be a solution. The connecting properties should accurately be updated, and mass conservation should be strictly controlled. Another numerical alternative is the AIM that reduces the nonlinear flow system. Furthermore, the space domain decomposition is efficient when the computational intensity over the model domain is highly heterogeneous.

A tight SIA-based coupling, which couples flow and RT operators iteratively or with implicit methods, is not frequent in MMF&RT codes because of the large nonlinearity of the flow problem. However, it could be a good alternative, since simplified flow formulations typically have restricting assumptions; and they are not suitable for general-purpose simulations. The tight SIA-based approach with high implicitness is critical for the modeling of reaction-driven hydro- and thermodynamic processes.


The modeling of complex thermodynamics is a necessary condition for MMF&RT. Multicomponent systems at supercritical conditions cannot be described using the ideal laws. Modeling mixtures requires implicitly updated fluid properties and relevant approximations of real fluid behavior. For example, an accurate modeling of compressible multicomponent real fluids can determine the evolution of fluid phases, and their displacement and influence on the matrix. The asymmetric approach of VLE allows for the modeling of various systems and gains from both constantly developing EOS and activity models. Based on the analysis of the MMF&RT methods and experience over the last four decades, the GIA or a tight SIA-based method with NPV or equivalent formulation can be considered a reliable option. The tight SIA can be preferable for reactive transport simulators. Due to the complexity of the MMF&RT structure, it is necessary to benchmark coupling methods. Although a few MMF benchmark problems exist (Pruess et al. 2004; Class et al. 2009; Nordbotten et al. 2012), they are CO2-oriented and have only one gas species in the gas phase. Multicomponent benchmarks (Sin et al. 2017b) are therefore required to better understand the structures and behavior of different coupling methods.

Coupled physical and chemical processes at the continuum scale also occur at the pore scale, and they can be modeled by approaches such as the pore network model, Lattice Boltzmann model, direct numerical simulation, and particle methods. The advancement of computational capabilities allows for the application of these methods in complex physics and geometries. The Shan–Chen multicomponent lattice Boltzmann model was efficiently applied to fluid–fluid and fluid–solid interfacial reactions in miscible two-phase flow (Chen et al. 2018; Li et al. 2018) and reactive transport with evolving porosity (Gao et al. 2017). For example, Li et al. (2018) illustrated the influence of fluid–solid interface relations (e.g., adhesive forces, specific surface area) on the effective permeability and capillary pressure. Recent studies also show possibilities to model physics across scales integrating over a general framework. Battiato et al. (2011) proposed a hybrid method, wherein an iterative procedure couples two models and defines concentrations and fluxes through the boundaries of pore and continuum scale domains. A further investigation of hybrid methods is necessary to better understand the physiochemical processes that occur across scales and contribute one to each other (Molins and Knabner 2019 this volume).


The increasing interest in the CCS field has initiated an intensive investigation of subsurface processes including numerous experiments and field studies, in addition to mathematical and numerical methods. Among various subjects, captured CO2 gas contains small amounts of impurities such as N2, O2, Ar, CO, SOx, NOx, CH4, H2S, and H2 (Kather 2009; IEAGHG 2011; Ussiri and Lal 2017). The type and concentration of co-injected impurities are dependent on the source of CO2 and the selected capture method. The presence of impurities modifies the solubility of CO2, and therefore influences the storage capacity (Li and Yan 2009a,b; Zirrahi et al. 2012; Ziabakhsh-Ganji and Kooi 2012, 2014). Moreover, The presence of impurities also affect the geochemistry and fluid dynamics of the system (Xu et al. 2007; Xiao et al. 2009; Corvisier et al. 2013, 2014, 2017; Jung et al. 2013; Sterpenich et al. 2013; Bacon et al. 2014; Lu et al. 2014; Ziabakhsh-Ganji and Kooi 2014; Pearce et al. 2015, 2016a,b; Wei et al. 2015; Lei et al. 2016; Lu et al. 2016; Waldmann and Rütters 2016; Todaka and Xu 2017).

Chromatographic effects have been observed in the field and laboratory following the injection of CO2 gas with impurities. Wei et al. (2015) reported the injection of a CO2-rich gas mixture with N2 and O2 during a Tongliao pilot experiment in China. Earlier arrivals of co-injected oxygen and nitrogen gases were observed when compared with the arrival of CO2. Inversely, a gas mixture of 98% CO2 and 2% H2S was injected into a depleted gas reservoir in Canada. The earlier arrival of CO2 and delayed breakthrough of H2S were recorded at producing wells. The chromatographic effect is significant at large scales (Wei et al. 2015; Paterson et al. 2013). A series of dynamic solubility experiments was carried out by Bachu and Bennion (2009), which demonstrated the partition of a gas mixture because of the differential solubility of gases. Numerical models were able to reproduce these observations. Gas partition also occurs when other gases are already present in the porous medium.

Convective dissolution, which is an important trapping mechanism, is also significantly dependent on the presence of impurities because of their influence on the density of the plume. The interplay between gas partition and density-driven flow was demonstrated in three cases of CO2 injection. and the modeling of MMF&RT allowed for the investigation of the critical effects of co-injected impurities on the evolution of the systems (Xu et al. 2007; Xiao et al. 2009; Lei et al. 2016; Todaka and Xu 2017).

Phase equilibrium and solubility

Batch geochemical simulations (0D simulations) offer a great opportunity to verify the ability of the abovementioned approaches to accurately simulate gas–water interactions in closed and homogeneous systems prior to the modeling of reactive transport processes at a larger scale. CO2–H2O and CO2–H2O–NaCl systems were studied intensively, and the consistent data obtained from the literature allow for the testing of the geochemical solvers and the associated databases. CO2 solubility and co-solubility data at 50 °C, 100 °C, and 150 °C (Figs. 56) were modeled using three geochemical solvers and one thermodynamical model:

The four models were in relatively good agreement with the CO2 solubility data for the entire range of the temperature and pressure. For the co-solubility data, Chess and Phreeqc gave results in good agreement with the data, while results from Duan and Gem–Selektor showed significant discrepancies with the selected data (Figs 5 and 6). The absolute average deviation (AAD) between the experimental solubility data and the model decreased from 4% using Chess, Gem–Selektor, and Duan to 2% using Phreeqc. AAD for co-solubility data decreased from 43% using Gem-Selektor, 22% using Duan, 14% using Phreeqc to 13% using Chess.

Geochemical solvers can be used to solve gas–water equilibria and the entire speciation. Consequently, Chess, Phreeqc, and Gem–Selektor can be used calculate the pH of CO2–H2O systems. The extension of the approach by Duan and Sun (2003), which was developed by Li and Duan (2007) can be used to calculate the pH of such systems. These four models were used to simulate pH measurements at 50 °C, 100 °C, and 150 °C (Fig. 7). The AAD between the experimental pH values and the model predictions decreased from 2% using the Chess, Gem–Selektor, and Phreeqc models to 1% using the model by Li and Duan (2007). Various gas–water systems and gas mixtures–water were successfully simulated using Chess (Figs. 8 and 9). Hajiw et al. (2018) presented a comparison between the geochemical approach and process homogeneous approaches using the same PR EOS, and then illustrated the good performances of the geochemical solvers applied to these systems. The AAD for these systems are 3% for Ar, 10 % for CO, 6% for CH4, 8% for H2, 8% for H2S, 4% for N2, 3% for O2, 11% for SO2, 7% for CO2–CH4, and 3% for CO2–N2. Other simulations were run using the abovementioned models to reproduce the CO2 solubility data measured in the presence in NaCl brines at 50 °C, 100 °C, and 150 °C at 150 bar (Fig. 10). The AAD was 12% using Gem–Selektor, 8% using Phreeqc, 4% using Duan, and 3% using Chess. The performances of the four models were good, although Gem–Selektor exhibited a slight overestimation of the CO2 solubility at high salinities, which can be attributed to the employed aqueous activity model.

Several pH measurements of CO2–H2O–NaCl systems at 50 °C and 1 M, 3 M, and 5 M of NaCl were also simulated (Fig. 11) The AAD between the experimental pH and the model decreased from 27% using Gem–Selektor model to 3% using the Chess, Phreeqc, and Li and Duan model. The overestimated CO2 solubility of Gem–Selektor probably leads to this rather important AAD for pH.

Overall, the benchmark presented here showed that that a variety of systems (gas mixtures, water, mixed brines, etc.) can be accurately simulated using geochemical solvers implemented in reactive-transport codes.

Gas chromatography

Dynamic solubility experiments were carried out by Bachu and Bennion (2009) to better understand the influence of a gas type (H2S, N2, SO2, CH4) and its concentration (2, 5, 30% of H2S) on gas separation at different temperature, pressure, and salinity conditions (61 °C and 13.5 MPa; 25 °C and 6 MPa; fresh water and 119 ppm of salinity). Moreover, the breakthrough of gases was recorded and analyzed. The experimental results revealed that the breakthrough time and profiles of the gas components are dependent on the solubility ratio between gases. Less soluble gases arrive earlier with higher breakthrough concentrations.

In these experiments, mixtures of H2S and CO2 were injected in a tube homogeneously packed with silicate-rounded sand. Two feed gas cases were used: 2% H2S and 98% CO2, and 30% H2S and 70% CO2. The gas was injected at a constant rate of 7.5 cm3/h into the tube at a constant pressure of 13.5 MPa, temperature of 61 °C, and brine salinity of 118950 ppm (Bachu and Bennion 2009). The experiments were modeled using HYTEC (Sin and Corvisier 2018) and CMG-GEM (Bachu et al. 2009).

The modeled gas saturation profiles and aqueous composition of both codes for the co-injected 2 % H2S and 98% CO2 differed only slightly. The fluid displacement modeled by the HYTEC lags the curves of the CMG-GEM (Figs. 12a,b). This can be attributed to the different solubility models used in the codes. Model results illustrate the gradual evolution of the differential solubility of H2S and CO2 into an increasing difference between the H2S and CO2 profiles, as shown in Figure 12b (the mole fraction of dissolved H2S is scaled). Since H2S is more soluble than CO2, the breakthrough of CO2 was prior to that of H2S (Fig. 13). The experimental and numerical results were in good agreement. They illustrated a later arrival of H2S. After the breakthrough, the H2S concentration increased and reached its concentration in the feed gas.

The same analysis of gas chromatography was applied for 95% CO2 mixtures: H2S, N2, SO2, and CH4 (Bachu and Bennion 2009). Four numerical models were run using HYTEC. The numerical and experimental results of breakthrough of gases were compared, as shown in Figure 14a: the mole fraction of the impurity gas phases H2S, N2, SO2, and CH4 were normalized by the relevant initial mole fraction of the injected gas (5%).

As can be seen from the figure, the less soluble N2 and CH4 arrived first and exhibited higher peaks than that CO2. The opposite was observed for H2S and SO2. The numerical results indicated the phenomenological behavior of fluids. The modeled height of the peaks was not precise. During the experiment with 5% SO2, it was difficult to detect SO2. Moreover, it was found in a very small amount in the gas phase, as it was almost completely dissolved in the brine.

Although the numerical results were slightly different from the experimental data (Figs. 14a,b), and there were several uncertainties in the data, the numerical models present a good representation of the gas chromatography, e.g., the gas partition process of mixtures due to the differential solubility. It should be noted that the gases are displaced with the common velocity, which is the Darcy's velocity of the gas mixture. During the passage, several molecules of a substance are dissolved in the liquid phase, and other molecules are carried by the general motion of the gas mixture. This is similar to the chromatography of solutes when the retardation is due to ion exchange. The accurate modelling of solubility is therefore necessary for the gas partition. The co-injected impurities can be classified in decreasing order of their solubilities (13.5 MPa, 61 °C), as follows:


The above corresponds to the breakthrough order of the gases (Fig. 15), a similar test with 5% Ar was performed.

Gas chromatography and density-driven flow

Co-injected gases can accelerate or decelerate the plume displacement, affect the volume capacity, react with the host rock, or be used as tracers for monitoring. They are often modeled as dissolved in formation water (Xu et al. 2007; Xiao et al. 2009; Lei et al. 2016; Todaka and Xu 2017). This assumption changes the gas dynamics, mass/molar fluxes, and neglects the gas partition in the plume. Another important physical phenomenon is convective mixing. The impact of impurities on convective flow was studied by Li and Jiang (2014) and Raad and Hassanzadeh (2016), who analysed effects of impurities on the density driven flow, onset time, and flux rate. However, with reference to the literature, there is no model of impure CO2 injection that considers the density of mixtures and convective dissolution at a large scale.

Impurities in the injected stream have an influence on the solubility of CO2, the gas dynamics of the plume, and the gravitational forces of the entire system; since the densities of mixtures can vary significantly with respect to their compositions. For example, in the presence of N2, O2 and Ar, the CO2-rich mixture is lighter than the pure CO2 at high pressure and temperature conditions. In contrast, the SO2 + CO2 mixture is significantly heavier than the pure CO2. Based on the analysis of the volumetric properties (Ahlers and Gmehling 2001; Al-Siyabi 2013; Li and Yan 2009b; Nazeri et al. 2017; Sin 2015), and under the consideration of previous observations of gas chromatography (Eqn. 44), three cases of gas injection were devised to study the gas partition and density driven flow. In the base case, the pure CO2 was injected into a homogeneous long confined aquifer with a height of 100 m and depth of 1100 m. The fluid and matrix properties were obtained from Sin and Corvisier (2018). The injection rate was constant: 10 kg/s or 0.03 m3/s for the pure CO2 scenario. Two models of impure CO2 were then designed with the same volumetric rate (but different mass rates): injection of 95% CO2, 4% N2 and 1% O2 mixture (air and CO2 model), and 95% CO2 and 5 % SO2 mixture (SO2 and CO2 model). The length of the aquifer was 10 km in the present models. The grid resolution was 38400 cells, with a length that varied from 2–1000 m. The axisymmetric geometry allowed for an adequate representation of the convective mixing (Sin 2015). The number of gas components in the models varied from 2–4 (H2O(g) included). Simulations were carried out using HYTEC (van der Lee et al. 2003; Lagneau and van der Lee 2010; Sin et al. 2017a).

The shape and general dynamics of the gas plume were similar in all cases, as expected (Fig. 16). The gas chromatography effect resulted in the accumulation of the less soluble N2 and O2 at the leading edge of the plume. In the SO2 + CO2 case, the delayed arrival of a more soluble SO2, which was concentrated near the injection well and lagging behind the CO2, was expected (Fig. 17). The gas partition occurred throughout the plume. However, it was most pronounced at the top of the aquifer due to the advection forces that continuously supply less soluble gas compounds (compare Figs. 17, 18, 19). Heterogeneous gas distribution contributes significantly to the thermodynamic properties such as density. The gas density maps reveal the dependence of the volumetric properties on the gas composition (Fig. 20).

The formation water that contained dissolved CO2 was heavier than the formation water without CO2. The liquid density of the formed current increased with an increase in the content of CO2, which was not the case for all the compounds, such as H2S. The difference in density creates instabilities that balance the gravity forces that can finally initiate the convective motion under the current at the onset time. Several studies were conducted, which involved the experimental and numerical analysis of the nature of the process in closed and open systems for homogeneous and heterogeneous media (Lindeberg and Wessel-Berg 1997; Elenius et al. 2015; Emami-Meybodi et al. 2015; Green and Ennis-King 2018; Riaz et al. 2006; Thomas et al. 2018; Wen et al. 2018). The presence of impurities changes the densities of the gas and aqueous phases, and the convective fluxes then differs depending on a type and quantity of impurities. In the present models, the small amount of SO2 (5%) resulted in a denser aqueous solution, and the increased convective motion accelerated the dissolution of gases.

Figures 21, 22, and 23 illustrate the effect of the gas partition on the distribution of the compounds dissolved in formation, which impacts the general evolution of the system. The liquid density of the SO2 + CO2 model was almost greater than that of the Air + CO2 model by a factor of 2, which indicates more CO2 should be dissolved. In addition, the gas plume of the Air + CO2 model was less dense than that of the SO2 + CO2 model. Consequently, the gas current that contained air propagated faster than the gas current of the SO2 + CO2 model (Fig. 24a). With respect to the CO2 inventory (Fig. 24b), the ratio of CO2(g) over the total CO2 was very similar for the pure CO2 and SO2 + CO2 models. However, the quantity of injected CO2 in the SO2 + CO2 model was significantly greater than that in the pure CO2 model due to the different densities of the injected streams (Fig. 24c).

The positive effect of SO2 on the convective dissolution was evident. However, it should be noted that the convective rate decreased afterwards, as the formation water underneath the plume became homogeneously saturated in dissolved gases over time. The accepted limit of the SO2 concentration in the injected stream should be defined with precautions with regard to its corrosive properties. SO2 dissolves in brines and increases the acidity resulting in the dissolution of carbonate minerals. The evolution of porosity and a significant impact on the storage integrity can be expected. The effects of SO2, CO2, and N2 on batch experiments of rock samples from the CO2 pilot Heletz was investigated (Hedayati et al. 2018). Geochemical modeling revealed an increase in porosity due to SO2 reactivity. With regard to N2, and O2, the partition of CO2, N2, and O2 characterized by a N2- and O2-rich front of the plume edge makes them applicable for the monitoring as demonstrated by Wei et al. (2015). In addition, SO2 is a redox reactive species that can transforms into sulfides, sulphate or precipitate in sulfide minerals.


The previous sections highlight efforts to model multiphase multicomponent reactive transport problems and several applications with respect to CCS and potential co-injected impurities. Until now, models required significant simplification for such applications. For example, as mentioned by Wolf et al. (2016), additional brine injections were considered in Toughreact, as impurities were dissolved and the gas phase was pure CO2 (Xu et al. 2007; Xiao et al. 2009; Lei et al. 2016; Todaka and Xu 2017). Consequently, in such simulations, the mass balances, fluxes, and chemical reactions were misrepresented, especially in the vicinity of injection wells due to the artificial addition of water and ions. Wolf et al. (2016) reported on the ability of recent developments to consider trace gas impurities in TOUGHREACT (Xu et al. 2014) with the physical flow properties of a pure CO2 plume with transported trace gas species. In the analysis of MMF&RT methods, a tight SIA-based coupling (such as STOMP) was highlighted. Bacon et al. (2014) conducted simulations of 99% CO2 + 1% H2S injection using STOMP, and the phase equilibrium was modeled using the asymmetric approach with the PR EOS similar to HYTEC (Sin et al. 2017a).

Despite all the intensive studies conducted on convective mixing, few studies considered co-injected impurities; and fewer considered geochemical reactions. For example, Thomas et al. (2016, 2018) recently carried out experiments on convective motion in alkaline and salt solutions. There are at least two major obstacles for this. First, a few simulators can model MMF&RT, which includes multicomponent gas mixtures and their properties, and carry out high-resolution simulations. Second, it is difficult to accurately predict thermodynamic properties, and experimental data for multicomponent systems are indeed critical. The modeling of such complex problems requires solid coupling methods. Recent MMF&RT developments offer large opportunities to model other systems than those previously presented. Several selected applications of MMF&RT and associated numerical difficulties are presented below.

Reservoir modeling is one of the major examples of MMF&RT applications. The modeling of multicomponent gas flow allows for a more in-depth understanding of long-term sour gas reservoir evolution. The knowledge of the compositional distribution within a gas reservoir has economical significance with respect to field development and resources production. Preferential leaching by underlying active aquifers is a significant factor that can cause a heterogeneous gas distribution. Soluble gases such as H2S and CO2 are dissolved and exported by the aquifer, thus developing diffusive gradients and a slow gas density-driven motion, which gradually develops over a long-term evolution. The impact of diffusive transport was studied only by Bonnaud et al (2012). Sin et al. (2017b) extended a problem and demonstrated the importance of modeling real fluids, the solubility of real gas mixtures, gas flows with composition-dependent thermodynamic properties at realistic temperature, and pressure conditions, which are characteristic for deep reservoirs. Based on this physical problem, a benchmark exercise was specifically designed to test multiphase multicomponent reactive transport codes (Sin et al. 2017b).

Carbon geo-sequestration also targets unconventional reservoirs such as shale and coalbed methane formations. Shale gas, which mainly consists of methane, is stored in reservoirs as free gas in fractures and as adsorbed gas on the surface of organic matter, which results in a nanoporous material. The CO2 with its linear molecular structure enters into small pores thus leading to adsorption and further diffusion in the organic matrix. Moreover, CO2 is preferentially adsorbed over CH4 that makes the shale reservoirs suitable for CO2 sequestration while enhancing gas recovery (EGR) (Kang et al. 2011). Godec et al. (2014b) estimated a potential storage of 740 Gt CO2 in shale gas worldwide. A recent pilot study of the injection of 510 t CO2 in the Chattanoga formation was successfully reported by Louk et al. (2017), which demonstrates the feasibility of CO2-EGR. The modeling of multicomponent gas transport in shale formations comprises viscous Darcy flow in micropores and slip flow in nanopores. Competitive adsorption, gas compressibility, and the description of real gas are also required. Numerous concepts were proposed to better understand complex mechanism, such as a new formulation of apparent permeability with a derived definition of tortuosity (Civan 2010), a double porosity approach (Javadpour et al. 2007), and a triple porosity model (Yan et al. 2016).

The IEA Greenhouse Gas R&D Programme (IEAGHG) recently re-assessed the status of CO2-enhanced coalbed methane (CO2-ECBM) and evaluated the worldwide potential storage of nearly 488 GtCO2 in “unmineable coal seams” (Godec et al. 2014a). Despite the lack of geological data on deep coal formation (field tests: USA, Godec et al. 2014a; New Mexico, Siriwardane et al. 2012; Jharia coalfields, India, Vishal et al. 2018), further investigations of the interactions between CO2 and coals are required: in particular, a better understanding and development of comprehensive models of coal swelling/shrinkage and permeability decreases/increases due to gas adsorption and desorption. Intensive experimental and numerical studies on the adsorption of gases and coal were carried out, e.g., the adsorption of pure CO2 and CH4, their mixtures, their influence on the coal matrix and surface area, the verification of analytical models such as the commonly used Langmuir and Toth isotherms (Bae and Bhatia 2006; Ottiger et al. 2008), and the competitive adsorption of CO2 and CH4 mixtures on dry and wet coal (Lee et al. 2013). The injection of CO2 and N2 as a mixture or alternating gases has attracted significant research attention. Moreover, N2 has a lower adsorption capacity than CO2 and CH4. Upon injection into coal, most of the N2 is contained in cleats, whereas desorbed CH4 migrates through the cleats. Given that only a small portion of N2 is adsorbed when compared with CO2, the matrix shrinks, and the permeability and porosity increase. Co-injecting N2 can limit the swelling of the coal matrix and enhance methane recovery. This was confirmed in field studies (Godec et al. 2014a, Oudinot et al. 2017). Stevenson et al. (1991) presented modeling multicomponent adsorption equilibria on the basis of the real adsorbed solution theory, and then provided binary interaction parameters for CO2, N2, and CH4. A recent study indicated a higher desorption of CH4 when injecting a CO2 + SO2 mixture (Luo et al. 2019). Another issue with respect to CO2-ECBM is matrix swelling/shrinkage, which is dependent on the effective stress, in addition to the gas type and gas adsorption/desorption; and therefore the gas compressibility and temperature. Moreover, it is a reversible and anisotropic property. Pan and Connell (2012) presented a comprehensive review of the permeability models and experimental data. In addition to sorption-induced swelling, effective stress, pressure-dependent permeability, and a high gas compressibility, the Klinkenberg effect has a significant influence in the gas flow in cleats (Fan et al. 2019; Wang et al. 2014; Zhu et al. 2007). Moreover, few studies have been conducted on CO2–brine–rock interactions (Wang et al. 2016).

Energy storage via gas storage in saline aquifers or salt caverns (called Power to Gas technology) is a potential application of the type of models presented in this contribution that take into account the influence of salinity on gas solubility. Hemme and van Berk (2017) proposed a numerical model using Phreeqc of CH4 storage in a salt cavern with a complete geochemical approach, whereas the mesh and flow and transport were significantly simplified. With respect to hydrogen storage, we suggest that the Klinkenberg effect for gas motion should be considered.

Within the framework of radioactive waste storage, at least two subjects may require multiphase multicomponent models: the atmospheric carbonation of concrete and the corrosion of metallic components that is responsible for H2 pressure build-up in the storage environment (see also Bildstein et al. 2019, this volume). The first topic is the subject of an ongoing benchmark exercise that combines the drying of concrete (i.e., desaturation), the flow and transport of CO2 from the atmosphere, and the carbonation. The second one has been the subject of a study showing oxic and anoxic corrosion of steel and H2 production, with an extent of oxidizing/reducing front depending on gas diffusion coefficients (De Windt et al. 2014).

The prediction of subsurface soil gas migration may also benefit from improvements in numerical modeling. Providing some developments to consider organic matter degradation and/or respiration from plants and microorganisms, such models could evaluate gas (major gases CO2 and O2, and inter gases as well) fluxes between the soil and atmosphere that could help monitor specific geochemical activities in soils (Alibert et al. 2018).


In this chapter, the existing mathematical and numerical approaches to multiphase multicomponent reactive transport and flow were analyzed, in addition to the necessity for MMF&RT modeling, especially for CO2-related problems.

Over the last two decades, there has been significant progress with respect to the understanding of the physics of multiphase flow. Novel measurement methods for capillary pressure and relative permeability combine experimental techniques and high-resolution numerical modeling. This allows for the accurate estimation of the characteristic properties of homogeneous and heterogeneous cores. Further studies are needed to obtain more experimental data for a wide range of representative rock types, to create a workflow for upscaling from the sub-core and core scales to the m-scale, and to accurately predict capillary trapping and hysteresis necessary for the evaluation of residual trapping during post-injection. Further investigations on relative permeability curves during CO2 exsolution are required, as they are significantly different from those of drainage processes (Falta et al. 2013; Zuo et al. 2017).

Recent developments in modeling thermodynamic systems allow for a better representation of the phase equilibrium. Non-ideal solutions and gas mixtures can currently be considered using (cubic) EOS in conjuction with activity models, e.g., the asymmetric approach. Moreover, the modeling of fluid properties of multicomponent systems is challenging due to the lack of experimental data.

Modeling MMF&RT is a complex problem. The description started from the basics of multiphase flow and geochemical reactions, highlighting the complexity added by the presence of two or more fluid phases. The formulation of multiphase flow is not unique. The decoupled pressure formulation using the global pressure simplifies the resolution of the set of equations, but this approach is suitable only for loosely coupled problems. The method is employed in several couplings with the reactive transport. The flow calculations are then less intensive; however, the iterative loop between the flow and reactive transport is required to cover the gaps between steps. The classic form of conservation equations for two-phase or multiphase flow may exhibit stability issues when selecting the natural primary variables originated by Coats in 1980. To avoid changing the set of primary variables, numerous alternatives were proposed by combining variables, expressing them in different forms, and defining conservation equations with respect to the mole, mass, or volume. The extended set of primary variables with complementarity conditions and the overall composition formulation were tested on CCS related problems and demonstrated high efficiencies in these cases. However, the Coats form is competitive and generally used for coupling with reactive transport.

The GIA is reliable for practical use when applied with efficient linear solvers and reduction techniques. The OSA can be found in most codes due to its flexibility, but retardations of information circulated through the operators can evolve into convergence issues. An obvious remedy is to model implicitly flow and phase equilibrium, iterating between the flow and transport, or solving reactive transport by GIA. It is now clear that the flow and equations of state/phase equilibrium should be solved simultaneously. The stability of MMF&RT is improved for the modeling of coupled processes with an evolving matrix across spatial and temporal scales. Overall, the recommendation is to use a GIA or a tight SIA-based coupling with natural primary variables or equivalent formulation. We consider that the tight SIA is preferable for reactive transport simulators.

0D simulations demonstrated a high accuracy for various systems (e.g., various gaseous compounds, gas mixtures, brine); thus verifying the numerical representation of these systems and allowing for the validation of different approaches selected in the tested geochemical solvers. Additional comparisons with well-constrained experimental results for multi-dimensional problems involving multiphase transport and flow would be of great interest to validate the coupling. A few benchmark studies for MMF&RT were conducted; however, the gas phase was only represented by CO2. Hence, MMF&RT benchmarks with multicomponent systems are still required.

Open access: Article available to all readers online.