Numerical Study of the Injection of Carbon Dioxide in a Homogeneous Porous Media

This work proposes a locally conservative and less restrictive algorithm to solve the problem dealt with in [1], i.e. a two-phase flow in a homogeneous porous medium (water and CO2), with mass absorption between the fluid phases and reaction between the CO2 phase and the rock. The latter is modeled by two non-linear hyperbolic equations that represent the transport of the flowing phases for a given velocity field (equations of saturation and concentration). From the numerical point of view, we use the operator splitting technique to properly treat the time scale of each physical phenomenon and a high-order non-oscillatory central-scheme finite volume method for nonlinear hyperbolic equations proposed by [2] that was extended for a system of equations with source terms to treat the equations that govern the saturation and concentration of phases. In addition, with respect to source terms, the mass flux between fluid phases was handled using the flash methodology, whereas kinetic theory was applied for reproducing the changes in porosity and permeability that are caused by the reaction of CO2 with the rock. The same physical trends observed in [1] were obtained in our numerical results which indicate a good predictive capability. Furthermore, this method avoids the difficulties that arise when adopting small time steps enforced by CFL stability restrictions. Finally, the results obtained show that the applicability of the KT method is beyond just a single nonlinear conservation law with the absence of source terms.

1. Introduction

It is known to all governments that the intensification of the greenhouse effect, due to the emission of gases which arises from industrial activities, is increasing the global temperature every year, dramatically affecting the Earth’s climate (see Houghton et al. [3] ). According to Ravagnani and Suslick II [4] among all polluting gases, it is carbon dioxide (CO2) the one that contributes the most to the intensification of the greenhouse effect. One of the possible ways to mitigate CO2 effects on climate change is called carbon sequestration, which means separation, transportation (supercritical state) and storage of CO2 in permeable rocks.

Models of geological CO2 sequestration should consider multiphase flows, equilibrium and non-equilibrium reactions, and partition phases in reservoirs. This is a prominent problem that has been investigated by several researchers: Meer [5] [6] [7] , Holt et al. [8] , Weir et al. [9] , Law and Bachu [10] , Lindberg [11] , Johnson et al. [12] , Ennis-King and Paterson [13] , Wellman et al. [14] , Pruess et al. [15] , Xu and Pruess [16] and Kumar et al. [17] . All these works address the CO2 sequestration problem but ignore the reaction of CO2 with the rock. Going into more details on the preceding references, Malik and Islam [18] conducted a study in which CO2 collected was injected into the Weyburn field. The purpose was to determine the optimal point of injection to maximize CO2 storage during an oil recovery procedure. Brinks and Fanchi [19] conducted a study of the coupling between CO2 motion and seismic response. Pruess et al. [15] and Xu and Pruess [16] used a detailed geochemical model to study the interactions between the fluid and the rock during a radial flow. Kumar et al. [17] conducted a thorough analysis of CO2 storage using a commercial simulator in a mesh of 64,000 elements. The authors concluded that the trapping of CO2 was an important mechanism for permanent storage. Ennis-King and Paterson [13] used a radial geometry to study in detail the factors that affect the geological sequestration of CO2, including the gravitational effect. Other authors considered the roles of the hysteresis of relative permeability (Spiteri et al. [20] ), dispersion (Calabrese et al. [21] ) and convective mixing (Ennis-King and Paterson [22] ) during storage of CO2. Finally, Obi and Blunt [1] applied the streamline method to solve the CO2 storage in an aquifer in the North Sea. They concluded that the CO2 distribution is dominated by advective transport due to the multiphase flow, and CO2 moves preferentially through the high permeability paths. There is no reaction; the flow of CO2 occurs until the value of residual saturation is reached. The precipitation leads to a decrease in porosity and permeability, while CO2 is stored in the solid phase. The storage efficiency of this mechanism is low, around 2%.

We present a locally conservative numerical methodology to simulate the two-phase flow (water and CO2) with mass absorption between the fluid phases and reaction between the CO2 phase and rock in a homogeneous porous medium. The transport of the flowing phases is modeled by two non-linear hyperbolic equations (equations of saturation and concentration). From the numerical point of view, we use the operator-splitting technique to properly handle the time scale of each physical phenomenon and a high-order non-oscillatory central-scheme finite volume method for nonlinear hyperbolic equations that govern the saturation and concentration of phases. Our work innovates by employing the KT scheme developed by Kurganov and Tadmor [2] instead of the streamline method used by Obi and Blunt [1] . Furthermore, the method was expanded for a system of equations with source terms in which the time scale was different for each physical phenomenon. This method avoids the difficulties that arise when adopting small time steps enforced by CFL stability restrictions (see Nessyahu and Tadmor [23] ). According to Kurganov and Lin [24] and Correa and Borges [25] , it happens because the underlying concept is to explore the local propagation speeds to obtain a more precise estimate of the width of Riemann fans, see Leveque [26] , and evolves the solution on a non-uniform dual mesh based on symmetrical non-smooth control volumes that include the Riemann fans. This procedure led to a substantial reduction of the numerical dissipation and, after it was projected back onto the original mesh, gave rise to an improved fully-discrete central scheme that exhibits a concise semi-discrete form. In addition, with respect to source terms, the mass flux between fluid phases was treated1 using the flash methodology (see [27] ) and the reaction of CO2 with rock2 which causes changes in porosity and permeability, was treated by applying principles of kinetic theory (see [14] ).

2. Mathematical Model

1The dissolution of CO2 in the aqueous phase.

2Which is known as precipitation.

Problem: Given the Darcy velocity ${v}_{Dt}:\Omega ×\left[0,T\right]\to {R}^{n},n=1,2,3$ . The porosity $\varphi :\Omega \to \left(0,1\right)$ , saturation $S:\Omega \to \left(0,1\right)$ and concentration of the CO2 phase $C:\Omega \to R$ are dictated by the following equations

$\frac{\partial \left({S}_{c}\varphi \right)}{\partial t}+\nabla \cdot \left({f}_{c}{v}_{Dt}\right)=-{T}_{d}^{c}$ (1)

$\frac{\partial \left(C\left(1-{S}_{c}\right)\varphi \right)}{\partial t}+\nabla \cdot \left[C\left(1-{f}_{c}\right){v}_{Dt}\right]={T}_{d}^{c}-{T}_{r}$ (2)

$\frac{\partial \varphi }{\partial t}=-\left(1-\frac{{m}_{s}}{{m}_{o}}\right){T}_{r}\text{.}$ (3)

Subject to boundary and initial conditions:

