Long-time evolution of sequestered CO$_2$ in porous media

CO$_2$ sequestration in subsurface reservoirs is important for limiting atmospheric CO$_2$ concentrations. However, a complete physical picture able to predict the structure developing within the porous medium is lacking. We investigate theoretically reactive transport in the long-time evolution of carbon in the brine-rock environment. As CO$_2$ is injected into a brine-rock environment, a carbonate-rich region is created amid brine. Within the carbonate-rich region minerals dissolve and migrate from regions of high concentration to low concentration, along with other dissolved carbonate species. This causes mineral precipitation at the interface between the two regions. We argue that precipitation in a small layer reduces diffusivity, and eventually causes mechanical trapping of the CO$_2$. Consequently, only a small fraction of the CO$_2$ is converted to solid mineral; the remainder either dissolves in water or is trapped in its original form. We also study the case of a pure CO$_2$ bubble surrounded by brine and suggest a mechanism that may lead to a carbonate-encrusted bubble due to structural diffusion.


I. INTRODUCTION
The sequestration of CO 2 in geological formations is widely considered to be an important approach for mitigating the rise of atmospheric CO 2 levels [1][2][3][4][5]. Deep saline aquifers and gas fields are primarily chosen for storage [1,4,6]. Supercritical CO 2 is injected into these porous media while displacing another fluid, brine [5,7]. The CO 2 then gradually dissolves; however, its long-term fate remains poorly understood [8]. Here we address the long time scale and study theoretically the evolution of the CO 2 after injection into the brine-rock system.
As the CO 2 is injected into the brine-rock environment, it initially becomes trapped, either by a physical mechanism in the presence of low permeability rocks, or by retention as a separate phase in the pore space due to interfacial tension [9]. The disordered structure of the void spaces forces the injected fluid along certain paths that create regions, or bubbles, of the injected fluid amid regions of the defending fluid, and vice versa [10][11][12][13]. This process is known as invasion percolation when, as for the supercritical CO 2 , the invading fluid is non-wetting. For two immiscible fluids, the fluid configuration is often determined by the structure of the rock and by surface tension effects. Further, the two-phase system can become unstable when reactant particles migrate from one phase to the other and change the chemical composition of each phase [14,15]. Within the high CO 2 phase, minerals dissolve; diffusion causes minerals and carbonate species to migrate from high concentration to low concentration regions, and a fraction of them precipitates. In nature, this process can be seen in hot springs, when a bubble of oxygen emerges from photosynthetic cyanobacteria [16,17]. The high gradient of CO 2 between the bubble and the surrounding water leads to loss of CO 2 in the vicinity of the bubble, drives up the saturation level, and eventually a crust is created around the bubble. In general, mineral precipitation on a small boundary layer at the interface may lead to lower diffusivity and slower kinetics. In the carbon sequestration process, this may cause a mechanical trapping of the CO 2 bubble and lower the solidification rate of the carbon minerals.
Here we develop theoretical understanding of this process of mechanical phase separation. We consider two scales: At the microscale, a single CO 2 bubble is surrounded by brine in the void space of a porous medium. In this case, the reactions occur at the interface between the bubble and the brine, as the CO 2 dissolves and reduces the pH in its vicinity. The macroscale averages over many such bubbles. In this case, a high concentration of the invaded CO 2 changes the properties of a macroscopic region. The region becomes more acidic and no precipitation occurs unless the carbonate species migrate to a different region.
The paper begins by addressing the macroscale problem and the mathematical background of the reactive diffusion equation. We then study the mobility change in the fluid-rock system due to the precipitate minerals. Finally, we consider the microscale case of a single CO 2 bubble amid brine and suggest a second mechanism that may also lead to separation and self-sealing.

