Numerical investigations of foam-assisted CO2 storage in saline aquifers

Abstract CO2-foam injection is a promising technology for reducing gas mobility and increasing trapping within the swept region in deep brine aquifers. In this work, a consistent thermodynamic model based on a combination of the Peng-Robinson equation of state (PR EOS) for gas components with an activity model for the aqueous phase is implemented to accurately describe the complex phase-behavior of the CO2-brine system. The phase-behavior module is combined with the representation of foam by an implicit-texture (IT) model with two flow regimes. This combination can accurately capture the complicated dynamics of miscible CO2 foam at various stages of the sequestration process. The Operator-Based Linearization (OBL) approach is applied to improve the efficiency of the highly nonlinear CO2-foam problem by transforming the discretized nonlinear conservation equations into a quasi-linear form based on state-dependent operators. We first validate our simulation results for enhanced CO2 dissolution in a small domain with and without the presence of a capillary transition zone (CTZ). Then a 3D unstructured reservoir is used to examine CO2-foam behavior and its effects on CO2 storage. Simulation studies show good agreement with analytical solutions in both cases with and without CTZ. Besides, the presence of a CTZ enhances the CO2 dissolution rate in brine. Foam simulations show that foams can reduce gas mobility effectively by trapping gas bubbles and inhibit CO2 from migrating upward in the presence of gravity, which in turn improves the sweep efficiency and opens the unswept region for CO2 storage. In the long run (post-injection), with the increasing effects of dissolution, the mechanism of residual trapping, due to the presence of foam, may not be significant. This work suggests a possible strategy to develop an efficient CO2 storage technology.