$\left\{\begin{array}{l}\varphi =\stackrel{¯}{\varphi },\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{on}\text{\hspace{0.17em}}{\Gamma }_{D}^{\varphi }\\ S\left(x\right)=s\left(x\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{on}\text{\hspace{0.17em}}{\Gamma }_{D}^{S}\\ C\left(x\right)=c\left(x\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{on}\text{\hspace{0.17em}}{\Gamma }_{D}^{C}\\ p\left(x,0\right)={p}_{0}\left(x\right)\\ \varphi \left(x,0\right)={\varphi }_{0}\left(x\right)\\ S\left(x,0\right)={S}_{0}\left(x\right)\\ C\left(x,0\right)={C}_{0}\left(x\right)\end{array}$

where ${f}_{c}$ is the fractional flow function. In order to understand this concept, let us start defining the total mobility as a function of the relative permeabilities of the fluid phases,

${\lambda }_{t}={\lambda }_{w}+{\lambda }_{c}=\frac{{K}_{Rw}}{{\mu }_{w}}+\frac{{K}_{Rc}}{{\mu }_{c}}$ (4)

where ${K}_{Rw}$ and ${K}_{Rc}$ are the coefficients of relative permeability of water and CO2, respectively, which vary from 0 to 1 as functions of the saturation. These parameters quantify the additional resistance that a moving fluid exerts on another. Furthermore, the fractional flow function of phase CO2, which is the ratio of the mobility of the phase CO2 and total mobility, is given by

${f}_{c}=\frac{{\lambda }_{c}}{{\lambda }_{t}}$ (5)

(see for further details [28] ).

It is worth noting the physical meaning behind the govern equations. The first one represents the transport of CO2 in its own phase, i.e., just modeling the transport of the single dissolved specie. The second one express the transport of the dissolved CO2 in the aqueous phase. Besides that, source terms that embodies dissolution (this source term appear at the first equation also) and precipitation phenomena can be found. Finally, the latter equation has implicitly the following hypothesis: (i instantaneous dissolution which systematically means that $C={K}_{d}$ if $C>0$ and ${T}_{d}^{c}=0$ if $S=0$ , (ii the precipitation has a forward a preferential direction, i.e.,the precipitation of the second mineral just occurs (if and only if) the precipitation of the primary mineral has happened, specifically, ${\text{CO}}_{2}+\text{XO}\to {\text{XCO}}_{\text{3}}$ , where XO is the primary mineral and XCO3 the second mineral. and (iii) the reaction rate is assumed to be proportional to the porosity as follows

$\frac{\partial \varphi }{\partial t}=\left\{\begin{array}{l}-\left(1-\frac{{m}_{s}}{{m}_{o}}\right){T}_{r}=-\psi \varphi \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{if}\text{\hspace{0.17em}}C>0\\ 0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{if}\text{\hspace{0.17em}}C=0\end{array}$

where y is a rate constant.

3. Algorithm

In the following subsections we are going to show a detailed description of the algorithm. As mentioned in the introduction, this method avoids the difficulties that arise when adopting small time steps enforced by CFL stability restrictions.

3.1. Evolution in Time

We begin with the description of the temporal evolution algorithm with fractional steps associated with the decomposition of operators for solving the system of equations of the proposed model. The given velocity field is used in the equations of convective transport that update the variables Sc e C. In this stage of the algorithm, the time step $\Delta {t}_{c}$ restricted by the CFL condition (see Kurganov and Tadmor [2] , Nessyahu and Tadmor [23] ). To simplify the algorithm, we make $\Delta {t}_{d}=\Delta {t}_{c}$ and execute the dissolution step with the new saturations and concentrations computed. After updating the variables Sc and C over the convective and the dissolution microstepping, we update the permeability and the porosity fields during the reactive step until the time step $\Delta {t}_{r}$ takes place (see Figure 1).

3.2. Predictor Problem

This procedure will be used for both the saturation and the concentration problems. However, for the sake of simplicity, the following development will be applied for the saturation S, which does not affect the generality of the method.

Given the porosity field $\varphi \left(x\right)$ and Darcy total velocity ${v}_{Dt}$ at time ${t}^{n}$ , find the saturation S such that:

$\varphi \frac{\partial S}{\partial t}+\nabla \cdot \left(f\right)=0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{em}\text{\hspace{0.17em}}\Omega ×\left({t}^{n},{t}^{n+1}\right)$ (6)

where $f\left(S\right)=\left[\frac{f\left(S\right)}{g\left(S\right)}\right]$ and satisfying the initial condition $S\left(x,{t}^{n}\right)={S}^{n}$ .

The KT method proposed by Kurganov and Tadmor [2] was developed for the nonlinear conservation laws, given by Equation (6), to a field of constant porosity f. The extension of the formulation for the case where $\varphi =\varphi \left(x\right)$ , i.e. the porosity is heterogeneous, was developed by Correa and Borges [25] . Theses schemes are the basis of the REA Algorithm (“Reconstruction, Evolution and Average”), which is described below:

Figure 1. Schematic representation that decouples the problem with respect to time.

Reconstruction: From the computed values ?€??€? of the average ${S}_{i,j}^{n}$ and approximations of the derivatives ${\left({S}_{x}\right)}_{i,j}^{n}$ and ${\left({S}_{y}\right)}_{i,j}^{n}$ obtained using the MinMod limiters, we will rebuild the piecewise polynomial approximation that, for the two-dimensional case, is a bilinear interpolant given by:

${L}_{i,j}^{n}\left(x\right)={S}_{i,j}^{n}+{\left({S}_{x}\right)}_{i,j}^{n}\left(x-{x}_{i}\right)+{\left({S}_{y}\right)}_{i,j}^{n}\left(y-{y}_{j}\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{para}\text{\hspace{0.17em}}x\in {C}_{i,j}\text{.}$ (7)

Evolution: The hyperbolic equation at step ${t}^{n}$ is solved over time in order to obtain the average ${S}_{{D}_{i,j}}^{n+1}$ at time ${t}^{n+1}$ in a dual mesh. It is important to use information about the value of local velocity to define the size of the dual mesh, which will characterize the width of the Riemann fans (Leveque [26] ).

Average: We will project the average in a dual mesh ${S}_{{D}_{i,j}}^{n+1}$ of each cell on the original mesh ${\tau }_{h}$ in order to obtain the average in each cell in the original mesh at time ${t}^{n+1}$ .

3.2.1. Reconstruct

The linear approximations by parts are constructed from the average of each volume, based on Equation (7). We will approach the derivatives using the MinMod limiters:

${\left({S}_{x}\right)}_{i,j}^{n}=\frac{1}{\Delta x}\text{MinMod}\left[\alpha \left({S}_{i,j}^{n}-{S}_{i-1,j}^{n}\right),\frac{1}{2}\left({S}_{i+1,j}^{n}-{S}_{i-1,j}^{n}\right),\alpha \left({S}_{i+1,j}^{n}-{S}_{i,j}^{n}\right)\right]$ (8)

${\left({S}_{y}\right)}_{i,j}^{n}=\frac{1}{\Delta y}\text{MinMod}\left[\alpha \left({S}_{i,j}^{n}-{S}_{i,j-1}^{n}\right),\frac{1}{2}\left({S}_{i,j+1}^{n}-{S}_{i,j-1}^{n}\right),\alpha \left({S}_{i,j+1}^{n}-{S}_{i,j}^{n}\right)\right]$ (9)

where

$\text{MinMod}\left(a,b\right)=\frac{1}{2}\left[\text{sgn}\left(a\right)+\text{sgn}\left(b\right)\right]\text{min}\left(|a|,|b|\right)\text{.}$ (10)

The parameter a depends on the CFL conditions and varies from $\alpha =1$ , which corresponds to the MinMod basic limiter, and $\alpha =2$ , which is usually the value taken as upper limit (see Nessyahu and Tadmor [23] ; Kurganov and Tadmor [2] ).

3.2.2. Evolution

We construct the dual mesh based on the information about the local propagation velocity on the faces of the volume. We begin with the definition of discontinuous saturations on the sides of the central cell

${s}_{l}^{+}\left(y\right)={L}_{i,j}^{n}\left({x}_{i-1/2},y\right),\text{ }{s}_{l}^{-}\left(y\right)={L}_{i-1,j}^{n}\left({x}_{i-1/2},y\right)$ (11)

${s}_{r}^{+}\left(y\right)={L}_{i+1,j}^{n}\left({x}_{i+1/2},y\right),\text{ }{s}_{r}^{-}\left(y\right)={L}_{i,j}^{n}\left({x}_{i+1/2},y\right)$ (12)

${s}_{d}^{+}\left(x\right)={L}_{i,j}^{n}\left(x,{y}_{j-1/2}\right),\text{ }{s}_{d}^{-}\left(x\right)={L}_{i,j-1}^{n}\left(x,{y}_{j-1/2}\right)$ (13)

${s}_{u}^{+}\left(x\right)={L}_{i,j+1}^{n}\left(x,{y}_{j+1/2}\right),\text{ }{s}_{u}^{-}\left(x\right)={L}_{i,j}^{n}\left(x,{y}_{j+1/2}\right)$ (14)

where the indices $l,r,d,u$ refer to each of central cell faces and the signals (+) and (−) represent the saturations on the right and left, respectively (see Figure 2).

Figure 2. Second order linear reconstruction in the cells ${C}_{i-1,j}$ and ${C}_{i,j}$ (schematic drawing from Correa and Borges [25] ).

The maximum local wave propagation velocity on each side $\gamma$ ( $\gamma =l,r,u,d$ ) can be computed from the following relationship:

${a}_{\gamma }=\frac{\underset{s\in \left[{s}_{\gamma }^{\mathrm{min}},{s}_{\gamma }^{\mathrm{max}}\right]}{\text{max}}\left(\frac{\text{d}f}{\text{d}S}\right)}{\text{min}\left({\varphi }_{\gamma },{\varphi }_{c}\right)}$ (15)

where: ${s}_{\gamma }^{\text{min}}=\text{min}\left[{s}_{\gamma }^{-}\left(x\right),{s}_{\gamma }^{+}\left(x\right)\right]$ , ${s}_{\gamma }^{\text{max}}=\text{max}\left[{s}_{\gamma }^{-}\left(x\right),{s}_{\gamma }^{+}\left(x\right)\right]$ , ${\varphi }_{\gamma }$ is the porosity of the volumes adjacent to the central volume ${\varphi }_{c}$ (see Nessyahu and Tadmor [25] ; Kurganov and Tadmor [2] ).

Assuming that this velocity is constant for each time step, we can state that the information which comes from the left side, for example, travels the half cell in the x direction in a time given by

$\Delta {t}_{l}=\frac{\Delta x}{2{a}_{l}}.$

Extending this reasoning to the other sides, we have a CFL constraint for the time step given by

$\Delta {t}_{c}\le \mathrm{min}\left\{\Delta {t}_{l},\Delta {t}_{r},\Delta {t}_{d},\Delta {t}_{u}\right\}.$ (16)

Taking a time step $\Delta {t}_{c}$ from instant ${t}^{n}$ satisfying constraint (??), we can de_ne the following reference points:

${x}_{ln}={x}_{i-1/2}-{a}_{l}\Delta \stackrel{¯}{t},\text{ }{x}_{lc}={x}_{i-1/2}+{a}_{l}\Delta \stackrel{¯}{t}$

${x}_{rn}={x}_{i+1/2}+{a}_{r}\Delta \stackrel{¯}{t},\text{ }{x}_{rc}={x}_{i+1/2}-{a}_{r}\Delta \stackrel{¯}{t}$

${y}_{dn}={y}_{j-1/2}-{a}_{d}\Delta \stackrel{¯}{t},\text{ }{y}_{dc}={y}_{j-1/2}+{a}_{d}\Delta \stackrel{¯}{t}$

${x}_{un}={y}_{j+1/2}+{a}_{u}\Delta \stackrel{¯}{t},\text{ }{y}_{uc}={y}_{j+1/2}-{a}_{u}\Delta \stackrel{¯}{t}$

where the indices c and n refer to points related to the central cell and its neighboring cells respectively. These points are shown in Figure 3 and are used to build new cells called dual mesh which, in turn, is represented by the following ordered pairs:

Figure 3. Dual cells and reference points (schematic drawing from Correa e Borges [25] ).

${C}_{l}=\left({x}_{ln},{x}_{lc}\right)×\left({y}_{j-1/2},{y}_{j+1/2}\right);$

${C}_{r}=\left({x}_{rc},{x}_{rn}\right)×\left({y}_{j-1/2},{y}_{j+1/2}\right);$

${C}_{d}=\left({x}_{i-1/2},{x}_{i+1/2}\right)×\left({y}_{dn},{y}_{dc}\right);$

${C}_{u}=\left({x}_{i-1/2},{x}_{i+1/2}\right)×\left({y}_{uc},{y}_{un}\right);$

${C}_{c}=\left({x}_{lc},{x}_{rc}\right)×\left({y}_{dc},{y}_{uc}\right).$

The evolution step is performed by integrating the conservation law, Equation (6), in the dual cells ${C}_{\gamma }$ ( $\gamma =l,r,u,d,c$ ) from time ${t}^{n}$ to $\stackrel{¯}{t}={t}^{n}+\Delta {t}_{c}$ , i.e.,

${\int }_{{C}_{\gamma }}\text{ }\varphi w\left(x,\stackrel{¯}{t}\right)\text{d}\Omega ={\int }_{{C}_{\gamma }}\text{ }\varphi w\left(x,{t}^{n}\right)\text{d}\Omega -{\int }_{{t}^{n}}^{\stackrel{¯}{t}}{\int }_{{\Gamma }_{{C}_{\gamma }}}f\left(w\left(x,t\right)\right)\cdot n\text{d}{\Gamma }_{{C}_{\gamma }}\text{d}\stackrel{¯}{t}$ (17)

where $w\left(x,\stackrel{¯}{t}\right)$ represents the linear reconstruction, Equation (7), and ${\Gamma }_{{C}_{\gamma }}$ denotes the boundary of the cell. Then, without loss of generality, we present the development to cell ${C}_{l}$ .

Now we define the discontinuous saturations at the midpoints of the sides

${S}_{l}^{±}={s}_{l}^{±}\left({y}_{j}\right),\text{ }{S}_{r}^{±}={s}_{r}^{±}\left({y}_{j}\right),$

${S}_{d}^{±}={s}_{d}^{±}\left({x}_{i}\right),\text{ }{S}_{u}^{±}={s}_{u}^{±}\left({x}_{i}\right).$

The first integral on the right hand side of Equation (17) can be computed as follows:

${\int }_{{C}_{l}}\text{ }\varphi w\left(x,{t}^{n}\right)\text{d}\Omega ={\int }_{{C}_{ln}}\text{ }{\varphi }_{l}{L}_{i-1,j}\left(x\right)\text{d}\Omega +{\int }_{{C}_{lc}}\text{ }{\varphi }_{c}{L}_{i,j}\left(x\right)\text{d}\Omega$ (18)

where ${C}_{ln}={C}_{l}\cap {C}_{i-1,j}$ and ${C}_{lc}={C}_{l}\cap {C}_{i,j}$ . Using the linear reconstructions (7) we can show that, according to Correa and Borges [25] ,

$\begin{array}{l}{\int }_{{C}_{l}}\text{ }\varphi w\left(x,{t}^{n}\right)\text{d}\Omega \\ =\frac{|{C}_{l}|}{2}\left\{{\varphi }_{l}\left[{S}_{i-1,j}^{n}+\frac{\Delta x}{2}{\left({S}_{x}\right)}_{i-1,j}^{n}\right]+{\varphi }_{c}\left[{S}_{i,j}^{n}-\frac{\Delta x}{2}{\left({S}_{x}\right)}_{i,j}^{n}\right]\right\}+O{\left(\Delta \stackrel{¯}{t}\right)}^{2}\\ =\frac{|{C}_{l}|}{2}\left\{{\varphi }_{l}{S}_{l}^{-}+{\varphi }_{c}{S}_{l}^{+}\right\}+O{\left(\Delta \stackrel{¯}{t}\right)}^{2}\end{array}$ (19)

where the length of the dual cell ${C}_{l}$ is given by $|{C}_{l}|=2{a}_{l}\Delta \stackrel{¯}{t}\Delta y$ .

For the approximation of the flux term, represented by the second integral on the right hand side of Equation (17), we applied the following approximation at time ${t}^{n}$ and the values of saturation at midpoint of the face, i.e.

${\int }_{{t}^{n}}^{\stackrel{¯}{t}}{\int }_{{\Gamma }_{{C}_{l}}}f\left(w\left(x,t\right)\right)\cdot n\text{d}\Gamma \text{d}\stackrel{¯}{t}\approx \Delta \stackrel{¯}{t}\Delta y\left\{f\left({S}_{lc}^{n}\right)-f\left({S}_{ln}^{n}\right)\right\}+O{\left(\Delta \stackrel{¯}{t}\right)}^{2}$ (20)

where ${S}_{lc}^{n}={L}_{i,j}^{n}\left({x}_{lc},{y}_{j}\right)$ and ${S}_{ln}^{n}={L}_{i-1,j}^{n}\left({x}_{ln},{y}_{j}\right)$ . Taking the average of Equation (17) over the cell ${C}_{l}$ we get

$\begin{array}{l}\frac{1}{|{C}_{l}|{\varphi }_{{C}_{l}}}{\int }_{{C}_{l}}\text{ }\varphi \left(x\right)w\left(x,\stackrel{¯}{t}\right)\text{d}\Omega \\ \approx \frac{1}{2{\varphi }_{{C}_{l}}}\left\{{\varphi }_{l}{S}_{l}^{-}+{\varphi }_{c}{S}_{l}^{+}\right\}-\frac{1}{2{a}_{l}{\varphi }_{{C}_{l}}}\left\{f\left({S}_{lc}^{n}\right)-f\left({S}_{ln}^{n}\right)\right\}+O\left(\Delta \stackrel{¯}{t}\right)\end{array}$ (21)

where ${\varphi }_{{C}_{l}}=\frac{{\varphi }_{c}+{\varphi }_{l}}{2}$ . Thus, we can define the average saturation over cell ${C}_{l}$ ,

obtained after the evolution on a time step $\Delta {t}_{c}$ as:

${S}_{{C}_{l}}=\frac{1}{2{\varphi }_{{C}_{l}}}\left\{\left({\varphi }_{c}{S}_{l}^{+}+{\varphi }_{l}{S}_{l}^{-}\right)-\frac{1}{{a}_{l}}\left[f\left({S}_{lc}^{n}\right)-f\left({S}_{ln}^{n}\right)\right]\right\}+O\left(\Delta \stackrel{¯}{t}\right).$ (22)

With a similar development, we obtain for the remaining dual cells the expressions for the average saturation:

${S}_{{C}_{r}}=\frac{1}{2{\varphi }_{{C}_{r}}}\left\{\left({\varphi }_{c}{S}_{r}^{-}+{\varphi }_{r}{S}_{r}^{+}\right)-\frac{1}{{a}_{r}}\left[f\left({S}_{rn}^{n}\right)-f\left({S}_{rc}^{n}\right)\right]\right\}+O\left(\Delta \stackrel{¯}{t}\right)$ (23)

${S}_{{C}_{d}}=\frac{1}{2{\varphi }_{{C}_{d}}}\left\{\left({\varphi }_{c}{S}_{d}^{+}+{\varphi }_{d}{S}_{d}^{-}\right)-\frac{1}{{a}_{d}}\left[g\left({S}_{dc}^{n}\right)-g\left({S}_{dn}^{n}\right)\right]\right\}+O\left(\Delta \stackrel{¯}{t}\right)$ (24)

${S}_{{C}_{u}}=\frac{1}{2{\varphi }_{{C}_{u}}}\left\{\left({\varphi }_{c}{S}_{u}^{-}+{\varphi }_{u}{S}_{u}^{+}\right)-\frac{1}{{a}_{u}}\left[g\left({S}_{un}^{n}\right)-g\left({S}_{uc}^{n}\right)\right]\right\}+O\left(\Delta \stackrel{¯}{t}\right)$ (25)

$\begin{array}{c}{S}_{{C}_{c}}={S}_{i,j}^{n}+\Delta \stackrel{¯}{t}\frac{|{C}_{i,j}|}{|{C}_{c}|}\left\{\frac{1}{2}\left[\left({a}_{l}-{a}_{r}\right){\left({S}_{x}\right)}_{i,j}^{n}+\left({a}_{d}-{a}_{u}\right){\left({S}_{y}\right)}_{i,j}^{n}\right]\begin{array}{c}\stackrel{}{}\\ \end{array}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\frac{f\left({S}_{rc}^{n}\right)-f\left({S}_{lc}^{n}\right)}{{\varphi }_{c}\Delta x}-\frac{g\left({S}_{uc}^{n}\right)-g\left({S}_{dc}^{n}\right)}{{\varphi }_{c}\Delta y}\right\}+O{\left(\Delta \stackrel{¯}{t}\right)}^{2}\end{array}$ (26)

where $|{C}_{i,j}|=\Delta x\Delta y$ and $|{C}_{c}|=\left[\Delta x-\left({a}_{l}+{a}_{r}\right)\Delta \stackrel{¯}{t}\right]\left[\Delta y-\left({a}_{d}+{a}_{u}\right)\Delta \stackrel{¯}{t}\right]$ . For more details on the development of the evolution step of the REA algorithm (see Kurganov and Tadmor [2] ; Correa and Borges [25] ).

3.2.3. Projection and Semidiscrete Formulation

Once the average saturations in the cells of the dual mesh are computed, we proceed by projecting them over the volume ${C}_{i,j}$ of the original mesh at time $\stackrel{¯}{t}$ , i.e.,

${\stackrel{¯}{S}}_{i,j}^{n}=\frac{1}{2|{C}_{i,j}|}\left\{|{C}_{l}|{S}_{{C}_{l}}+|{C}_{r}|{S}_{{C}_{r}}+|{C}_{d}|{S}_{{C}_{d}}+|{C}_{u}|{S}_{{C}_{u}}+2|{C}_{c}|{S}_{{C}_{c}}\right\}-{I}_{{C}_{\beta }\cap {C}_{\theta }}$ (27)

where the term

${I}_{{C}_{\beta }\cap {C}_{\theta }}=\underset{\beta =l,r}{\sum }\underset{\theta =d,u}{\sum }|{C}_{\beta }\cap {C}_{\theta }|\frac{{S}_{{C}_{\beta }}+{S}_{{C}_{\theta }}}{2}$ (28)

corrects the overlap that occurs in the regions of intersection of the dual cells ${C}_{\beta }$ and ${C}_{\theta }$ . It can be seen that $|{x}_{\beta c}-{x}_{\beta }|$ and $|{y}_{\beta c}-{y}_{\beta }|$ are $O\left(\Delta \stackrel{¯}{t}\right)$ which implies ${I}_{{C}_{\beta }\cap {C}_{\theta }}=O{\left(\Delta \stackrel{¯}{t}\right)}^{2}$ . Knowing this, using the expressions for the average saturations in the dual mesh (Equations (22)-(26)) and remembering that $|{C}_{\beta }|$ and $|{C}_{\theta }|$ represent the volumes of the dual cell, we obtain as a result

${\stackrel{¯}{S}}_{i,j}^{n}=\frac{\Delta \stackrel{¯}{t}}{\Delta x}\left\{{a}_{l}{S}_{{C}_{l}}+{a}_{r}{S}_{{C}_{r}}\right\}+\frac{\Delta \stackrel{¯}{t}}{\Delta y}\left\{{a}_{d}{S}_{{C}_{d}}+{a}_{u}{S}_{{C}_{u}}\right\}+\frac{|{C}_{c}|}{|{C}_{i,j}|}{S}_{{C}_{c}}+O{\left(\Delta \stackrel{¯}{t}\right)}^{2}.$ (29)

The term of the projection on the central cell can be written from Equation (26) as

$\begin{array}{c}\frac{|{C}_{c}|}{|{C}_{i,j}|}{S}_{{C}_{c}}={S}_{i,j}^{n}\frac{|{C}_{c}|}{|{C}_{i,j}|}+\Delta \stackrel{¯}{t}\left\{\frac{1}{2}\left[\left({a}_{l}-{a}_{r}\right){\left({S}_{x}\right)}_{i,j}^{n}+\left({a}_{d}-{a}_{u}\right){\left({S}_{y}\right)}_{i,j}^{n}\right]\begin{array}{c}\stackrel{}{}\\ \end{array}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\frac{f\left({S}_{rc}^{n}\right)-f\left({S}_{lc}^{n}\right)}{{\varphi }_{c}\Delta x}-\frac{g\left({S}_{uc}^{n}\right)-g\left({S}_{dc}^{n}\right)}{{\varphi }_{c}\Delta y}\right\}+O{\left(\Delta \stackrel{¯}{t}\right)}^{2}.\end{array}$ (30)

$\frac{|{C}_{c}|}{|{C}_{i,j}|}{S}_{i,j}^{n}\left\{1-\left({a}_{l}+{a}_{r}\right)\frac{\Delta \stackrel{¯}{t}}{\Delta x}-\left({a}_{d}+{a}_{u}\right)\frac{\Delta \stackrel{¯}{t}}{\Delta y}\right\}+O{\left(\Delta \stackrel{¯}{t}\right)}^{2}$ (31)

enabling us to write, by replacing (31) in (30) and using the definition of saturations on the faces,

$\begin{array}{c}\frac{|{C}_{c}|}{|{C}_{i,j}|}{S}_{{C}_{c}}={S}_{i,j}^{n}-\frac{\Delta \stackrel{¯}{t}}{\Delta x}\left\{\left({a}_{l}{S}_{l}^{+}+{a}_{r}{S}_{r}^{-}\right)+\frac{1}{{\varphi }_{c}}\left[f\left({S}_{rc}^{n}\right)-f\left({S}_{lc}^{n}\right)\right]\right\}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\frac{\Delta \stackrel{¯}{t}}{\Delta y}\left\{\left({a}_{d}{S}_{d}^{+}+{a}_{u}{S}_{u}^{-}\right)+\frac{1}{{\varphi }_{c}}\left[g\left({S}_{uc}^{n}\right)-g\left({S}_{dc}^{n}\right)\right]\right\}+O{\left(\Delta \stackrel{¯}{t}\right)}^{2}\text{.}\end{array}$ (32)

Substituting into Equation (29) we have:

$\begin{array}{c}\frac{{\stackrel{¯}{S}}_{i,j}-{S}_{i,j}^{n}}{\Delta \stackrel{¯}{t}}=\frac{1}{\Delta x}\left\{{a}_{l}\left({S}_{{C}_{l}}-{S}_{l}^{+}\right)+{a}_{r}\left({S}_{{C}_{r}}-{S}_{r}^{-}\right)-\frac{1}{{\varphi }_{c}}\left[f\left({S}_{rc}^{n}\right)-f\left({S}_{lc}^{n}\right)\right]\right\}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{1}{\Delta y}\left\{{a}_{d}\left({S}_{{C}_{d}}-{S}_{d}^{+}\right)+{a}_{u}\left({S}_{{C}_{u}}-{S}_{u}^{-}\right)-\frac{1}{{\varphi }_{c}}\left[g\left({S}_{uc}^{n}\right)-g\left({S}_{dc}^{n}\right)\right]\right\}+O\left(\Delta \stackrel{¯}{t}\right).\end{array}$ (33)

This is the basic equation to obtain the semidiscrete formulation. To this end, we take the limit when $\Delta \stackrel{¯}{t}\to 0$ , which gives:

$\begin{array}{l}\underset{\Delta \stackrel{¯}{t}\to 0}{\mathrm{lim}}\frac{{\stackrel{¯}{S}}_{i,j}-{S}_{i,j}^{n}}{\Delta \stackrel{¯}{t}}\\ =\frac{1}{\Delta x}\underset{\Delta \stackrel{¯}{t}\to 0}{\mathrm{lim}}\left[{a}_{l}\left({S}_{{C}_{l}}-{S}_{l}^{+}\right)+{a}_{r}\left({S}_{{C}_{r}}-{S}_{r}^{-}\right)\right]-\frac{1}{{\varphi }_{c}\Delta x}\underset{\Delta \stackrel{¯}{t}\to 0}{\mathrm{lim}}\left[f\left({S}_{rc}^{n}\right)-f\left({S}_{lc}^{n}\right)\right]\\ -\frac{1}{\Delta y}\underset{\Delta \stackrel{¯}{t}\to 0}{\mathrm{lim}}\left[{a}_{d}\left({S}_{{C}_{d}}-{S}_{d}^{+}\right)+{a}_{u}\left({S}_{{C}_{u}}-{S}_{u}^{-}\right)\right]-\frac{1}{{\varphi }_{c}\Delta y}\underset{\Delta \stackrel{¯}{t}\to 0}{\mathrm{lim}}\left[g\left({S}_{uc}^{n}\right)-g\left({S}_{dc}^{n}\right)\right].\end{array}$ (34)

Noting that,

$\left\{\begin{array}{l}\underset{\Delta \stackrel{¯}{t}\to 0}{\mathrm{lim}}{S}_{lc}^{n}={S}_{l}^{+},\text{ }\underset{\Delta \stackrel{¯}{t}\to 0}{\mathrm{lim}}{S}_{\mathrm{ln}}^{n}={S}_{l}^{-}\\ \underset{\Delta \stackrel{¯}{t}\to 0}{\mathrm{lim}}{S}_{rc}^{n}={S}_{r}^{-},\text{ }\underset{\Delta \stackrel{¯}{t}\to 0}{\mathrm{lim}}{S}_{rn}^{n}={S}_{r}^{+}\\ \underset{\Delta \stackrel{¯}{t}\to 0}{\mathrm{lim}}{S}_{dc}^{n}={S}_{d}^{+},\text{ }\underset{\Delta \stackrel{¯}{t}\to 0}{\mathrm{lim}}{S}_{dn}^{n}={S}_{d}^{-}\\ \underset{\Delta \stackrel{¯}{t}\to 0}{\mathrm{lim}}{S}_{uc}^{n}={S}_{u}^{-},\text{ }\underset{\Delta \stackrel{¯}{t}\to 0}{\mathrm{lim}}{S}_{un}^{n}={S}_{u}^{+}\end{array}$

$\underset{\Delta \stackrel{¯}{t}\to 0}{\mathrm{lim}}\frac{{\stackrel{¯}{S}}_{i,j}-{S}_{i,j}^{n}}{\Delta \stackrel{¯}{t}}=\frac{\text{d}{S}_{i,j}}{\text{d}t}\text{.}$ (35)

Now substituting Equations (22)-(25) into expression (35) and making the necessary simplifications we arrive at the semidiscrete formulation:

$\begin{array}{c}\frac{\text{d}{S}_{i,j}}{\text{d}t}=\frac{1}{\Delta x}\left\{\frac{{a}_{l}{\varphi }_{l}}{{\varphi }_{c}+{\varphi }_{l}}\left({S}_{l}^{-}-{S}_{l}^{+}\right)+\frac{{a}_{r}{\varphi }_{r}}{{\varphi }_{c}+{\varphi }_{r}}\left({S}_{r}^{+}-{S}_{r}^{-}\right)-\frac{1}{{\varphi }_{c}+{\varphi }_{l}}\left[f\left({S}_{l}^{+}\right)-f\left({S}_{l}^{-}\right)\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\frac{1}{{\varphi }_{c}+{\varphi }_{r}}\left[f\left({S}_{r}^{+}\right)-f\left({S}_{r}^{-}\right)\right]-\frac{1}{{\varphi }_{c}}\left[f\left({S}_{r}^{-}\right)-f\left({S}_{l}^{+}\right)\right]\right\}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{1}{\Delta y}\left\{\frac{{a}_{d}{\varphi }_{d}}{{\varphi }_{c}+{\varphi }_{d}}\left({S}_{d}^{-}-{S}_{d}^{+}\right)+\frac{{a}_{u}{\varphi }_{u}}{{\varphi }_{c}+{\varphi }_{u}}\left({S}_{u}^{+}-{S}_{u}^{-}\right)-\frac{1}{{\varphi }_{c}+{\varphi }_{d}}\left[g\left({S}_{d}^{+}\right)-g\left({S}_{d}^{-}\right)\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\frac{1}{{\varphi }_{c}+{\varphi }_{u}}\left[g\left({S}_{u}^{+}\right)-g\left({S}_{u}^{-}\right)\right]-\frac{1}{{\varphi }_{c}}\left[g\left({S}_{u}^{-}\right)-g\left({S}_{d}^{+}\right)\right]\right\}.\end{array}$ (36)

According to Kurganov and Tadmor [2] , this formulation can be written in terms of numerical fluxes by means of the expression:

$\frac{\text{d}{S}_{i,j}}{\text{d}t}=\frac{{H}_{r}-{H}_{l}}{\Delta x}+\frac{{H}_{u}-{H}_{d}}{\Delta y}$ (37)

where the numerical fluxes are given by

${H}_{\beta }=\frac{{a}_{\beta }{\varphi }_{\beta }}{{\varphi }_{c}+{\varphi }_{\beta }}\left({S}_{\beta }^{+}-{S}_{\beta }^{-}\right)-\frac{1}{{\varphi }_{c}}\left[\frac{{\varphi }_{c}\varphi \left({S}_{\beta }^{+}\right)+{\varphi }_{\beta }\varphi \left({S}_{\beta }^{-}\right)}{{\varphi }_{c}+{\varphi }_{\beta }}\right]\text{,}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{para}\text{\hspace{0.17em}}\beta =r,u$ (38)

${H}_{\beta }=\frac{{a}_{\beta }{\varphi }_{\beta }}{{\varphi }_{c}+{\varphi }_{\beta }}\left({S}_{\beta }^{+}-{S}_{\beta }^{-}\right)-\frac{1}{{\varphi }_{c}}\left[\frac{{\varphi }_{\beta }\varphi \left({S}_{\beta }^{+}\right)+{\varphi }_{c}\varphi \left({S}_{\beta }^{-}\right)}{{\varphi }_{c}+{\varphi }_{\beta }}\right]\text{,}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{para}\text{\hspace{0.17em}}\beta =l,d$ (39)

with $\phi =f$ for $\beta =r,l$ and $\phi =g$ for $\beta =u,d$ , (see Equation (6)). We use the fourth order Runge-Kutta method to solve the system of ordinary differential Equations (37).

3.3. Corrector Problem

We describe the formulation of the corrector step used for both the saturation and concentration equations. For brevity, the procedure will be presented only for the saturation.

Knowing the solution obtained in the predictor step ${S}^{*\left(n+1\right)}$ and the field $\varphi \left(x,t\right)$ , find ${S}_{c}^{**}\left(x,{t}_{n+1}\right)$ for each time subinterval $\Delta {t}_{c}=\left({t}_{n},{t}_{n+1}\right)$ , $n=0,\cdots ,N-1$ with ${t}_{0}=0$ , ${t}_{n}=T$ , such that:

$\varphi \frac{\partial {S}_{c}^{**}}{\partial t}=-{S}_{c}^{**}\frac{\partial \varphi }{\partial t}$ (40)

with the initial condition

${S}_{c}^{**}\left(x,{t}_{n}\right)={S}_{c}^{*\left(n+1\right)}\text{.}$ (41)

Integrating Equation (40) over the interval $\Delta {t}_{c}$ for each point $x\in \Omega$ , we obtain:

${S}^{**n+1}=\frac{{\varphi }^{n}{S}^{*n+1}}{{\varphi }^{n+1}}.$ (42)

3.4. Dissolution

According to Obi and Blunt [1] , based in [27] , in order to compute the number of moles of CO2 in the liquid phase ${T}_{d}^{c}$ and update the values of saturation ${S}_{c}$ and concentration C at time ${t}^{n+1}$ , we can use the algorithm

${\left({T}_{d}^{c}\right)}^{n+1}={\varphi }^{n+1}\left[\left(1-{S}_{c}^{n+1}\right){C}^{n+1}+{m}_{c}{S}_{c}^{n+1}\right]$ (43)

where

$\text{if}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\left({T}_{d}^{c}\right)}^{n+1}\ge {\varphi }^{n+1}{K}_{d}\to \left(\begin{array}{l}{C}^{n+1}={K}_{d}\hfill \\ {S}_{c}^{n+1}=\frac{{\left({T}_{d}^{c}\right)}^{n+1}/{\varphi }^{n+1}-{K}_{d}}{{m}_{c}-{K}_{d}}\hfill \end{array}$

$\text{if}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\left({T}_{d}^{c}\right)}^{n+1}<{\varphi }^{n+1}{K}_{d}\to \left(\begin{array}{l}{C}^{n+1}={\left({T}_{d}^{c}\right)}^{n+1}/{\varphi }^{n+1}\hfill \\ {S}_{c}^{n+1}=0\hfill \end{array}$

The dissolution constant is defined as

${K}_{d}=\frac{{x}_{{\text{CO}}_{\text{2}}}^{L}{m}_{w}}{1-{x}_{{\text{CO}}_{\text{2}}}^{L}\left(1-\frac{{m}_{w}}{{m}_{c}}\right)}.$ (44)

where ${m}_{w}$ is the molar density of water without CO2 dissolved, ${m}_{c}$ is the molar density of CO2 without water dissolved and ${x}_{{\text{CO}}_{\text{2}}}^{L}$ is the mole fraction of CO2 which can dissolve in the aqueous phase.

3.5. Reaction

In order to solve Equation (3), we compute the source term of reaction ( ${T}_{r}$ ) and update the porosity field f. Moreover, an empirical law will be used to update the permeability field K.

In agreement with [1] we compute the total number of moles of CO2 per unit volume in both phases ( ${m}_{tt}$ ) using the following algorithm

${m}_{tt}={\varphi }^{n}\left[\left(1-{S}_{c}^{n}\right){C}^{n}+{m}_{c}{S}^{n}\right].$ (45)

The total number of moles of CO2 per unit volume which effectively reacts with the rock is defined as

${T}_{r}=\frac{{m}_{w}}{1-\frac{{m}_{w}}{{m}_{c}}}.$ (46)

To calculate the variation of porosity, we use the following procedure

$\left(\begin{array}{ll}{\varphi }^{n+1}={\varphi }^{n}\left[\mathrm{exp}\left(-K\Delta {t}_{r}\right)\right]\hfill & \text{if}\text{\hspace{0.17em}}{m}_{tt}\ge {T}_{r}{\varphi }^{n}\left[1-\mathrm{exp}\left(-K\Delta {t}_{r}\right)\right]\hfill \\ {\varphi }^{n+1}={\varphi }^{n}-\frac{{m}_{tt}}{{T}_{r}}\hfill & \text{if}\text{\hspace{0.17em}}{m}_{tt}<{T}_{r}{\varphi }^{n}\left[1-\mathrm{exp}\left(-K\Delta {t}_{r}\right)\right]\hfill \end{array}$

Finally, the permeability is updated from the empirical formulation (Wellman et al. [14] ):

${K}^{n+1}={K}^{n}{\left(\frac{{\varphi }^{n+1}}{{\varphi }^{n}}\right)}^{3.4}.$ (47)

4. Results and Discussion

For a given velocity field ${v}_{Dt}\left(0,t\right)=0.18\text{\hspace{0.17em}}\text{m}/\text{year}$ , an injection of CO2 was performed for a period of 20 years ( ${S}_{c}\left(0,t\right)=0.8$ ) and the saturation profile was recorded over a period of 180 years after the injection ceased. The following physical parameters were considered: temperature of 353 K; water density of 1050 kg×m−3; CO2 density of 710 kg×m−3; water viscosity of 5 ´ 10−4 Pa×s; CO2 viscosity of 6 ´ 10−5 Pa×s; CO2 molar density ( ${m}_{c}$ ) of 16140 moles×m−3; effective molar density of the rock ( ${T}_{r}$ ) of 26.5 moles×m−3; porosity of $\varphi =0.15$ ; absolute permeability of $K=100\text{\hspace{0.17em}}\text{mD}$ ; water residual saturation of ${S}_{rw}=0.2$ ; CO2 residual saturation of ${S}_{rc}=0.2$ . We use a reservoir of dimensions $2000\text{\hspace{0.17em}}\text{m}×1500\text{\hspace{0.17em}}\text{m}×200\text{\hspace{0.17em}}\text{m}$ and a mesh of 400 elements. We adopt the point of injection of 750 m from the origin and pressure difference between this point and the right edge is assumed to be 29 MPa. During 20 years, the CO2 was injected in a homogeneous porous media and the saturation profiles were recorded after the injection stopped for three cases: 1) disregarding dissolution and reaction, the flow is due to advection and dispersion only, non reactive rock, such as quartz, is considered (see Figure 4), 2) allowing dissolution but no reaction (see Figure 5), and 3) both effects considered (see Figure 6).

The Buckley-Leverett profile (see [29] ) coincides with the CO2 saturation

Figure 4. Saturation profile during 20 years of injection, and 180 years after injection has stopped desconsidering dissolution and precipitation.

Figure 5. Saturation profile during 20 years of injection, and 180 years after injection has stopped considering dissolution.

Figure 6. Saturation profile during 20 years of injection, and 180 years after injection has stopped considering dissolution and precipitation.

recorded during injection, which shows the quality of the algorithm compared to the analytical solution (see Figure 7). Moreover, after the end of the injection3, the saturation profile was recorded for 180 years and the same physical trends presented in [1] were observed. First, comparing the effect of the dissolution with the inert case, in considering dissolution, the saturation profile length was slightly shorter. The reason for it is due to the fact that the CO2 was transferred to the aqueous phase (see Figure 4 and Figure 5). Lastly, when the dissolution and reaction were considered simultaneously, an expressive impact on the saturation profile could be noticed4. Such effect was due to the fact that the conditions of this simulation provided sufficient time for a meaningful amount of precipitation (see Figure 6).

5. Conclusions

320 years.

4During 180 years after the injection has stopped.

In this paper, we presented a locally conservative numerical method for the simulation of geological storage of CO2. The innovative aspects of the new methodology are:

• The time scale of each physical phenomenon was properly treated using the operator splitting technique;

• An extension of the KT method was develop for a system of equations with source terms and applied to solve numerically the carbon sequestration model adopted;

• The advantages of the KT method were inherited by the extended version, i.e., this method avoids the difficulties that arise when adopting small time steps enforced by CFL stability restrictions.

In addition, the impact of the dissolution and precipitation of CO2 during the

Figure 7. Comparison in between the analytical solution of Buckley-Leverett (see [29] ) and numerical solution obtained by the proposed algorithm.

injection in a homogeneous porous medium was analyzed. The results obtained are in agreement with the analytical solution for the injection phase and with the saturation profiles available in the literature for the phase after the injection ceased.

Acknowledgements

The author would like to acknowledge CAPES (Coordination for the Improvement of Higher Education Personnel) for their financial support; the professors José A. R. Parise and Ivan F. M. Menezes for their valuable contributions to this work; and Marisa R. Shanstrom and Fernando Saint-Martin Soares for their editorial assistance.

Conflicts of Interest

The authors declare no conflicts of interest.

Cite this paper

Sica, L. (2018) Numerical Study of the Injection of Carbon Dioxide in a Homogeneous Porous Media. Open Journal of Fluid Dynamics, 8, 115-132. doi: 10.4236/ojfd.2018.81009.

 [1] Obi, E.-O.I. and Blunt, M.J. (2006) Streamline-Based Simulation of Carbon Dioxide Storage in a North Sea Aquifer. Water Resources Research, 42, 1-13. https://doi.org/10.1029/2004WR003347 [2] Kurganov, A. and Tadmor, E. (2000) New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection-Diffusion Equations. Journal of Computational Physics, 160, 241-282. https://doi.org/10.1006/jcph.2000.6459 [3] Houghton, J.T., Ding, Y., Griggs, D.J., Noguer, M., Van Der Linden, P.J. and Xiasu, D. (2001) Climate Change 2001: The Scientific Basis Contribution of Working Group I. Technical Report, Intergovernmental Panel on Climate Change (IPCC). Cambridge University Press, Cambridge. [4] Ravagnani, A. and Suslick S. (2008) II. Modelo dinamico de sequestro geol′ogico de CO2 em reservat′orios de petr′oleo. Revista Brasileira de Ciencias, 38, 39-60. [5] Meer, L.V. (1993) The Conditions Limiting CO2 Storage in Aquifers. Energy Conversion Management, 34, 959-966. https://doi.org/10.1016/0196-8904(93)90042-9 [6] Meer, L.V. (1995) The CO2 Storage Efficiency in Aquifers. Energy Conversion Management, 36, 513-518. https://doi.org/10.1016/0196-8904(95)00056-J [7] Meer, L.V. (1996) Computer Modelling of Underground CO2 Storage. Energy Conversion Management, 37, 1155-1160. https://doi.org/10.1016/0196-8904(95)00313-4 [8] Holt, T., Jensen, J, and Lindeburg, E. (1995) Underground storage of CO2 in Aquifers and Oils Reservoir. Energy Conversion Management, 36, 535-538. https://doi.org/10.1016/0196-8904(95)00061-H [9] Weir, G., White, S., and Kissling, W. (1995) Reservoir Storage and Containment of Greenhouse Gases. Energy Conversion Management, 36, 531-534. https://doi.org/10.1016/0196-8904(95)00060-Q [10] Law, D. and Bachu, S. (1996) Hydrogeological and Numerical Analysis of CO2 Disposal in Deep Sedimentary Aquifers in the Alberta Sedimentary Basin. Energy Conversion and Management, 37, 1167-1174. https://doi.org/10.1016/0196-8904(95)00315-0 [11] Lindberg, E. (1997) Escape of CO2 from Aquifers. Energy Convers. Manage, 38, 235-240. https://doi.org/10.1016/S0196-8904(96)00275-0 [12] Johnson, J., Steefel, J., Nitao, J. and Knauss, K. (2000) Reactive Transport Modeling of Subsurface CO2 Sequestration: Identification of Optimal Target Reservoir and Evaluation of Performance Based on Geochemical, Hydrologic and Structural Constraints. 8th International Forum of the International Energy Foundation, Las Vegas, July 2000, 23-28. [13] Ennis-King, J. and Paterson, L. (2002) Engineering Aspects of Geological Se-questration of Carbon Dioxide. In SPE Asia Pacific Oil and Gas Conference and Exhibition, Melbourne, 8-10 October 2002, 8-10. https://doi.org/10.2118/77809-MS [14] Wellman, T., Grigg, R., McPherson, B., Svec, R. and Lichtner, P. (2003) Evaluation of CO2-Brine-Reservoir Rock Interaction with Laboratory flow Test and Reactive Transport Modeling. In International Symposium on Oilfield Chemistry, Houston, 5-7 February, 5-7. https://doi.org/10.2118/80228-MS [15] Pruess, K., Xu, T., Apps, J. and Garcia, J. (2003) Numerical Modeling of Aquifer Disposal of CO2. SPE Journal, SPE-83695-PA, 49-60. https://doi.org/10.2118/83695-PA [16] Xu, T., Apps, J.A. and Pruess, K. (2003) Reactive Geochemical Transport Simulation to Study Mineral Trapping for CO2 Disposal in Deep Saline Arenaceous Aquifers. Journal of Geophysical Research, 108, 1-12. [17] Kumar, A., Ozah, M., Pope, G., Bryant, S., Sepehmoori, K. and Lake, L. (2005) Reservoir Simulation of CO2 Storage in Deep Saline Aquifers. SPE Journal, 10, 336-348. https://doi.org/10.2118/89343-PA [18] Malik, Q.M. and Islam, M.R. (2000) CO2 Injection in the Weyburn Field of Canada: Optimization of Enhanced Oil Recovery and Greenhouse Gas Storage with Horizontal Wells. SPE/DOE Improved Oil Recovery Symposium, Tulsa, 3-5 April 2000, SPE-59327-MS. [19] Fanchi, J. (2001) Geological Sequestration: Modeling and Monitoring Injected CO2. SPE/EPA/DOE Exploration and Production Environmental Conference, San Antonio, 26-28 February 2001, SPE-66749-MS. [20] Spiteri, E., Juanes, R., Blunt, M. and Orr Jr., F. (2005) Relative Permeability Hysteresis: Trapping Models and Application to Geological CO2 Sequestration. SPE Annual Technical Conference and Exhibition, Dallas, 9-12 October 2005, SPE-96448-MS. [21] Calabrese, M., Masserano, F. and Blunt, M. (2005) Simulation of Physical-Chemical Processes during Carbon Dioxide Sequestration in Geological Structures. SPE Annual Technical Conference and Exhibition, Dallas, 9-12 October 2005, SPE-95820-MS. https://doi.org/10.2118/95820-MS [22] Ennis-King, J. and Paterson, L. (2005) Role of Convective Mixing in the Long-Term Storage of Carbon Dioxide in Deep Saline Formations. SPE Journal, 10, 349-356. https://doi.org/10.2118/84344-PA [23] Nessyahu, H. and Tadmor, E. (1990) Non-Oscillatory Central Differencing for Hyperbolic Conservation Laws. Journal of Computational Physics, 87, 408-463. https://doi.org/10.1016/0021-9991(90)90260-8 [24] Kurganov, A. and Lin, C.T. (2007) On the Reduction of Numerical Dissipation in Central-Upwind Schemes. Communications in Computational Physics, 2, 141-163. [25] Correa, M.R. and Borges, M.R. (2013) A Semi Discrete Central Scheme for Scalar Hyperbolic Conservation Laws with Heterogeneous Storage Coefficient and Its Application to Porous Media Flow. International Journal for Numerical Methods in Fluids, 73, 205-224. https://doi.org/10.1002/fld.3794 [26] Leveque, R.J. (2002) Finite Volume Methods for Hyperbolic Problems. Cambridge Texts in Applied Mathematics, Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511791253 [27] Smith, J.M., Van Ness, H. and Abbott, M. (2004) Introduction to Chemical Engineering Thermodynamics. McGraw-Hill Education, New York. [28] Chen, Z., Huan, G. and Ma, Y. (2006) Computational Methods for Multiphase Flows in Porous Media. Computational Science and Engineering, Society for Industrial and Applied Mathematics. [29] Bear, J. (1972) Dynamics of Fluids in Porous Media. American Elsevier Publishing Company, Amsterdam.