II. MACROSCOPIC REACTIVE TRANSPORT
The injection of supercritical CO 2 (scCO 2 ) into a porous medium, initially occupied with another fluid, i.e., brine (salty water), generates two different regions in the system [14,18]. The first region is where the CO 2 displaced  [14]. Microscale: CT image of material distribution in Frio sandstone, adapted from [21]. (Upper right) Region 1 consists of CO 2 bubbles (black) surrounded by brine (green) and grains of the rock (red). (Lower right) Region 2 is the brine-rock system.
most of the existing brine. In this region, the void space in the porous medium is filled with bubbles of CO 2 and water saturated with carbonate species. The CO 2 dissolves rapidly in the water to reach equilibrium which results in low pH, and also a high concentration of dissolved minerals [19,20]. The second region is the intact brine-rock system, which is characterized by a low concentration of CO 2 and higher pH. Fig. 1 depicts the complexity of the CO 2 -brine-rock system. The existence of two phases and concentration gradients drives the components to migrate from one phase to another. As they diffuse, they react to reach a local equilibrium. Although pure scCO 2 may be clogged due to mechanical or capillary trapping [1,22], in the presence of water it can dissolve into its ionic forms, i.e. bicarbonate and carbonic acid, until it reaches a local thermodynamical equilibrium [23]. These carbonate species may diffuse more easily through the brine. Within the high CO 2 phase, minerals dissolve because of the acidic environment. They then migrate from high concentration to low concentration regions and a fraction of them precipitates. The evolution of each component can be described by the reactive diffusion equation Here, the C i , i = 1, .., 5, are the concentration of the CO 2 , HCO − 3 , CO 2− 3 , H + and the mineral Ca 2+ , respectively. D i is the isotropic diffusion coefficient for the ith component, and R i is the reaction rate defined by the carbonate system [23,24] and the dissolution and precipitation of calcite mineral. The latter can be expressed as [25][26][27], where M represents the density of precipitated calcite. The rate coefficient is defined by [26] where A(r, t) is the reactive surface area, and k +i are the rate constants for the forward reactions. The saturation ratio [25] is defined by the ion activity product divided by the solubility constant of calcite. Details of the reactions and the values of the reaction constants can be found in ref. [26,28]. Whether the calcite will precipitate or dissolve is determined by the value of Ω; when Ω < 1 dissolution occurs; when Ω > 1 mineral precipitates, and when Ω = 0 the system is at equilibrium.

III. NUMERICAL SIMULATION
In our model, we assume that there is no diffusion of the solid mineral as it nucleates and precipitates at the surface of the rock. Also, the capillary forces between the scCO 2 and the water prevent the mixing between the two phases; therefore the diffusion constant of CO 2 is also set to zero. Since the dissolution of calcite is considerably fast compared to dissolution of other minerals, we assume that, after injection, calcite dissolves quickly until it reaches equilibrium. Only new calcite that has been precipitated can be dissolved again. Thus, the reactive surface area is set to zero, as long as there is no previous precipitation of calcite.
We solve the reactive diffusion equations (Eqs. (1) and (2)) for each component using Galerkin finite elements method on quadratic triangular grid with a 4 th order Runge-Kutta integration scheme. The simulation starts when the carbonate system, the pH, and the dissolved calcium mineral are at local equilibrium in each phase. In region 1, the carbonate-rich brine-rock region (see Fig. 1), we set the pH to 6 and the total dissolved inorganic carbon to 10 −2 mol/liter. In region 2, we set pH= 7 and the carbonate concentration at 10 −3 mol/liter. The simulation starts with the concentration of each one of the carbonate species and the dissolved calcium mineral at local equilibrium in each region.
We consider a carbonate-rich circular shape (region 1) is surrounded by brine (region 2). We find that the accumulation of calcite, shown in Fig. 2a, occurs on a small boundary layer close to the interface. The precipitation profile is approximated by the typical diffusion length l ∼ √ Dt, shown in Fig. 2b.   3 shows qualitatively how the interface curvature alters the mineral precipitation. We find that the accumulated crust is highly dependent on the curvature of the interface due to the nature of a diffusion process through a curved interface [29]; negative curvature with respect to the inner bubble of the CO 2 phase has maximum mineral precipitation at the interface.
Collectively, these results indicate that precipitation of mineral occurs on a small boundary layer at the interface and a carbonate crust, created in a small boundary layer, can lead to a mechanical separation of the two phases.

IV. SYSTEM SIZE DEPENDENCE
Mineral precipitation along the interface changes the porosity of the rock, decreases the effective diffusivity, and may clog existing voids [1,22]. The decrease in the mobility of the ions in the solution could bring further solidification of minerals to a halt. Whether a pathway between pores will be clogged or not is highly dependent on the local density of the accumulated mineral, which depends also on the amount of the carbonate species, and the size and the shape of the domain noted as region 1. To study how a change in the size alters the maximum mineral density, we initiate a 1-dimensional domain of size d of region 1, and run the simulation until equilibrium. During the simulation, the mineral density reaches a maximum and then decreases as the solution at the interface become more acidic. We (1). Changing the equation into dimensionless variables, x → ξd, t → d 2 D τ , leads to and from Eq. 2, the density of the precipitated mineral becomes where A(τ ) and B(τ ) correspond to the concentration of Ca 2+ and CO 2− 3 . If A(τ ) and B(τ ) are purely diffusion controlled, i.e. if the reaction term in Eq. (5) can be neglected, the accumulated precipitation scales like d 2 . In practice, the exponent may be lower than 2 because of the contribution of the reaction terms that change the equilibrium state in our system.