Introduction
Currently, due to various anthropogenic activities, the concentration of carbon dioxide (CO 2 ) in the atmosphere is having significant and observable effects on the environment. It's believed to be a major contributor to global climate change, such as rising sea level and ocean acidification (IPOC and IPCC, 2014;NASA, 2018). Carbon capture, utilization and storage (CCUS) in subsurface geological formations have been proved to be one viable and promising solution for this environmental issue (Pruess et al., 2004;Raziperchikolaee et al., 2013;Alcorn et al., 2019;Ajoma et al., 2020). Deep saline aquifers have been considered as more ideal sites for CO 2 injection and long-term storage. Compared to other target geological formations, such as depleted oil and gas reservoirs and coal-bed methane, saline aquifers are ubiquitous worldwide and have the largest potential storage capacity, which makes them feasible for large scale long-term sequestration (Gale, 2004;Bachu et al., 2007;Li et al., 2018).
Typically, the presence of an impermeable seal at the top of a formation can hinder CO 2 from moving upward, trapping CO 2 in aquifers (Malik and Islam, 2000;Jessen et al., 2005;Vitoonkijvanich et al., 2015). However, since gas phases generally have higher mobility due to lower viscosity compared to the reservoir fluid, the injected CO 2 will migrate along the top of the reservoir dominated by gravity forces (Hesse et al., 2008). Along this process, CO 2 may leak into the atmosphere if it reaches faults or abandoned wells (Celia and Nordbotten, 2009). This effect also causes very poor sweep efficiency of CO 2 (i.e., lowering storage capacity).
These issues can be overcome or minimized by reducing gas mobility and increasing trapping within the pore space of the swept region. Simultaneous water and gas (SWAG) injection or water alternating gas (WAG) injection can improve CO 2 sweep efficiency (Caudle and Dyes, 1958;Bedrikovetsky, 2003). Laboratory studies have shown that SWAG and WAG injection reduce CO 2 mobility and improve its sweep efficiency. Streamline-based simulation results show co-injection of water at a volume ratio of 15% increases the storage efficiency around 9.0%, compared to 3.0% when only pure gas is injected, while there is a significant improvement of sweep efficiency (Qi et al., 2009).
Foam injection is a promising technology for gas-mobility control in the petroleum industry and aquifer remediation (Rossen, 1996). Recently, the foam enhanced oil recovery (EOR) technique is being extended to CO 2 storage, thus reducing greenhouse gas emissions (Vitoonkijvanich et al., 2015;Izadi and Kam, 2018;Alcorn et al., 2019). Foam is an agglomeration of gas bubbles separated from each other by thin liquid films, which can improve the sweep efficiency of injected gases by mitigating or reducing the effect of low gas viscosity and reservoir layers (Bikerman, 1973;Schramm, 1994;Rossen, 1996). Currently, foam is used in diversion of acid in well-stimulation treatments, diversion of gas in EOR processes and diversion of treatment fluids in soil remediation processes (Rossen and Wang, 1999;Lake et al., 2014). Foam-assisted CO 2 injection (i.e., adding surfactant to generate CO 2 foams in situ) provides insights to maximize the potential of CO 2 storage as well.
Fundamentally, capillary effects and the drag on foam films reduce gas mobility considerably (e.g., by 10 ~ 10 4 times), through trapping gas bubbles (e.g., 90-99% of gas) and increasing the flow resistance of flowing bubbles (Kil et al., 2011). The reduction in gas mobility improves the sweep efficiency remarkably and opens otherwise unswept formation for CO 2 storage. More CO 2 thus is trapped in the pore space rather than migrate upward. The stress on the overburden rock is relaxed, reducing the risk of cracking it. As injection stops, nearly 100% of injected gas in the swept zone is trapped in-situ (as a discontinuous phase) by lamellae (Kil et al., 2011), as long as foam remains stable. The dispersion of CO 2 in liquid increases the contact area of CO 2 with rock and water and thus improves storage capacity.
Prior to foam deployment, one needs to understand the following key issues. The first one is how to predict the behavior of the injected CO 2 stream. In the post-injection period, the footprint of injected CO 2 plays an important role in the security and permanence of CO 2 storage (Li et al., 2018). The key underlying mechanism is how foam can overcome the instability at the interface between the displacing and displaced phases caused by poor mobility ratio (leading to fingering or channeling) and density contrast (leading to gravity segregation). The second important phenomenon is the residual trapping of CO 2 during the migration through the saline aquifer; then dissolution starts to play a significant role at longer timescales. We need an accurate model to represent the major physical and chemical processes induced by CO 2 foam injection into potential disposal reservoirs, such as miscible and immiscible displacement, partitioning of CO 2 among different fluid phases and thermal effects (Pruess et al., 2004). Last, but not least, the nonlinearity of this coupled process challenges conventional simulation, which often translates into an extreme computational cost. It is essential to establish a robust and accurate simulation technique which can model these processes in a realistic and quantitative fashion.
In this work, therefore, we study the coupling of CO 2 sequestration with foam injection (co-injecting CO 2 and surfactant solution). For an accurate description of this phase behavior, a recently developed thermodynamic model based on a combination of a cubic Equation of State (EOS) with an activity model has been implemented (Ziabakhsh-Ganji and Kooi, 2012). This model combines a classic fugacity formulation for the supercritical gas phase and an activity model combined with Henry's law constants for the aqueous brine. This implementation makes the thermodynamic model more accurate than conventional cubic EOS. The implicit-texture (IT) model (CMG-STARS, 2012) used in this study assumes that foam generation and destruction reach a local steady-state instantaneously and represents the effect of foam bubbles implicitly by introducing a mobility-reduction factor. This mobility-reduction factor, used to rescale gas mobility with foam, is a function of water saturation, oil saturation, surfactant concentration, capillary number and salinity.
To accurately simulate these highly nonlinear coupled foam-assisted CO 2 storage processes, a new approach, named Operator-Based Linearization (OBL), where performance, flexibility and robustness can be combined, is introduced to reduce the nonlinearity of complex physical problems (Khait and Voskov, 2018). The OBL approach transforms the discretized mass-conservation equations to space-dependent and state-dependent operators. While space-dependent operators are treated conventionally, the state-dependent operators are approximated by discrete representation on a uniform mesh in parameter-space. These state-dependent operators rely on current local physical properties (e.g., density, viscosity, relative permeability), which represent the most nonlinear part of the governing equations. The continuous representation of these operators is achieved through the multilinear interpolation, which provides a unique tool for approximate representation of the exact physics of the problem. Then the implementation of fully-implicit simulation code is significantly simplified with the OBL methodology. The discretized PDE and the property evaluation are completely separated from each other. That helps to easily implement advanced numerical approaches, e.g., share-memory parallel implementation on CPU or GPU, which can be combined with high flexibility of the simulation code, e.g., the direct implementation of all properties in Python . The OBL approach also provides an opportunity to control the nonlinearity in physics by changing the resolution of parameter space. This paper is structured as follows. First, we briefly describe our numerical and thermodynamic models. Then we validate our simulation capabilities against analytical solutions, mainly focusing on the enhanced CO 2 dissolution. Furthermore, we investigate the behavior of the CO 2 plume with brine-assisted (co-injecting CO 2 and brine) and foam-assisted (co-injecting CO 2 and surfactant solution) CO 2 injection, including the plume footprint, the amount of CO 2 dissolved and residually trapped, storage capacity and efficiency using an unstructured 3D reservoir with homogeneous properties. We conclude the paper by summarizing the main conclusions.

Governing equations
In this section, we briefly consider the governing equations and nonlinear formulation for two-phase multi-component isothermal flow in porous media: where subscript j ∈ {w, n} denotes the wetting phase (brine) and the nonwetting phase (supercritical CO 2 ). ϕ is porosity, s j is phase saturation, ρ j is phase molar density, x cj is component mole fraction in a phase.
u is Darcy velocity, J is Fick's diffusion flux. In addition, the multiphase extension of Darcy's law is applied to describe the flow of the two-phase system: where k is permeability tensor, k rj is relative permeability, μ j is phase viscosity, p j is pressure in phase j, g is the vector of gravitational acceleration, and D is the depth. p c is capillary pressure, which relates the pressures of the two phases. p n is the non-wetting phase, p w is the wetting phase. Capillary pressure is a function of saturation, often expressed as p c (s w ). s j = vj/ρ j ∑ n j j=1 vj/ρ j and v j is the molar fraction of phase j.
J cj is the diffusion-dispersion tensor of component c in phase j, which is described by where x is mass fraction, D is diffusion coefficient. Generally, nearly all foam models alter only the transport properties of gas and assume that liquid properties remain the same function of saturations as in the absence of foam, which is in accordance with laboratory investigations (Friedmann et al., 1991;Rossen, 1996;Dholkawala et al., 2007;Lotfollahi et al., 2016). In the presence of foam, gas is trapped by stationary lemallae to reduce gas mobility. In the implicit-texture foam model (CMG-STARS, 2012) used here, foam reduces gas mobility by modifying gas relative permeability with a mobility-reduction factor (FM) as shown below: where k f rg and k rg are gas relative permeability with and without foam, respectively; fmmob is defined as the maximum-attainable gas-mobility reduction, and F 1 through F 6 are functions accounting for the effects of physical factors on gas mobility (e.g., surfactant concentration, water saturation, oil saturation, oil composition, capillary number, and salinity). In this project, we consider only two functions, F 1 and F 5 , capturing the effects of surfactant concentration and water saturation on foam strength. The details are shown in the appendix.
A finite-volume discretization on a general structured mesh approximated in space using a two-point flux approximation (TPFA) and backward Euler approximation in time is applied. This introduces strong nonlinearity into the system of the governing equations, especially in the presence of complicated physics. We need to linearize the problem, which requires determining all the partial derivatives with respect to these nonlinear unknowns and assembling the Jacobian and residuals. After the linearization step, the Newton-Raphson method is adopted to solve the linearized system of equations on each nonlinear iteration. In conventional simulation, the Jacobian should be assembled with accurate numerical property values and their derivatives with respect to nonlinear unknowns. This process requires either various interpolations (for properties such as relative permeabilities of different phases) or solution of a highly nonlinear system in combination with the chain rule and inverse theorem, which could increase the computational cost.

OBL approach
Following the OBL approach, all variables in the Eq. (1) fully defined by the physical state ω can be grouped together and represented by the state-dependent operators (Khait and Voskov, 2018). The discretized mass-conservation equation in operator form is Here, where ω and ω n are nonlinear unknowns in the current and previous timestep, respectively; L(i) is the set of neighbors of the control volume l; θ(ξ, ω, u) is the source term. V, ϕ 0 , and c r are initial volume, porosity and rock compressibility, respectively. Φ l j is the phase pressure difference between neighbor cells. Γ l and Γ l d are the space-dependent part of convective and diffusive transmissibility, respectively. In this study, both effects of gravity and capillarity are considered. The phasepotential-upwinding (PPU) strategy is applied to compute the numerical flux (Khait and Voskov, 2018).
This representation allows us to decouple a computation of nonlinear physics from conventional discretization terms. Instead of performing complex evaluations of properties and their derivatives with respect to nonlinear unknowns in the course of simulation, we can parameterize operators in physical space at the preprocessing stage or adaptively with a limited number of supporting points (Khait and Voskov, 2018). Then during the simulation, a multi-linear interpolation is applied to evaluate the operators in the current timestep, which improves the performance of the linearization stage. Meanwhile, this approach can reduce the nonlinearity of the physical problem due to the application of piece-wise representation of operators (Voskov, 2017).