V. EFFECTIVE DIFFUSIVITY
The diffusivity in a porous medium depends on several factors which are related to the pore geometry and the effective porosity accessible by diffusion [30,31]. For simplicity, we consider the permeability reduction to be linear with the porosity loss and mineral precipitation, and the effective diffusivity D e to be linear with the porosity φ: where D 0 is the bulk diffusivity corresponding to the initial porosity, φ t is the porosity available for transport, θ = M is the maximum density of the precipitated mineral and θ c is a critical accumulated density in which D e vanishes.  Fig. 4 with D e = D 0 . Red squares: The results allowing the porosity to decrease and the effective diffusivity D e to decrease as the precipitation grows. Note the enhanced precipitation when the critical porosity, θ c , is reached.
At this point, particles cannot cross between the regions, and the system seeks a local chemical equilibrium.
As the accumulated mineral reaches the critical porosity, it changes locally the permeability and mechanical separation may occur due to diffusivity loss and self-sealing. We also observe that while the diffusion coefficient decreases, it creates even more precipitation in a smaller area close to the interface. These results, shown in Fig. 5, are in agreement with the diffusion length scale of l ∼ √ D e t from the numerical approximation and the mineral precipitation of M ∼ 1/D e from Eq. (6). In addition, when the system is clogged, i.e. when the effective diffusivity goes to zero, the system seeks a local chemical equilibrium.
From the numerical results, we can identify two trapping mechanisms that occur as particles migrate from one region to the other. For small carbonate-rich regions, the CO 2 dissolves completely into the brine and the low pH does not allow a significant precipitation of the carbonate minerals. Thus, the CO 2 is trapped in its dissolved forms. For larger regions, solidification of minerals stops as clogging occurs, and the developed crust separates the two regions. In both cases the total amount of the precipitated carbonate minerals is small compare to concentrations of the other carbonate species.

VI. THE MICROSCALE: A SINGLE CO2 BUBBLE
In this section, we discuss a mechanism that leads to a carbonate-encrusted bubble, in which a single bubble of pure CO 2 surrounded by brine develops a crust at the interface. Unlike the upscaled case discussed above, where each region consists of the same ingredients but with a different concentration, here we have a pure bubble (no upscaling) with only CO 2 , either in a liquid phase or gas phase. The problem is pictured in Fig. 6. In this case, no chemical reactions occur inside the CO 2 bubble; reactions instead occur only at the interface with the brine. The CO 2 then dissolves into the brine, creates charged carbonate species and reduces the pH in its vicinity. Most of the components diffuse slowly away from the interface, however, the protons in water diffuse faster due to structural diffusion, also known as the Grotthoss mechanism [32]. The concentration gradient with unequal diffusivities generates an electrical field that slows down the rapidly diffusing protons, but also attracts positively charged minerals toward the interface. (The dynamics is described by the Nernst-Planck equations [33,34]). This leads to supersaturation at the interface and mineral precipitation. This process suggest that an isolated bubble can remain stable due to a self-sealing mechanism.

VII. DISCUSSION AND CONCLUSION
We have shown that in a system of reactive fluids, a gradient in the concentration between regions leads to a supersaturation at the interface and to precipitation and porosity loss. A local change in the porosity reduces permeability and may cause mechanical separation. This process is important in understanding the long-term sequestration of carbon dioxide in subsurface geological formation. We predict that a small domain of a region consisting mostly of the carbonate species in an acidic environment will dissolve completely into the brine. Larger domains are more likely to be self-sealed, and only a fraction of the carbonate species will be precipitated. For a single CO 2 bubble in brine solution, we suggest a mechanism that causes self-sealing due to the pH gradient at the interface of the bubble and structural diffusion.
Our results suggest that only a small fraction of the injected CO 2 is converted to a solid mineral. The remainder stays in its dissolved ionic form or is trapped in its original form. Whether a domain will go into dissolution trapping or mechanical separation depends on the concentration gradients, the properties of the rock and the porosity available for transport.