Thermodynamic model
The model describes thermodynamic equilibrium between a nonaqueous phase (i.e., a multi-component mixture which can be in gas, supercritical or condensed conditions) and an aqueous phase (i.e., liquid which includes dissolved hydrocarbon and gases). Due to the instantaneous local equilibrium assumption, phase-behaviour calculations are decoupled from flow and transport. In a multi-phase system, an exact thermodynamic equilibrium is required at every nonlinear iteration in the molar formulation Here z c = ∑ j x cj p j s j / ∑ j p j s j is overall composition and f cj (p, T, x j ) is the fugacity of component c in phase j. The set of thermodynamic relations described by Eq. (14) to Eq. (17) must be simultaneously solved for the conditions of pressure, temperature and composition in each grid block in the nonlinear loop.
In this work, a fugacity-activity model is used to solve for thermodynamic equilibrium based on the idea originally proposed by Kritchevsky and Iliinskaya (1945). In this approach, the fugacity of the gas phase is expressed in terms of the fugacity coefficient (f g c = pψ c y c ) and the aqueous phase in terms of activity (f w c = h c κ c x c ). In thermodynamic equilibrium (f g c = f w c ), the phase-equilibrium constant of each component K c can be obtained: where p is the total pressure in the system, ψ c the fugacity coefficient of the gas phase, h c Henry's constant, κ c activity coefficient, x c and y c the molar fraction of each component in aqueous phase and gas phase, respectively. Eq. (18) is used to calculate K values for different gas components. The equilibrium constant for the water component is calculated with a separate relation proposed by Spycher et al. (2003): where K 0 H2O is the equilibrium constant of H 2 O at the reference pressure of 1 bar, T is temperature in Kelvins, V H2O molar volume of H 2 O. More detailed description can be found in Spycher et al. (2003).
Phase calculations are performed on all phases and phase partitioning is calculated using multistage negative flash as described by Iranshahr et al. (2010) with successive substitution iteration. In order to initiate the negative-flash procedure, composition-independent ideal K-values provide an initial guess of phase fractions. Then, based on the output of the first iteration (phase fractions and composition of each phase), fugacity coefficients are updated to obtain new K-values. Once the thermodynamic system is solved, the thermophysical properties associated with the mass-conservation equations, such as phase density and phase viscosity, can be determined. The accuracy of this thermodynamic model vs. experimental results has been validated in Morshuis (2019).

Enhanced dissolution
We begin by validating our simulation approach through studying the detailed behavior of gravity induced instabilities and the associated dissolution rate in small domains. Elenius et al. (2012Elenius et al. ( , 2014 investigated the full problem of two-phase flow with gravity currents and convective dissolution in the absence and presence of the capillary transition zone (CTZ), and these results can be used as a benchmark for verification of our simulation approach. In this work, we take two small models, as shown in Fig. 1. One represents a scenario where the CTZ is negligible, and another one is with a realistic capillary transition zone. All the parameters which are used in the simulations and the simplifications in these models can be found in Elenius et al. (2015). Fig. 2 demonstrates the CO 2 concentration for the simulation after 200 years with single-phase brine and the simulation with a stagnant CTZ. Obviously, the stagnant CTZ enhances the concentration of CO 2 in the fingers away from the interface, leading to a faster propagation of the fingers, compared with the no-flux top boundary case. This is consistent with the findings of Elenius et al. (2012Elenius et al. ( , 2014. Therefore, we can infer that the presence of the stagnant CTZ, to some extent, can improve the storage efficiency by enhancing dissolution rate. Following the definition of dissolution rate in Elenius et al. (2015), we calculate the rate of CO 2 mass transfer to the (single-phase) brine region across the interface per area (length) of the top interface: where h and c are the thickness and mean concentration of the singlephase brine region respectively. Elenius et al. (2014) also provided a semi-analytical solution for the dissolution rate with the effect of the capillary transition zone: and at negligible effect of the transition zone: where K is permeability, Δρ w density difference between brine and brine with dissolved CO 2 , g gravitational acceleration, X max maximum solubility, μ w water viscosity, and d the exponent of the relative-permeability function which is obtained by fitting the water relative permeability.  Initial position of region with brine (blue, X = 0 kg/kg) and two-phase conditions (red, X = 0.03 kg/kg, corresponding to x = 0.0125 mol CO 2 /mol brine). In (a), nowflow conditions are applied for all boundaries, and the concentration and pressure are fixed at the top of the domain by specifying a large pore volume; in (b), CO 2 is provided by means of the CTZ, but the entire two-phase region has a very large pore volume to maintain the initial saturation profile and the high CO 2 concentration. For further details, see Elenius et al. (2015). the dissolution rate is reduced with time until the nonlinear onset time is reached. It also shows that the presence of CTZ can reduce the onset time. After the nonlinear onset time, fingers start growing and the rate increases due to convection. For both the single-phase and the two-phase with a CTZ simulations, the dissolution rate stabilizes close to the analytical solution.
As shown in Fig. 3, the dissolution is reduced at late times when CO 2 fingers approach the bottom of the aquifer. Here we use the stagnant CTZ to investigate the behavior of fingers at late times. CO 2 starts to dissolve in brine and fill up the domain gradually (Fig. 4a). But the dissolution rate is reduced at late time mainly because of the merging of fingers and the increase of overall CO 2 concentration. After 3000 years, CO 2 concentration is already rather high, though it is still below the solubility limit anywhere in the single-phase brine region. Our simulation results with the CTZ show a similar t peel = 350 years, i.e., the time at which the dissolution rate starts to decrease, which is consistent with Slim (2014)'s findings. After t peel , Slim also found the dissolution decreases from a constant value to a value proportional to 1/(t + g) 2 without a CTZ (1/t 2 in Elenius et al. (2015)). Here we fit the coefficient g based on our simulation results with the CTZ, and g ≈ − 1000 gives a good match (Fig. 4b). These results validate the greater accuracy of our enhanced dissolution model.

Model description
When CO 2 is injected into a formation saturated with brine, it migrates upwards due to gravity and forms a nearly horizontal layer overlying the brine phase. After a short time, CO 2 starts to dissolve in the brine, as a result of molecular diffusion and density-driven convection and in part is trapped in situ as residual gas. Many researchers have found that foam-assisted CO 2 injection can increase sweep efficiency by mitigating gravity segregation processes (Vitoonkijvanich et al., 2015;Izadi and Kam, 2018). Therefore, it can increase the storage capacity due to the larger swept area and the increasing residual gas saturation.
In order to simulate this process, we consider a 3D homogeneous horizontal reservoir with unstructured mesh and fine mesh size as shown    21) and (22) in the presence and absence of the CTZ. The subfigure inside shows the mass flux at early times.
in Fig. 5. The height and the radius of the model is 30 m and 400 m, respectively. There are 30 layers and the average number of elements in radial direction is 192. The top and the bottom surfaces of the reservoir are no-flow boundaries. We also assume for simplicity that surfactant is already present in the water phase throughout the porous medium and the adsorption of surfactant is neglected. Other parameters, such as rock and fluid properties, are listed in Table 1. Although the scale of this model is just a few hundred meters, it provides an accurate representation of CO 2 sequestration with realistic thermodynamics conditions. As shown in Elenius et al. (2015), the proposed mesh resolution (around meter scale) provides a numerically converged solution for enhanced dissolution phenomena, which is studied here in a fully 3D setting.
To simplify the problem, we neglect any chemical reactions imposed in the brine by interactions with CO 2 phase, such as CO 2 -rock mineral reactions and CO 2 -brine dissociation. The temperature assumed to be constant during the simulation. The simulation domain, a 5 • sector of the cylinder, is initially saturated with formation brine with no dissolved CO 2 . The injection well fully perforating the entire vertical interval is located at the left boundary and constant pressure is assumed at right boundary with no-flow conditions along the rest of interfaces. A fixed gas injection rate of 4.0 m 3 /day, corresponds to 0.06 Mt/year for full domain (normalized to 360 • ). The injection well is closed after one year of injection.
Another simplification is the model of gas trapping due to the presence of foam. Gas trapping is an important mechanism in the foamassisted CO 2 storage process, especially after injection. Friedmann et al. (1991) measured trapped gas fractions in the range 75% to 90% over a wide range of velocities. Tang and Kovscek (2006) found a significant decrease in trapped gas with increasing gas velocity. Jones et al. (2018) also found in micro-models that as the superficial velocity increases, the fraction of trapping gas decreases. There are no complete models to describe the amount of trapped gas due to the injection of foam. In our study, for simplicity, we assume the residual (i.e., trapped) gas saturation rises by 0.1 in the presence of foam. This assumption is not rigorously correct because, as noted, the trapped gas saturation with foam is larger. Such low value, to some extent, can represent a reduction in gas trapping due to depletion of surfactant in long term. In addition, in the upper layer where foam is collapsed or cannot be generated, the residual saturation does not change. During the simulation, only one set of relative-permeability curve is used. However, gas saturation is much larger than S gr and the only effect of this assumption is a modest reduction in k rg .
Foam-assisted CO 2 storage simulations for a brine aquifer are performed with the Delft Advanced Research Terra Simulator (DARTS) which is capable of modeling complex flow and transport related to various energy applications (Khait and Voskov, 2017;Kala and Voskov, 2020;Wang et al., 2020). A combination of Peng-Robinson (Peng and Robinson, 1976) and Kritchevsky-Illiinskaya (Kritchevsky and Iliinskaya, 1945) equations of state is deployed in this study because it could provide more reasonable results for the vapor-liquid equilibrium properties as well as the volumetric properties of CO 2 mixtures (Li and Yan, 2009). The empirical correlation used to determine the brine solution density was developed by Spivey et al. (2004). Garcia (2001) provided a correlation for the density of brine with dissolved CO 2 . The aqueous viscosity is computed by the correlations developed by Mao and Duan (2009) (brine solution) and Islam and Carlson (2012) (brine with dissolved CO 2 ). The density and viscosity of non-aqueous phases are determined by Peng and Robinson (1976) and Lee et al. (1966), respectively.

Results and discussion
In this section, we present the results of brine-assisted and foamassisted CO 2 injection into a homogeneous reservoir, including the behavior of the CO 2 plume in injection and subsequent post-injection processes.
The injected CO 2 exists as supercritical fluid under the selected reservoir conditions. Fig. 6 illustrates the saturation of the supercritical CO 2 after 1 year injection. During the brine-assisted CO 2 injection, gas segregates with water and migrates upwards quickly because of the low density and viscosity of CO 2 compared with the formation brine. In the meantime, it displaces the formation brine and thereby increases the contact area for CO 2 storage. The plume, however, sweeps only the nearwell region and then rises to the upper layer. Thus the storage efficiency, especially in the near-well region, is rather low due to the limited swept region.
Foam injection can significantly enlarge the swept area by reducing gas mobility. When CO 2 and surfactant are co-injected into the formation, foam can be generated in the near-well region; then gas mobility is reduced remarkably (max. 100 times in this study) and much more space will be open for CO 2 storage, see Fig. 6(b) for details. The plume front in foam injection moves slowly and uniformly, which reduces the risk of leakage, especially during CO 2 EOR processes where wells distance is  Fig. 6. Saturation of supercritical CO 2 after 1 year injection. The white dashed line is the CO 2 plume front.
X. Lyu et al. limited. Under steady-state, an analytical model for uniform co-injection of water and gas in homogeneous, horizontal reservoirs can be used to predict the segregation length (Stone, 1982). In this study, less than 10% pore volume (0.06 PV) of gas is injected. No obviously separated regions, therefore, can be distinguished with a sharp boundary compared with the previous research (Stone, 1982;Rossen et al., 2010). However, in this transient displacement process, foam exhibits its capacity to hinder gas rising upwards and increase the sweep area. Fig. 6(b) shows that the segregation point where water and gas separate completely, is more than 100 m from the injection well. At early time, foam may reduce the dissolution rate due to the reduced contact area between CO 2 and brine in the upper layers. However, in the long run, the dissolution increases because the free gas after segregation as well as collapsed foam still migrates upwards to overlie the brine phase in the upper layer, thus increasing the contact area. With a fixed injection rate, the required injection pressure for foam is much higher, around 125.4 bar; while the injection pressure is only 93.8 bar for co-injecting water and gas.
Figs. 7 and 8 display the saturation of the supercritical CO 2 with time. In both cases, the mobile CO 2 forms a nearly horizontal layer overlying the brine phase. As shown in Fig. 6(a), when injection ceases, the front of CO 2 plume approaches the right boundary. Therefore, the CO 2 plume arrives at the right boundary in a short time in the post-injection period. With the dissolution of CO 2 in the upper part of reservoir, the leading tip retracts and disappears gradually (Figs. 7(a) and 8 (a)). After foam injection, gravitational force dominates the flow, and gas migrates upwards and accumulates there. Once gas saturation is high enough (i.e., water saturation is lower than the limiting water saturation) in the upper layer, foam collapses and gas mobility increases dramatically. Foam cannot be re-generated there, which makes the override zone thin in the foam-assisted post-injection process (Figs. 7(b) and 8 (b)). Foam-injection retards the late-time dissolution rate. However, the residual trapped CO 2 phase with foam-assisted injection is much greater than that of brine-assisted injection, in terms of the swept area and saturation of immobile gas. Foam increases the swept area and during the post-injection process, the residual gas saturation increases through foam trapping gas bubbles. The enlarged swept area provides higher capacity for trapping of CO 2 .
In our simulation of one year of injection, there is no override zone ahead of the foam zone until gas injection ceases. At this time gas migrates upward from the foam zone and forms an override zone that extends radially outward. Over time, as the override zone grows, gas saturation within that zone falls to residual gas saturation. Below the override zone (in dark red in Fig. 7(b)), there is a second zone (two grid blocks deep) with residual gas. This zone is created during the advance of the override zone, due to lower mobility of gas at intermediate gas  saturations. This effect is magnified by numerical dispersion at the displacement front (Lyu et al., 2021). Later, residual gas in both zones can dissolve into water connected to the top of the aquifer, much as in the capillary transition zone in Section 3. Residual CO 2 in both override zones dissolves into brine gradually over time, as shown in Fig. 7(a). There is also large zone of trapped residual CO 2 near the well, where foam remains stable (i.e., at lower water saturation). In practice, one could increase the injection pressure to expand the swept area (Rossen et al., 2010), subject to limitations on injection pressure.
Figs. 9 and 10 show the mole-fraction distribution of CO 2 with time. CO 2 fingers move downwards and grow gradually in both cases. The fingers between the override zone and bottom brine form earlier in brine-assisted CO 2 injection because override happens rapidly (Figs. 9 (a) and 10 (a)). Finally, the average CO 2 concentration in the whole domain (excluding the residual trapped region) in brine-assisted injection is higher than that with foam-injection. Once the tips of fingers reach the bottom boundary of the domain, CO 2 fingers start to expand in the horizontal direction and merge with others. The number of fingers therefore is reduced, resulting from the mutual interaction between the fingers during the diffusion process. Note that the brine-assisted and foam-assisted CO 2 injection shows similar behavior, including the migration and dissolution of the CO 2 plume. The injection of foam is mainly applied to prevent CO 2 from migrating upwards and reduce the breakthrough time during the injection period: the effects of foam on CO 2 plume migration and dissolution at the upper layers at later time are negligible.
In order to observe how the leading tip propagation changes with time, we show the results in foam-assisted injection (Fig. 11) where the leading tip stops before it reaches the right boundary. As mentioned above, foam does not affect the migration of the CO 2 override zone, so  this result can represent the behavior of the CO 2 plume for either brineassisted or foam-assisted CO 2 co-injection strategies in the post-injection period as long as the domain is large enough. The plume speed decreases with time until the plume stops and retracts after about approximately 150 years, 370 m away from the injection point. The presence of the CTZ causes a reduction in tip speed. Our results show a similar trend to those of Elenius et al. (2015). This interaction between the speed of the leading tip and convective mixing also can be observed from the distribution of dissolved CO 2 under the plume, see Fig. 8. Fig. 12 displays the global mass transfer into the single-phase brine region, which is defined as the amount of CO 2 entering the single-phase region per unit time: R = dM CO2 /dt. Both injection strategies show similar results: R increases at early time and later decreases with time. As shown in Elenius et al. (2015), the global mass-transfer decreases gradually after the t peel , which is different from our simulation results. In our simulation, the thickness of the domain is just 30 m, which causes the fingers reaching the bottom boundary in a very short time (around 150 years). Once the fingers arrive at the bottom, the dissolution rate starts to decrease, also seen in Fig. 3. With foam injection, R increases faster at early time and reaches a slightly lower peak. On the one hand, once the injection ceases, foam sweeps much more area, increasing CO 2 trapping, leading to a higher dissolution rate over a short period. On the other hand, the increased residual gas reduces the amount the CO 2 which can dissolve into brine.
In this work, all properties are dependent on pressure, temperature and molar composition of each component. Therefore, Eqs. (21) and (22), are not necessarily valid. However, in the post-injection process, the variation of pressure is slight (~3 bar), and we assume constant temperature. Therefore we still can use Eq. (21) to approximate the enhanced dissolution rate due to the presence of the CTZ. Note all the properties in Eq. (21) are average: for instance, we calculate all water densities in all elements of the mesh and divide it by the total number of elements to get the corresponding water density. Here, Δρ w = 5.75 kg/ m 3 , X max = 0.017 mol/mol (0.0415 kg/kg), ρ w (X max ) = 982.6 kg/m 3 , and μ w = 0.86 cp. We then obtain the average dissolution rate with the CTZ, F ave = 0.254 kg/(m 3 year) (Eq. (21)). We compare this analytical dissolution rate with our simulation results. In brine-assisted CO 2 injection, R max = 680 kg/year, corresponding to F max = 0.325 kg/(m 3 year) (F max = R max /(A × ϕ)). This dissolution rate is 27.9% larger than that of analytical solution.
Considering the trapping mechanisms and time scale in this research, we estimate the effectiveness of CO 2 geological storage, and three trapping indices are used to represent the contribution of residual trapping and dissolution trapping mechanism,

Residual trapping index (RTI) =
Total mass of residually trapped CO 2 (kg) Total mass of injected CO 2 (kg) , Dissolution trapping index (DTI) = Total mass of dissolved CO 2 (kg) Total mass of injected CO 2 (kg) , Total trapping index (TTI) = RTI + DTI. Fig. 13 shows the variation of the trapping indices of different injection strategies over time. The CO 2 plume moves further from the well and enlarges the contact area between the plume and formation brine after shutting off the well. Thus enables much more efficient dissolution of CO 2 into the aqueous phase at the two-phase interface; DTI increases accordingly. The capacity for dissolving CO 2 in brine-assisted CO 2 injection is much greater while the amount of residually trapped CO 2 is lower. The variation of RTI is opposite to that of DTI in both cases and less significant in brine-assisted CO 2 injection. However, residual trapping plays a more important role in foam-assisted injection, with a greater trapping index (0.32). After 1000 years, around 92.5% of CO 2 is dissolved into brine after co-injecting brine and CO 2 compared to 62.3% of dissolved CO 2 with foam-injection. In total, 94.3% of CO 2 is trapped in foam-assisted CO 2 injection, increased by around 1.5% compared to brine-assisted CO 2 injection. The efficiency of CO 2 storage, expressed by the ratio of the volume of CO 2 accessible or occupied by CO 2 in a given pore volume of a porous medium to that volume, is different in these two scenarios, though the total trapping index is close. The storage efficiency of foam-assisted CO 2 injection is about 23.4% which is around 8 times than that of brine-water co-injection (3.0%), due to the enlarged swept area by foam.
As mentioned above, foam can mitigate gravity override during CO 2 injection and reduce the risk of leakage or breakthrough. At early time, foam can improve the amount of trapped CO 2 , but in the long run, with   the increasing ability of dissolution, the mechanism of residual trapping may play a less-important role. More-accurate modeling is required to predict the foam characteristics in CO 2 storage processes.
In this study, we use a simple foam model to investigate the effect of foam co-injection to CO 2 trapping. This model does not capture all the characteristics, but it still represents some of the most important mechanisms of foam-assisted CO 2 injection. For practical applications, foam generation and coalescence should be included into the physical model, and gas trapping should be represented more completely. There are other essential issues, such as the cost of surfactant, the depletion of surfactant over time, and the foam injectivity, to be considered. These factors will be taken into account in the future research.

Conclusions
In this work, we develop and validate a realistic phase-behavior model for simulation of CO 2 sequestration in aquifers. The consistent thermodynamic model, based on a combination of a classic cubic equation of state (EOS) for gas components with an activity model for the aqueous phase, can accurately predict the complex phase behavior of the CO 2 plume in brine. An advanced numerical performance provided by the Operator-Based Linearization scheme allows us to perform fullphysics simulation in a 3D sector model. The CO 2 sequestration physics is complemented with a foam model which provides us the ability to investigate the effect of foam co-injection on CO 2 trapping. The following conclusions can be made: • The dissolution rate caused by the gravitational instabilities is enhanced further in the presence of a capillary transition zone (CTZ). Our numerical results show good agreement with the analytical solution in the simplified 2D setting.
• Foam injection can mitigate gravity override during gas injection by reducing gas mobility. This process increases the amount of residual trapped CO 2 by 32.0% in this study. In addition, the presence of foam reduces the amount of flowing gas, thus reducing the risk of leakage. With a more realistic treatment of dissolution fingers in 3D model, the predicted average dissolution rate is almost 30% larger than that predicted by the analytical model. • The final total trapping index in both cases are close in 1000 years, indicating that in the long run (post-injection), with the increasing ability of dissolution, the mechanism of increased residual trapping, due to the presence of foam, may not be significant.

Declaration of Competing Interest
The authors report no declarations of interest.
The widely used implicit-texture foam model, also named CMG-STARS model (Cheng et al., 2000;CMG-STARS, 2012), is applied to investigate the effects of water saturation (S w ) and surfactant concentration (W s ) on foam stability. For simplicity, here we list only those factors used in these calculations and simulations:  where k 0 rg (S w ) is gas relative permeability without foam. fmmob, fmdry, epdry, fmsurf, and epsurf are foam model parameters. As shown in Eq. (A.1), gas mobility is reduced in the presence of foam by scaling foam-free gas relative permeability in this model, as illustrated in Fig. A.1. Foam can form whenever water, gas, and surfactant meet in sufficient quantities.
The Brooks-Corey relations for relative-permeability model are as follows: We take the entry pressure to be p e = 0.2 bar. C is equal to 0.0109 to exclude hysteresis in our simulation. The detailed calculations of phase density and viscosity can be found in Morshuis (2019). The model can accurately represent phase properties within the given pressure and temperature.