Modeling of Carbon Dioxide Leakage from Storage Aquifers

Long-term geological storage of CO2 in deep saline aquifers offers the possibility of sustaining access to fossil fuels while reducing emissions. However, prior to implementation, associated risks of CO2 leakage need to be carefully addressed to ensure safety of storage. CO2 storage takes place by several trapping mechanisms that are active on different time scales. The injected CO2 may be trapped under an impermeable rock due to structural trapping. Over time, the contribution of capillary, solubility, and mineral trapping mechanisms come into play. Leaky faults and fractures provide pathways for CO2 to migrate upward toward shallower depths and reduce the effectiveness of storage. Therefore, understanding the transport processes and the impact of various forces such as viscous, capillary and gravity is necessary. In this study, a mechanistic model is developed to investigate the influence of the driving forces on CO2 migration through a water saturated leakage pathway. The developed numerical model is used to determine leakage characteristics for different rock formations from a potential CO2 storage site in central Alberta, Canada. The model allows for preliminary analysis of CO2 leakage and finds applications in screening and site selection for geological storage of CO2 in deep saline aquifers.


Introduction
It is generally accepted that increase in concentrations of greenhouse gases (GHG) such as carbon dioxide (CO 2 ), methane (CH 4 ) and nitrous oxide (N 2 O) is the main cause of global warming [1].It has been reported that the anthropogenic emissions have increased atmospheric concentrations of CO 2 by around 35% during the last 200 years, from 280 ppm in 1800 to 380 ppm in 2005 [1].The rate of anthropogenic emissions to the atmosphere is also continuing to increase.In 1950, anthropogenic emissions of CO 2 were around 6 Gt of CO 2 ; in 2000, this number was around 22 Gt CO 2 /year; and in 2015 the number was estimated to be approximately 37 Gt CO 2 /year [2].From this amount, approximately 45% of the emitted anthropogenic carbon remains in the atmosphere, and the other 55% goes into the ocean and land biomass [2].Therefore, to mitigate the effects on climate change, CO 2 emissions to the atmosphere should be reduced, and the concentration of CO 2 needs to be stabilized at around 450 ppm by 2050 [1].One attractive option for reducing CO 2 emissions is carbon capture and storage (CCS) [1].CCS involves capturing CO 2 from energy related sources such as large industrial plants, transporting it to a suitable storage location and storing it in deep geological formations, before it is emitted to the atmosphere.
Among the proposed CCS options, sequestration in deep saline aquifers appears to be the most promising option due to the relatively larger storage capacity and worldwide distribution of these formations [1].Deep saline aquifers are permeable layers (usually found in sedimentary basins), saturated with salt water and bounded from below and above by less permeable rocks.Proper analysis Fluids 2018, 3, 80 3 of 31 dimensionless gravitational and capillary numbers and examined their effects on storage capacity of reservoirs.They showed that a low ratio of gravitational to viscous forces (low gravitational number) and to a lesser extent a high ratio of capillary to viscous forces (high capillary number) leads to a higher CO 2 storage capacity.
Possible leakage pathways could be in the form of natural interruptions and breaches through the confining strata such as open faults, fractures and erosional channels, or they could be in the form of manmade or artificial pathways such as abandoned wells or activated faults and fractures [13,14].
Figure 1 shows pathways for CO 2 migration for a possible leakage.Inactive faults can be reactivated because of local pressure changes during CO 2 injection.Likewise, closed fractures may open by exceeding bottom hole pressure from the minimum in-situ stress.There is an active debate on whether injection of CO 2 will trigger large earthquakes and reactivate faults and therefore threaten the seal integrity of CO 2 storage reservoirs [15,16].The main geomechanical aspects affecting the integrity of the CO 2 storage sites have been reviewed by different researchers [17,18].For large scale injection operations, domains in the order of thousands of square kilometers need to be analyzed to guarantee the safety of storage [19].
Fluids 2018, 3, x FOR PEER REVIEW 3 of 32 and to a lesser extent a high ratio of capillary to viscous forces (high capillary number) leads to a higher CO2 storage capacity.
Possible leakage pathways could be in the form of natural interruptions and breaches through the confining strata such as open faults, fractures and erosional channels, or they could be in the form of manmade or artificial pathways such as abandoned wells or activated faults and fractures [13,14].
Figure 1 shows pathways for CO2 migration for a possible leakage.Inactive faults can be reactivated because of local pressure changes during CO2 injection.Likewise, closed fractures may open by exceeding bottom hole pressure from the minimum in-situ stress.There is an active debate on whether injection of CO2 will trigger large earthquakes and reactivate faults and therefore threaten the seal integrity of CO2 storage reservoirs [15,16].The main geomechanical aspects affecting the integrity of the CO2 storage sites have been reviewed by different researchers [17,18].For large scale injection operations, domains in the order of thousands of square kilometers need to be analyzed to guarantee the safety of storage [19].Leakage of CO2 can be assessed by means of numerical and analytical models.Analytical solutions can only be derived with sufficient restrictive assumptions.Various analytical solutions have been published in the literature, especially for the leakage along existing wells [20][21][22].Additionally, theoretical models are developed for investigating diffusive leakage of brine from an aquifer into overburden and underburden formations during geological storage of CO2 [23,24].Although several efforts have been devoted to assess the risk of leakage through existing wells, leakage through faults and fractures has been less addressed in the literature, and it is partly due to the difficulty in characterizing the features of such pathways.Useful numerical simulations for CO2 leakage through existing/activated faults have been published in the literature [25][26][27][28].Recently, Kang et al. [29] used a simplified model to derive expressions for pressure drawdown and interface upconing around a leaky fault.They described a vertical two-phase flow along a fault and concluded that, for estimation of leakage, structures and properties of the fault must be included in the model.
Even with large potential capacity of saline aquifers, a major obstacle to implement CO2 sequestration in such formations is the lack of insights into the risks associated with it.In relatively short time scales, CO2 injected into saline aquifer is in free phase and prone to leakage.Therefore, for Leakage of CO 2 can be assessed by means of numerical and analytical models.Analytical solutions can only be derived with sufficient restrictive assumptions.Various analytical solutions have been published in the literature, especially for the leakage along existing wells [20][21][22].Additionally, theoretical models are developed for investigating diffusive leakage of brine from an aquifer into overburden and underburden formations during geological storage of CO 2 [23,24].Although several efforts have been devoted to assess the risk of leakage through existing wells, leakage through faults and fractures has been less addressed in the literature, and it is partly due to the difficulty in characterizing the features of such pathways.Useful numerical simulations for CO 2 leakage through existing/activated faults have been published in the literature [25][26][27][28].Recently, Kang et al. [29] used a simplified model to derive expressions for pressure drawdown and interface upconing around a leaky fault.They described a vertical two-phase flow along a fault and concluded that, for estimation of leakage, structures and properties of the fault must be included in the model.
Even with large potential capacity of saline aquifers, a major obstacle to implement CO 2 sequestration in such formations is the lack of insights into the risks associated with it.In relatively short time scales, CO 2 injected into saline aquifer is in free phase and prone to leakage.Therefore, for selection of storage sites, it is important to know, in case of a leaky fault or fracture, how much CO 2 will leak and how long it will take to happen.The combined effect of different mechanisms and their evolution with time and space can only be evaluated through sophisticated numerical simulations.The objective of this paper is to develop a simple model which allows for quick and relatively easy screening of storage reservoirs by investigating processes and driving forces that are active on CO 2 plume in relatively short time scales after CO 2 injection.

Model Description and Governing Equations
To study migration of CO 2 along a water-saturated fault zone, a one-dimensional (1D) model is developed.To provide conservative estimates of the leakage potential, we ignore the effect of dissolution trapping and the system is assumed to be isothermal.
CO 2 (non-wetting phase) and formation brine (wetting phase) are assumed to be immiscible and the interphase transfer between phases (mass transfer between the fluids, i.e., dissolution of CO 2 in brine and evaporation of brine into CO 2 ) is neglected.
Mass conservation equation and the extended Darcy formula are the governing equations of two immiscible fluids flow in porous layer for a linear system which can be written separately for each phase.The modified Darcy's law for multi-phase flow, which relates the volumetric flux of each phase to the respective phase pressure, can be written as [5]: where → v α is the Darcy velocity [LT −1 ] of phase α, k is the absolute (intrinsic) permeability [L 2 ] of the formation, k rα is the relative permeability of phase α, µ α is the viscosity [ML −1 T −1 ] of phase α, p α is the pressure [ML −1 T −2 ] of phase α, ρ α is the density [ML −3 ] of phase α, g is the gravitational [LT −2 ] acceleration, z is the vertical coordinate [L] and positive upward, and λ α = k rα /µ α is the mobility of phase α.
Mass balance equation for each phase is given by [5]: where φ is formation porosity and q α represents source/sink term [T −1 ].By substitution of Equation (1) into Equation (2), a system of partial differential equations is obtained for each phase.Additional closure equations are required to solve the system of the partial differential equations.First, the sum of phase saturations equals unity.
Second, the pressure difference between the two phases equals capillary pressure p c [5]: Finally, to solve the system of partial differential equations, constitutive relations are required, which include capillary pressure and relative permeability as functions of saturation.In this study, the following assumptions have been made to simplify the modeling of two-phase flow.The porous layer (leakage pathway) is considered homogenous, isotropic and one-dimensional tilted with an angle θ from the horizontal axis.Formation brine (w) and gas/supercritical phase (nw) are single-component fluids and are assumed to be incompressible.The flow is considered isothermal, sources and sinks are not present and the dynamic viscosity is also assumed to be constant.These assumptions allow development of a simple and fast model for analysis of leakage and have been previously used in many studies.Nevertheless, modeling of CO 2 leakage from leakage pathways in the absence of these assumptions requires detailed geological and numerical models.Simulation of CO 2 leakage including all details, such as 3D features, pathway heterogeneities, phase behavior and phase change, geomechanical effects, thermal, and geochemical effects, is very time consuming and difficult, if not impossible.In addition, simple models often provide a first order estimate of CO 2 leakage necessary for site screening and risk assessment.A schematic of the problem is shown in Figure 2. phase change, geomechanical effects, thermal, and geochemical effects, is very time consuming and difficult, if not impossible.In addition, simple models often provide a first order estimate of CO2 leakage necessary for site screening and risk assessment.A schematic of the problem is shown in Figure 2. Using the above-mentioned assumptions and adding the definition of capillary pressure (Equation ( 4)), CO2 flux can be written as: With substitution in mass balance equation, saturation equation for CO2 phase can be written as: where v v v Equation ( 6) is a second-order nonlinear partial differential equation, which is of the general form for describing flow of two immiscible incompressible fluids in a one-dimensional porous media.To develop the dimensionless form of Equation ( 6), the following dimensionless variables are defined.Using the above-mentioned assumptions and adding the definition of capillary pressure (Equation ( 4)), CO 2 flux can be written as: With substitution in mass balance equation, saturation equation for CO 2 phase can be written as: where v t = v nw + v w is constant total velocity and is defined as the summation of CO 2 and brine velocities.In addition, λ t = λ nw + λ w and ∆ρ = ρ w − ρ nw are total mobility and density difference, respectively.Equation ( 6) is a second-order nonlinear partial differential equation, which is of the general form for describing flow of two immiscible incompressible fluids in a one-dimensional porous media.To develop the dimensionless form of Equation ( 6), the following dimensionless variables are defined.
Fluids 2018, 3, 80 6 of 31 In the above equations, S is the normalized gas saturation varying between 0 and 1, S wr is the irreducible (or residual) saturation of brine, and S cc is the critical gas saturation.Fluid saturations are often normalized with respect to the ranges of values occurring in the problem under consideration. 1 − S wr and S cc are the maximum and minimum saturations that occur during the migration of CO 2 , respectively.Dimensionless variables z D and t D are dimensionless length and time (0 ≤ z D ≤ 1, t D ≥ 0), respectively.Parameters l cr and t cr are the characteristic values for length and time, respectively.Choosing the correct characteristic scale is of great importance and requires a well-developed understanding of the system being analyzed since it can change the dimensionless numbers (see Section 3.1) by orders of magnitudes [11].In this study, the total length of the system is used to scale the migration length (l cr = L) and t cr is the time at which the migrating CO 2 travels the length of the system by pure advection.
With substitutions, the dimensionless velocity of CO 2 (v D ) is obtained as: and the gas saturation equation (Equation ( 6)) therefore can be written as: To determine solution of Equation ( 12), proper saturation functions for capillary pressure and relative permeability are required.In this study, the following relative permeability relations as function of saturation are used in the form of Brooks-Corey equation [30].
In the above equation, k max rnw is the relative permeability of CO 2 at irreducible brine saturation, and k max rw is the relative permeability to brine at the critical CO 2 saturation.For the drainage cycle, k max rw = 1 and S cc = 0. Parameters nc and nb are the constants (also known as Corey's coefficients), which affect the shape of the relative permeability curves.Increasing values of these coefficients cause the relative permeability curve to become concave.The capillary pressure is represented by a logarithmic function, which has been conventionally used for such problems [31].
where p d is the capillary entry pressure.Different dimensionless numbers can be established to define the balance between driving forces taking place during and after CO 2 injection.For the CO 2 injection application, Kopp et al. [11] used dimensionless capillary and gravitational numbers and qualitatively examined their effects on storage capacity.Here, we use a similar form of dimensionless numbers for our system.
Substitution of Equations ( 13)- (15) into Equations (11) and (12) gives: ) where v DD , v DG and v DV are capillary, gravity and viscous contribution to dimensionless CO 2 velocity, respectively.N Ca is the dimensionless capillary number and N Gr is the dimensionless gravity number, which relate forces acting on the system and are analyzed based on the reservoir parameters as well as fluid properties as follow.D(S), G(S) and V(S) are capillary, gravity and viscous functions, respectively, and are defined as: where Equation ( 17) is used to simulate CO 2 migration in a linear porous media and considers viscous, capillary and gravity forces simultaneously.Equation ( 17) is a nonlinear second order partial differential equation, which is also known as nonlinear convection-diffusion equation.A similar form of this equation has been developed before for the application of water flooding in oil reservoirs [32].The first term in Equation ( 17) is the accumulation term; the second term includes the contribution of capillary forces in the displacement process and has a diffusive nature; the third term contains the gravitational contributions to the flow and has an advective character; and the fourth term is the contribution of viscous forces and same as the third term has an advective character.
The problem described by Equation ( 17) is completed by setting proper initial and boundary conditions.In this study, we assume that the entire domain is initially fully saturated with brine, therefore only primary drainage process is considered here, and the residual CO 2 is not a concern.At the inlet boundary, it is considered that CO 2 saturation is always at its maximum value.This assumption seems reasonable when CO 2 has already displaced the formation brine beneath a leakage pathway and a residual saturation of the wetting phase has been established.At the outlet boundary it is assumed that the changes in CO 2 saturation is negligible.This boundary condition infers negligible diffusive (capillary) flux at the outer boundary compared to gravity and viscous fluxes.The initial and boundary conditions are then given by: The model described above forms the basis of this study.Numerical discretization and detailed validation of the developed model is presented in Appendices A and B, respectively.

Results and Discussion
The developed model in the previous section is utilized to understand processes and driving forces active during leakage of CO 2 into the shallower formations.This section is organized as follows.First, the effect of dimensionless numbers and interaction of forces on the displacement process and leakage is investigated.Next, the developed model is used to assess leakage in seven formations located in Alberta, Canada, which are considered potential options for geological storage of CO 2 .

Interaction of Gravity, Capillary, and Viscous Forces
To investigate the effect of forces acting on the fluid flow and rate of leakage, we define seven different datasets.Detailed comparison and properties of each set is presented in Table 1.The following properties are assumed to be the same for all the seven datasets.Length of the domain (L) is considered to be 500 m.Density difference between CO 2 and brine (∆ρ = ρ w − ρ nw ) is assumed to be 300 kg/m 3 and porosity (φ) is 0.2.Brine and CO 2 viscosity are assumed to be 0.86 and 0.06 mPa•s, respectively.Absolute permeability is taken to be 100 mD and residual water saturation (S wr ) is assumed to be 0.3.Parameters used for these datasets are only for comparison and they are not related to a specific reservoir.The last three rows in Table 1 show the calculated mobility ratio, N Ca and N Gr .Dataset 1 is considered as base case.In each dataset, the changing parameter is bolded in Table 1.Two analysis groups have been established.First, the effect of changes in N Ca and N Gr is investigated, while D(S), G(S) and V(S) are kept constant (Datasets 2-4).Second, changes of relative permeability parameters of CO 2 and brine system are examined, which include the effect of mobility ratio and Corey exponents nb and nc (Datasets 5-7).Changes in gas saturation profile and gas dimensionless velocity are examined for each dataset.Breakthrough time is considered to be the time when migrating CO 2 appears at the outlet boundary.
Figure 3 shows the parameters map of N Gr − N Ca space.A reference equilibrium state between gravitational and viscous forces is indicated with dotted line at N Gr = 1 and between capillary and viscous forces at N Ca = 1.In Figure 3, the lower left sector represents domination of viscous forces over capillary and gravitational forces.As we move from this sector to the right or top, viscous forces lose their influence in favor of capillary (lower right sector) or gravitational (upper left sector) forces, respectively.Finally, in the upper right sector, capillary and gravitational forces dominate over viscous forces.
follows.First, the effect of dimensionless numbers and interaction of forces on the displacement process and leakage is investigated.Next, the developed model is used to assess leakage in seven formations located in Alberta, Canada, which are considered potential options for geological storage of CO2.

Interaction of Gravity, Capillary, and Viscous Forces
To investigate the effect of forces acting on the fluid flow and rate of leakage, we define seven different datasets.Detailed comparison and properties of each set is presented in Table 1.The following properties are assumed to be the same for all the seven datasets.Length of the domain ( L ) is considered to be 500 m.Density difference between CO2 and brine ( w nw
Absolute permeability is taken to be 100 mD and residual water saturation ( wr S ) is assumed to be 0.3.Parameters used for these datasets are only for comparison and they are not related to a specific reservoir.The last three rows in Table 1 show the calculated mobility ratio, Ca N and Gr N .Dataset 1 is considered as base case.In each dataset, the changing parameter is bolded in Table 1.Two analysis groups have been established.First, the effect of changes in Ca N and Gr N is investigated, while () DS , () GS and () VS are kept constant (Datasets 2-4).Second, changes of relative permeability parameters of CO2 and brine system are examined, which include the effect of mobility ratio and Corey exponents nb and nc (Datasets 5-7).Changes in gas saturation profile and gas dimensionless velocity are examined for each dataset.Breakthrough time is considered to be the time when migrating CO2 appears at the outlet boundary.
Figure 3 shows the parameters map of A reference equilibrium state between gravitational and viscous forces is indicated with dotted line at 1 Gr N = and between capillary and viscous forces at 1 Ca N = .In Figure 3, the lower left sector represents domination of viscous forces over capillary and gravitational forces.As we move from this sector to the right or top, viscous forces lose their influence in favor of capillary (lower right sector) or gravitational (upper left sector) forces, respectively.Finally, in the upper right sector, capillary and gravitational forces dominate over viscous forces.Capillary, gravity and viscous functions are plotted for all seven datasets in Figure 4a-c.The shape and value of these functions are influenced by Corey exponents in relative permeability equation (nc and nb), and mobility ratio.Capillary function D(S) (Figure 4a) and gravity function G(S) (Figure 4b) have a bell shape, whereas shape of V(S) is the same as the fractional flow equation in Buckley-Leveret problem.
In the first four datasets, dimensionless functions D(S), G(S) and V(S) are same.In Datasets 5-7, relative permeability parameters are changed to evaluate their effect on flow process.Migration of CO 2 plume depends on balance of these forces, which determines the security of storage.Capillary, gravity and viscous functions are plotted for all seven datasets in Figure 4a-c.The shape and value of these functions are influenced by Corey exponents in relative permeability equation (nc and nb), and mobility ratio.Capillary function () DS (Figure 4a) and gravity function () GS (Figure 4b) have a bell shape, whereas shape of () VS is the same as the fractional flow equation in Buckley-Leveret problem.
In the first four datasets, dimensionless functions () DS , () GS and () VS are same.In Datasets 5-7, relative permeability parameters are changed to evaluate their effect on flow process.Migration of CO2 plume depends on balance of these forces, which determines the security of storage.for Datasets 1 to 7 from Table 1.  1.
Figure 5 shows the normalized gas saturation and contribution of different forces in CO 2 dimensionless velocity versus dimensionless length at time t D = 0.2 for Datasets 1 to 4.
Fluids 2018, 3, x FOR PEER REVIEW 11 of 32 depicted in Figure 5b, increasing dimensionless velocity ( v D ) to values higher than one is more pronounced in this case due to the simultaneous effect of gravity and capillary.Early leakage of brine may be used as an indicator of subsequent leakage of CO2.However, in cases such as Dataset 4, significant back flow of brine may occur, which can lead to less brine leakage prior to CO2 leakage.This effect could be emphasized especially for higher values of

Ca N and
Gr N , where it could result in CO2 leakage to shallower depths without significant leakage of brine.For Dataset 1 (Table 1), N Ca and N Gr have the same magnitude and less than one, and position of this dataset in the parameter space map is in the lower left sector of Figure 3.This is translated to dominance of viscous forces over capillary and gravity forces.As shown in Figure 5, higher viscous forces lead to a CO 2 plume with sharp front.As a result, plume migration is slowest among other datasets (Datasets 2-4).In this case, gravitational forces are rather weak, leading to a low buoyancy-induced upward migration of CO 2 .
Fluids 2018, 3,80 In Dataset 2, the effect of capillary entry pressure, p d , is investigated.By changing the capillary entry pressure (p d ) from 0.2 MPa to 1.5 MPa, capillary number (N Ca ) increases up to 1.161, while the gravity number (N Gr ) remains the same as the one used in Dataset 1.This leads to a shift to the lower right corner of the parameter map, as shown in Figure 3. Increasing capillary number will increase the influence of capillary forces over viscous forces.Stronger capillary forces lead to a more diffusive-like front propagation.An earlier breakthrough is a result of higher capillary forces, as shown in the saturation distributions and velocity profiles for Datasets 1 and 2 in Figure 5a,b (and later in Figure 6a).As shown in Figure 5c, as N Ca increases in Dataset 2, contributions of capillary forces on CO 2 plume velocity (v DD ) increases.As the effect of capillary forces increases, the inlet dimensionless velocity (v D ) of the migrating CO 2 takes values greater than one.This is because the dimensionless velocity of the migrating CO 2 is given by v D = 1 − v b /v t .At high capillary forces, there is backward flow of brine as CO 2 migrates upward.This leads to v b values to be less than zero that results in v D values greater than one.In other words, a dimensionless velocity greater than one indicates backward flow of formation brine from the leakage pathway.This effect is particularly important in the early times and, as the plume evolves, the effect of capillary diminishes.We further study the effect of relative permeability parameters on the displacement process.Relative permeability depends on several factors such as pore size characteristics, wettability of fluids and phase saturation [33].Neglecting these dependencies will affect the interpretation of short and long-term fate of the injected CO2 in deep saline aquifers.Therefore, detailed laboratory measurements are necessary to predict the dependency of the relative permeability to saturation of phases present in the medium.In our model, shape of the relative permeability curve is determined by the Corey exponents and endpoint relative permeability to each phase.1) to show the effect of varying N Ca and N Gr ; and (b) Datasets 1 and 5-7 (Table 1) to show the effect of relative permeability parameters M, nc and nb.
The effect of gravitational forces is examined using Dataset 3 while N Ca is kept similar to Dataset 1.As the density difference or the tilt angle (θ) increase the effect of gravity on CO 2 plume migration increases.The gravity number (N Gr ) varies up to 1.139 for a tilt angle of θ = π/2 (or a vertical leakage pathway).By increasing N Gr , the effect of gravity forces on dimensionless velocity (v DG ) increases and therefore a more buoyant flow regime is established leading to a faster plume evolution for Dataset 3, compared to Dataset 1.Such a flow regime may be preferred since it creates an extended contact between the injected CO 2 and brine leading to higher dissolution of CO 2 in brine.However, increasing N Gr enhances the effect of gravity segregation and leads to an earlier breakthrough of CO 2 and possibly an enhanced risk of leakage, which is not desirable.
Next, we study the effect of total Darcy velocity by changing v t = 3 × 10 −6 m/s to a lower value of v t = 3 × 10 −7 m/s.Decreasing v t causes an increase in both N Ca and N Gr .However, since velocity appears in both numbers, their ratio will stay the same.This will lead to a move to the upper right corner of the parameter map shown in Figure 3.By moving to this sector, viscous forces lose their influence in favor of both capillary and gravity forces.The combined effect of increasing capillary and gravity forces will cause fastest evolution of CO 2 plume compared to other datasets.In Figure 6a, the average gas saturation is plotted against dimensionless time for Datasets 1-4.In terms of leakage, Dataset 1 results in a more compact and less diffusive front with delayed breakthrough of the leaked CO 2 compared to the other cases.However, in such a case, a sudden release of CO 2 is expected.
On the other hand, when capillary forces are dominant (Dataset 2), the leaking CO 2 forms a diffusive front and consequently an earlier breakthrough of the gradually leaking CO 2 is expected.This may be practically important since it provides more time to take remedial actions.Similarly, increasing gravity effect in Dataset 3 results in an earlier breakthrough time.It is worth mentioning that the effect of increasing N Gr by one order of magnitude on accelerating the breakthrough time, is much less than that of increasing N Ca .Finally, for Dataset 4, the breakthrough time is the lowest among the other datasets due to domination of both capillary and gravity forces.In addition, as depicted in Figure 5b, increasing dimensionless velocity (v D ) to values higher than one is more pronounced in this case due to the simultaneous effect of gravity and capillary.Early leakage of brine may be used as an indicator of subsequent leakage of CO 2 .However, in cases such as Dataset 4, significant back flow of brine may occur, which can lead to less brine leakage prior to CO 2 leakage.This effect could be emphasized especially for higher values of N Ca and N Gr , where it could result in CO 2 leakage to shallower depths without significant leakage of brine.
We further study the effect of relative permeability parameters on the displacement process.Relative permeability depends on several factors such as pore size characteristics, wettability of fluids and phase saturation [33].Neglecting these dependencies will affect the interpretation of short and long-term fate of the injected CO 2 in deep saline aquifers.Therefore, detailed laboratory measurements are necessary to predict the dependency of the relative permeability to saturation of phases present in the medium.In our model, shape of the relative permeability curve is determined by the Corey exponents and endpoint relative permeability to each phase.
In this section, effects of mobility ratio (which could be due to changes on k max rnw or viscosity ratio), nb and nc on plume migration is investigated through Datasets 5-7.These parameters affect shapes of D(S), V(S) and G(S) functions.For simplicity, and to illustrate the effects, we considered nc = nb = 2 for Datasets 1-5.The normalized gas saturation and contribution of different forces on CO 2 dimensionless velocity, versus dimensionless length is plotted in Figure 7 at time t D = 0.2 for Datasets 1 and 5-7.
for these cases.Therefore, in these cases, brine will be displaced along with CO2 to the shallower formation.
Figure 6b shows the average normalized saturation versus dimensionless time for Datasets 1 and 5-7.The results showed that a lower endpoint relative permeability of CO2 and therefore lower mobility ratio postpone the breakthrough time.Additionally, increasing nc increases the breakthrough time.It is important to note that the overall effect of the parameters that appear in both capillary and gravity numbers (

Ca N and
Gr N ) as well as capillary, gravity and viscous functions ( () DS , () GS , and () VS ), determine the controlling mechanisms.These parameters (density, viscosity, capillary pressure, relative permeability, etc.) are dependent on the in-situ conditions of temperature, pressure, salinity and the interaction of fluids (CO2 and brine) and rock.Therefore, for each case, dimensionless numbers and functions should be evaluated to determine the behavior of the system.To better show  In Dataset 5, k max rnw is taken to be 0.065 (one order less than k max rnw in Dataset 1).For a constant total velocity, reducing CO 2 endpoint relative permeability by one order of magnitude causes a decrease in M, N Ca and N Gr by the same order.As shown in Figure 3, this will cause a shift to the lower left sector.Reducing N Ca and N Gr will make the displacement process highly viscous dominant.Additionally, unlike previous datasets, a change in mobility ratio will cause a change in capillary, gravity and viscous functions.As shown in Figure 4a,b, a decrease in mobility ratio from M = 10 to M = 1 will dramatically increase the values of D(S) and G(S) functions.However, since the process is highly viscous dominant (due to small values of N Ca and N Gr ), the contribution of these forces on plume velocity is negligible.Influence of reducing M on V(S) is shown in Figure 4c (comparing Dataset 5 with Dataset 1).From the frontal advance theory, a higher average saturation of CO 2 behind front is obtained for Dataset 5.As expected, the results shown in Figure 7 reveal that a relative permeability curve with low k max rnw results in less CO 2 propagation velocity and a significant delay of the breakthrough time.
In Datasets 6 and 7, the effect of Corey exponents (as a measure of wettability) on behavior of CO 2 plume is investigated.Changes in Corey exponents will not change the dimensionless numbers of N Ca and N Gr and therefore position of various forces in the parameter space shown in Figure 3 will be the same as Dataset 1.However, the Corey exponents affects shapes of functions D(S), V(S) and G(S).In Dataset 6, nc is kept the same as the one in Dataset 5 (nc = 2), and nb is increased to 4. As nb increases, the relative permeability to brine, D(S) and G(S) decrease.As shown in Figure 4c, by increasing nb, the average saturation behind the front reduces.Therefore, increasing nb will result in a faster plume evolution, and a decrease in time of breakthrough.This effect is better shown by comparing Datasets 6 and 1 in Figure 7. Increasing nc, from 2 to 4 in Dataset 7 results in a reduction in relative permeability of CO 2 .Similar to Dataset 6, values of D(S) and G(S) will decrease, but the average saturation behind the front (Figure 4c) increases and therefore a more efficient displacement and a delayed breakthrough time is achieved in this scenario as compared to Dataset 6.
Since the process is viscous dominant in Datasets 1 and 5-7, contribution of gravity and capillary forces on CO 2 plume velocity is negligible, and no countercurrent displacement of brine is observed for these cases.Therefore, in these cases, brine will be displaced along with CO 2 to the shallower formation.
Figure 6b shows the average normalized saturation versus dimensionless time for Datasets 1 and 5-7.The results showed that a lower endpoint relative permeability of CO 2 and therefore lower mobility ratio postpone the breakthrough time.Additionally, increasing nc increases the breakthrough time.
It is important to note that the overall effect of the parameters that appear in both capillary and gravity numbers (N Ca and N Gr ) as well as capillary, gravity and viscous functions (D(S), G(S), and V(S)), determine the controlling mechanisms.These parameters (density, viscosity, capillary pressure, relative permeability, etc.) are dependent on the in-situ conditions of temperature, pressure, salinity and the interaction of fluids (CO 2 and brine) and rock.Therefore, for each case, dimensionless numbers and functions should be evaluated to determine the behavior of the system.To better show the effect of different forces, breakthrough time is calculated for different values of total Darcy velocity for Datasets 1 and 5 where mobility ratios are M = 10 and M = 1, respectively.By increasing the total Darcy velocity, N Ca and N Gr will decrease, reflecting the effect of increasing viscous forces.For instance, for M = 1, by increasing v t from 1 × 10 −9 m/s to 1 × 10 −3 m/s, values of N Ca and N Gr will both change from as high as 46.5 to 0.00005.Therefore, for lower values of v t , capillary and gravity forces dominate over viscous forces.Figure 8a shows effects of gravity forces.Continuous lines display the calculated breakthrough time when all forces are included, and the dashed lines show breakthrough time when gravity effects are neglected.The results demonstrate that gravity effects lead to an earlier breakthrough time for small values of v t due to increasing buoyancy effects.However, as v t increases (viscous dominated regime), the effect of gravity becomes negligible.Next, we study the effect of capillary forces on breakthrough time.As shown in Figure 8b, capillary forces have a significant effect on breakthrough time for small values of v t .In other words, neglecting capillary effects results in overestimation of breakthrough time.However, the influence of capillary forces becomes negligible by increasing v t .
Fluids 2018, 3, 80 becomes negligible.Next, we study the effect of capillary forces on breakthrough time.As shown in Figure 8b, capillary forces have a significant effect on breakthrough time for small values of v t .
In other words, neglecting capillary effects results in overestimation of breakthrough time.However, the influence of capillary forces becomes negligible by increasing v t .

Case Study on CO2 Storage Aquifers in Western Canada
In this section, the developed model is used to study potential CO2 storage sites in central Alberta, Canada.Province of Alberta has the largest greenhouse gas emissions in Canada, with annual emissions close to 274.1 Mt CO2 in 2015 [34].The province is underlain by the Alberta basin, which is suitable for CO2 geological sequestration in all parts except for its shallow northeastern corner [35].The largest concentration of CO2 sources in Alberta is in the Edmonton region where four coal-fired power plants are located near Wabamun Lake, west of Edmonton, which have combined annual CO2 emissions of more than 30 Mt CO2 [36]. Figure 9 shows the location of major CO2 emission sources in central Alberta, Canada.CO2 sequestration in deep saline aquifers in proximity of these power plants is a promising option for reducing CO2 atmospheric emissions since the deep coal seams and oil and gas reservoirs in local area do not have sufficient capacity for sequestration of CO2.Because CO2 sequestration in geological media and especially in deep saline aquifers is a recently growing field, no relevant measured data have been published regarding displacement characteristics of CO2-brine systems at in situ conditions until early 2000.Relevant data were only available for CO2-oil systems for enhanced oil recovery purposes and a handful of measurements for CO2-oil-brine ternary systems.
To fill the knowledge gap, Bennion and Bachu conducted a series of experiments to measure relative permeability and capillary pressure characteristics at in-situ conditions for CO2-brine systems [33,[37][38][39][40][41][42].They studied sandstone, carbonate, shale and anhydrite rock formations in the Alberta basin in central Alberta.These formations are general representatives of the in-situ temperature, pressure, and salinity in entire Alberta basin, and likely for all on-shore North American sedimentary basins.

Case Study on CO 2 Storage Aquifers in Western Canada
In this section, the developed model is used to study potential CO 2 storage sites in central Alberta, Canada.Province of Alberta has the largest greenhouse gas emissions in Canada, with annual emissions close to 274.1 Mt CO 2 in 2015 [34].The province is underlain by the Alberta basin, which is suitable for CO 2 geological sequestration in all parts except for its shallow northeastern corner [35].The largest concentration of CO 2 sources in Alberta is in the Edmonton region where four coal-fired power plants are located near Wabamun Lake, west of Edmonton, which have combined annual CO 2 emissions of more than 30 Mt CO 2 [36].Figure 9 shows the location of major CO 2 emission sources in central Alberta, Canada.CO 2 sequestration in deep saline aquifers in proximity of these power plants is a promising option for reducing CO 2 atmospheric emissions since the deep coal seams and oil and gas reservoirs in local area do not have sufficient capacity for sequestration of CO 2 .Because CO 2 sequestration in geological media and especially in deep saline aquifers is a recently growing field, no relevant measured data have been published regarding displacement characteristics of CO 2 -brine systems at in situ conditions until early 2000.Relevant data were only available for CO 2 -oil systems for enhanced oil recovery purposes and a handful of measurements for CO 2 -oil-brine ternary systems.
To fill the knowledge gap, Bennion and Bachu conducted a series of experiments to measure relative permeability and capillary pressure characteristics at in-situ conditions for CO 2 -brine systems [33,[37][38][39][40][41][42].They studied sandstone, carbonate, shale and anhydrite rock formations in the Alberta basin in central Alberta.These formations are general representatives of the in-situ temperature, pressure, and salinity in entire Alberta basin, and likely for all on-shore North American sedimentary basins.
Table 2 presents a summary of in-situ conditions for the cored intervals from three sandstone and three carbonate formations in the Wabamun Lake area, southwest of Edmonton, Alberta [33,37].The location of the wells from which core samples were taken is shown in Figure 9.For the Wabamun Group, two samples with low and high permeability were evaluated which results in a total of seven different rock sets.Temperatures in Table 2 are evaluated based on the depth of the samples by considering a geothermal gradient of 25 • C/km [37].
The stratigraphic downhole model for strata in the Wabamun Lake area is available in the literature [37].Since site specific data for leakage pathways are not available, we used the same data reported for the storage formations to represent the characteristics of the leakage pathway.Nevertheless, the analysis provided is general and allows application of the developed scaling to a leakage pathway when site specific data are available.Table 2 presents a summary of in-situ conditions for the cored intervals from three sandstone and three carbonate formations in the Wabamun Lake area, southwest of Edmonton, Alberta [33,37].The location of the wells from which core samples were taken is shown in Figure 9.For the Wabamun Group, two samples with low and high permeability were evaluated which results in a total of seven different rock sets.Temperatures in Table 2 are evaluated based on the depth of the samples by considering a geothermal gradient of 25 °C/km [37].* For easy identification in the figures, letters have been assigned to different sample units.Other parameters used in this study are presented in Table 3.All the parameters in Table 3, except the ones with asterisk, are taken from measurements of Bennion and Bachu [37,40].Bennion and Bachu [37] measured relative permeability parameters at reservoir conditions for rock samples of Table 2 and reported absolute permeability, CO 2 and brine viscosity, end-point relative permeability of CO 2 (k max rnw ), residual brine saturation (S wr ), and generated drainage relative permeability curves for CO 2 -brine systems.In another paper, Bennion and Bachu [40] reported the fitted Corey exponents based on their measured relative permeability data for rock samples of Table 2. Corey exponents are listed in Table 3.In another work, Bennion and Bachu [38] reported capillary pressure, interfacial tension and pore size distribution characteristics on a series of carbonate and sandstone formations from Wabamun Lake area in Alberta, together with those formations given in Table 2.We used the reported CO 2 -brine capillary pressure curves, from their study, and curve-fitted to Equation (15).The evaluated capillary entry pressures (p d ) are reported in Table 3.
In Table 3, density of brine is calculated based on Rowe and Chou [43] correlation using temperature, pressure and salt mass fraction for different formations.CO 2 density is calculated using Soave-Redlich-Kwong (SRK) equations of state [44].Total Darcy velocity (v t ) in Table 3 is estimated based on v t = k∆ρg/µ nw , which reflects the maximum buoyancy velocity for CO 2 and is calculated using Darcy's law based on the parameters for each formation.The calculated values in Table 3 are in the range of reported Darcy velocities for CO 2 in the literature [45,46].In addition, although different velocities could happen in the reservoir, this can be a good measure for comparison of different formations.Finally, we considered a length of 500 m and a tilt angle of π/2 in this study.Generally, there are no data on the geometry and direction of the leakage pathways for the formations.We considered identical length and tilt angles for the study of all formations to provide a fair comparison between different storage formations.The calculated values for M, N Ca and N Gr are presented in the last three rows of Table 3.The estimated values of mobility ratio have a wide range, starting at 1.06 for Cooking Lake formation to 8.15 for Wabamun Low Perm formation.Estimated values of N Ca and N Gr are shown in a N Gr − N Ca space in Figure 10.A reference equilibrium state between gravitational and viscous forces is indicated with dotted line at N Gr = 1, and between capillary and viscous forces at N Ca = 1.As shown in Figure 10, all rock samples fall into the lower left sector where viscous forces dominate over capillary and gravity forces.It is seen that N Ca varies three orders of magnitude, while variations of N Gr are more gradual and fall within one order of magnitude for different formations.
As shown in Figure 10, Cooking Lake carbonate has the lowest capillary and gravity numbers of N Ca = 0.0004 and N Gr = 0.0687.The highest effect of gravity forces is in Basal Cambrian sandstone with N Gr = 0.5444, and the highest effect of capillary forces is in Nisku carbonate with N Ca = 0.5260.
It is worth mentioning that we used maximum buoyance velocity in the previous calculations.However, velocity of a migrating CO 2 plume may vary within orders of magnitude.To study the effect of total velocity used in our scaling analysis, we perform an analysis to investigate the effect of velocity on the position of various storage formations in the N Gr − N Ca parameter map.Using different values for the total Darcy velocity leads to different capillary and gravity numbers, but the ratio of these forces remains the same.Figure 11 shows the estimated values for N Ca and N Gr when total Darcy velocity is varied from 10 −10 m/s to 10 −4 m/s for all the formations shown in Table 3.This figure indicates that the ratio of capillary to gravity numbers stays the same as the velocity is changing while various formations have different ratios.In the case of Cooking Lake formation, for instance, while gravity forces are at equilibrium with viscous forces (i.e., N Gr = 1), capillary forces are greatly dominated by viscous forces (N Ca = 0.006).On the other hand, for Nisku formation, when N Gr = 1, capillary number is N Ca = 3.33, which shows the effects of capillary forces are much higher for Nisku formation than that for Cooking Lake formation.It is worth mentioning that we used maximum buoyance velocity in the previous calculations.However, velocity of a migrating CO2 plume may vary within orders of magnitude.To study the effect of total velocity used in our scaling analysis, we perform an analysis to investigate the effect of velocity on the position of various storage formations in the  3.This figure indicates that the ratio of capillary to gravity numbers stays the same as the velocity is changing while various formations have different ratios.In the case of Cooking Lake formation, for instance, while gravity forces are at equilibrium with viscous forces (i.e., ).On the other hand, for Nisku formation, when  3. , which shows the effects of capillary forces are much higher for Nisku formation than that for Cooking Lake formation.

Other than
Ca N and Gr N , relative permeability parameters are important in estimations of rate of leakage.Figure 12a-c shows capillary, gravity and viscous functions for all formations shown in Table 3.For all of the rock samples except Nisku and Ellerslie formations, the effect of capillary and gravity forces is reduced due to the small values of  3. The numbers in brackets is the ratio of N Gr /N Ca for different formations.
Other than N Ca and N Gr , relative permeability parameters are important in estimations of rate of leakage.Figure 12a-c shows capillary, gravity and viscous functions for all formations shown in Table 3.For all of the rock samples except Nisku and Ellerslie formations, the effect of capillary and gravity forces is reduced due to the small values of N Ca and N Gr , and therefore, viscous function, V(S), is the dominant term for estimating saturation and velocity profiles.Mobility ratio as well as Corey exponents of nc and nb determine the behavior of capillary, gravity and viscous functions.As indicated in previous section, a lower mobility ratio and a higher Corey exponent for CO 2 increase the time of breakthrough and the average normalized saturation behind the displacement front.Among the rock samples, Cooking Lake formation has the lowest mobility ratio of M = 1.06 and highest Corey exponent for CO 2 (nc = 5.6), which results in the highest average saturation behind the displacement front (Figure 12c).In addition, because of small mobility ratio, the values of D(S) and G(S) are relatively high for this case, though their effect on velocity profile is insignificant.Despite similar Corey exponents (nc and nb) of Wabamun Low Perm and Cooking Lake formation, mobility ratio is higher and equal to M = 8.15 for Wabamun Low Perm.Further, by comparing these formations (Figure 12a,b), it can be inferred that Wabamun Low Perm has smaller D(S) and G(S) values in comparison to Cooking Lake formation, and the average saturation behind the displacement front is lower for Wabamun Low Perm as well (see Figure 12c).For Basal Cambrian sandstone, mobility ratio is relatively high (M = 6.44), however, the average saturation behind the displacement front is also high, since Corey exponent for CO 2 is nc = 5 (Figure 12c).Therefore, contribution of D(S) and G(S) are negligible due to the value of mobility ratio (Figure 12a,b).Nisku formation has the lowest Corey exponent for CO 2 , nc = 1.1, and even though the mobility ratio for this formation is relatively small and equal to M = 2.09, it has the lowest average saturation behind the displacement front (Figure 12c).Additionally, as shown in Figure 12a,b, contribution of D(S) and G(S) are relatively high for this dataset.Finally, the behavior of capillary, gravity and viscous functions for Wabamun High Perm and Ellerslie formations are very similar, since their mobility ratios and Corey exponents are almost equal.
Figure 13 shows the normalized gas saturation and contribution of different forces on CO 2 dimensionless velocity, for all rock samples at time t D = 0.2.It can be immediately observed that, due to the domination of viscous forces over capillary and gravity forces (small values of N Ca and N Gr ), plume evolution and the dimensionless velocity greatly resemble the viscous function V(S), for all formations, except Nisku and Ellerslie, for which capillary forces are relatively high.Nisku formation has the fastest plume evolution and a countercurrent flow of brine is observed at the inlet due to the effect of capillary forces (Figure 13b,c).Same as Nisku formation, Ellerslie formation does not develop a sharp font, and contribution of capillary forces on velocity causes a small countercurrent flow of brine at the inlet.CO 2 front is relatively sharp for the rest of the formations.The front propagation for Cooking Lake formation is slow compared to other formations, which is due to the small mobility ratio and high nc, as well as small N Ca and N Gr .Shapes of the CO 2 plume and velocity profiles are similar for Basal Cambrian and Wabamun Low Perm formations.This is because of their position on N Gr − N Ca map (Figure 10) and similar values of mobility ratio and nc (Table 3).Again, in the above analysis, we used maximum buoyance velocity.However, velocity of a migrating CO 2 plume may vary within orders of magnitude, which can change the position of the formations in the N Gr − N Ca parameter map.Nevertheless, the scaling analysis provided is general and allows application to specific storage site when site specific data are available.
The average normalized gas saturation versus dimensionless time for all the formations are shown in Figure 14, which are in good agreement with the saturation and velocity profiles discussed earlier.For instance, for the Cooking Lake formation, the breakthrough of CO 2 is delayed due to the less diffusive nature of the displacement front, whereas in the Nisku formation, the effect of capillary forces leads to a diffusive shape front and therefore an earlier breakthrough is expected with lower CO 2 saturation.
for this formation is relatively small and equal to M = 2.09, it has the lowest average saturation behind the displacement front (Figure 12c).Additionally, as shown in Figure 12a,b, contribution of () DS and () GS are relatively high for this dataset.Finally, the behavior of capillary, gravity and viscous functions for Wabamun High Perm and Ellerslie formations are very similar, since their mobility ratios and Corey exponents are almost equal.for all formations, except Nisku and Ellerslie, for which capillary forces are relatively high.Nisku formation has the fastest plume evolution and a countercurrent flow of brine is observed at the inlet due to the effect of capillary forces (Figure 13b,c).Same as Nisku formation, Ellerslie formation does not develop a sharp font, and contribution of capillary forces on velocity causes a small countercurrent flow of brine at the inlet.CO2 front is relatively sharp for the rest of the formations.The front propagation for Cooking Lake formation is slow compared to other formations, which is due to the small mobility ratio and high nc, as well as small  2.
Since the discussed values are dimensionless, further discussion on the real time of breakthrough and the cumulative amount of leaked CO 2 requires site specific data such as residual brine saturation and porosity of the formations.These parameters are not identical for all the formations and taking these parameters into account would influence the interpretation of results, as described in the following.
The real time of breakthrough and cumulative amount of leakage are calculated based on the dimensionless time of breakthrough and the average normalized saturation in Figure 14 and Table 4.The dimensionless time of breakthrough is converted to the real time of breakthrough using Equations ( 9) and (10), and the cumulative amount of leakage per m 2 of the cross-sectional area of the migration pathway is calculated using S × (1 − S wr ) × L × φ, where S wr is the brine residual saturation and S is the average normalized saturation, S = map (Figure 10) and similar values of mobility ratio and nc (Table 3).Again, in the above analysis, we used maximum buoyance velocity.However, velocity of a migrating CO2 plume may vary within orders of magnitude, which can change the position of the formations in the Gr Ca NN − parameter map.Nevertheless, the scaling analysis provided is general and allows application to specific storage site when site specific data are available.CO2 saturation.
Since the discussed values are dimensionless, further discussion on the real time of breakthrough and the cumulative amount of leaked CO2 requires site specific data such as residual brine saturation and porosity of the formations.These parameters are not identical for all the formations and taking these parameters into account would influence the interpretation of results, as described in the following.The real time of breakthrough and cumulative amount of leakage are calculated based on the dimensionless time of breakthrough and the average normalized saturation in Figure 14 and Table 4.The dimensionless time of breakthrough is converted to the real time of breakthrough using Equations ( 9) and (10), and the cumulative amount of leakage per m 2 of the cross-sectional area of the migration pathway is calculated using The maximum Darcy buoyancy velocity considered for each formation has a great impact on estimating the breakthrough time.Cooking Lake, Wabamun High Perm, and Nisku formations have the highest Darcy buoyancy velocities and thus lower estimated breakthrough times of 51, 56, and 40 days, respectively.Darcy buoyancy velocity decreases for Viking, Ellerslie, Basal Cambrian and Wabamun Low Perm formations, respectively, and as can be seen in Figure 15, time of breakthrough wr S  = −  , gives an estimate of the available pore volume for CO2.
In addition, time of breakthrough is highest for Wabamun Low Perm formation, (~308 years), since the estimated Darcy velocity is lowest for this case compared to the other formations, (

 =
).A possible leakage from a storage aquifer through a leakage pathway can be remediated specially if detected early enough.On the other hand, chances of contribution of other trapping mechanisms and therefore reducing the amount of leakage is higher for cases with higher time of leakage, such as Wabamun Low Perm or Basal Cambrian formations.
Celia et al. [47] assessed risk of brine and CO2 leakage through abandoned wells in the same area that we discussed here (as depicted in Figure 9).They simulated 50 years of CO2 injection using a semi-analytical modeling approach over a study area of 2500 km 2 .They concluded that the behavior of the system is dependent on the interplay of formation and fluid properties, maximum injectivity of the formation and the properties of the leaky wells.Generally, lower number of oil and gas wells penetrate the caprocks of deeper formations and this will clearly result in a tradeoff between depth of injection and risk of leakage.However, their simulations imply that the number and design of the injection wells and injectivity of the formations play a critical role on the assessment of leakage.Due to lack of data, they assigned the well permeabilities based on the probability distribution.With their assumption of a single vertical injection well and considering the limitation of injectivity, they concluded that Nisku and the Basal Cambrian formations are the two viable options for CO2 injection.The maximum Darcy buoyancy velocity considered for each formation has a great impact on estimating the breakthrough time.Cooking Lake, Wabamun High Perm, and Nisku formations have the highest Darcy buoyancy velocities and thus lower estimated breakthrough times of 51, 56, and 40 days, respectively.Darcy buoyancy velocity decreases for Viking, Ellerslie, Basal Cambrian and Wabamun Low Perm formations, respectively, and as can be seen in Figure 15, time of breakthrough is delayed for these formations, correspondingly.Porosity times maximum saturation, φ * = (1 − S wr ) × φ, gives an estimate of the available pore volume for CO 2 .
In addition, time of breakthrough is highest for Wabamun Low Perm formation, (~308 years), since the estimated Darcy velocity is lowest for this case compared to the other formations, (v t = 1.3 × 10 −9 ), and at the same time, the amount of leaked CO 2 is lowest for this case, due to the small pore volume available for fluids flow (φ * = 0.03).A possible leakage from a storage aquifer through a leakage pathway can be remediated specially if detected early enough.On the other hand, chances of contribution of other trapping mechanisms and therefore reducing the amount of leakage is higher for cases with higher time of leakage, such as Wabamun Low Perm or Basal Cambrian formations.
Celia et al. [47] assessed risk of brine and CO 2 leakage through abandoned wells in the same area that we discussed here (as depicted in Figure 9).They simulated 50 years of CO 2 injection using a semi-analytical modeling approach over a study area of 2500 km 2 .They concluded that the behavior of the system is dependent on the interplay of formation and fluid properties, maximum injectivity of the formation and the properties of the leaky wells.Generally, lower number of oil and gas wells penetrate the caprocks of deeper formations and this will clearly result in a tradeoff between depth of injection and risk of leakage.However, their simulations imply that the number and design of the injection wells and injectivity of the formations play a critical role on the assessment of leakage.Due to lack of data, they assigned the well permeabilities based on the probability distribution.With their assumption of a single vertical injection well and considering the limitation of injectivity, they concluded that Nisku and the Basal Cambrian formations are the two viable options for CO 2 injection.

Conclusions
In this study, a numerical model is developed to investigate the influence of driving forces on CO 2 plume evolution and breakthrough time during a leakage from a water saturated pathway.The results of this study can be used for estimation of leakage for different potential CO 2 storage sites where formation properties of leakage pathways are available.
Dimensionless terms of D(S), G(S) and V(S) along with dimensionless capillary and gravity numbers are introduced to express the effect and ratio of capillary, gravity and viscous forces.A leakage process is characterized by a two-phase system, in which the non-wetting and less viscous CO 2 displaces the resident brine.Results from this study show that increasing the effect of gravity forces (and therefore N Gr ) leads to an earlier breakthrough of CO 2 .Similar to N Gr , by increasing the effect of capillary forces (N Ca ) a diffusive-like front is established, which also results in an earlier breakthrough.At high values of N Ca , a backward flow of brine is observed as CO 2 migrates upward.For high N Ca , the saturation distribution is diffusive in nature, which results in a gradually leaking CO 2 .In terms of operations, a sharp front leakage may be more problematic than a gradual leakage (diffusive dominated flow) since it is more difficult to take remedial actions to prevent leakage in case of a sudden release of CO 2 .Additionally, it is seen that the effect of increasing N Gr by one order of magnitude on accelerating the breakthrough time, is much less than that of increasing N Ca .
The relative permeability-saturation relationships have shown to be of great influence for plume evolution.Results show that a low relative permeability of CO 2 (low end-point relative permeability of CO 2 or high Corey exponent for CO 2 ) results in a more compact plume evolution.A compact front propagation results in a delayed breakthrough but has the risk of a sudden release of a large volume of leaked CO 2 .
The developed numerical model was applied to different potential storage formations in Western Canada.In the absence of site specific data for leakage pathways, data reported for the corresponding storage formations were used.Nevertheless, the presented analysis provided is general and allows application of the developed scaling to specific leakage pathway when data are available.For these formations, it is seen that N Ca varies within three orders of magnitude, while variations of N Gr are more gradual and fall within one order of magnitude.In addition, due to small values of N Ca and N Gr , plume evolution and the dimensionless velocity greatly resemble the viscous function, V(S), for all formations, except Nisku and Ellerslie, for which capillary forces are relatively high.It is worth mentioning that, for scaling analysis, we used the maximum buoyancy velocity for different formations, however the velocity of a migrating CO 2 plume may vary within orders of magnitude which may change the position of various storage formations in the N Gr − N Ca parameter map.Comparing different formations in Alberta basin revealed that in determining time of breakthrough and the cumulative amount of leaked CO 2 , parameters such as residual brine saturation and porosity of the medium are critical and significantly contribute to the interpretation of results.
The numerical models presented in this paper predict displacement process for a one-dimensional homogenous media.However, the processes induced by CO 2 leakage may be quite complex and may operate on a broad range of space and timescales.For instance, viscous fingering may significantly alter the shape of CO 2 -water displacement front, creates channels with high CO 2 effective permeability and thus earlier leakage of CO 2 .In addition, the leaked CO 2 may undergo phase changes due to temperature variations, which have not been considered in this study.Therefore, further research should consider temperature variation as well as two-and three-dimensional models with detailed representation of leakage pathways heterogeneities for a more realistic representation of the problem.
where subscript i = 1, 2, 3, . . .denotes the block index (position) and ∆z D is the block size.Superscripts n represent the time step index.
For the evaluation of coefficients at the grid-cell boundaries, two approaches have been used.Arithmetic mean is the formulation for evaluation of diffusion terms for D(S) i+1/2 and D(S) i−1/2 .
whereas V(S) and G(S) on the interface blocks are approximated with their values on the upstream block.With these, Equation (A1) could be written as: The numerical procedure that has been used here is depicted in Figure A2.As shown in this figure, Equation (A3) is written for all the grid blocks at every time step which results in a coupled set of algebraic equations in the following matrix form of Ax = b.where α = ∆t D (∆z D ) 2 , β = ∆t D ∆z D and d i = N Ca αD i−1/2 , d i+1 = N Ca αD i+1/2 .These sets of equations need to be solved simultaneously to determine the saturation distribution at each time and for every grid block.After initialization, all the coefficients of matrix A are evaluated at the new time level (S n+1,ν ) where ν is iteration step.Same as matrix A, the coefficients of vector b are evaluated at the new time level, except the saturation at the current time level (S n i ).Now that all the coefficients are defined, saturation at the new time step (S n+1,ν+1 ) can be evaluated iteratively until the convergence condition is achieved, which in this case is ∑(S n+1,ν+1 − S n+1,ν ) 2 ≤ 1 × 10 −10 .
Figures A4 and A5 compare the analytical and numerical solutions for advection-diffusion equation for N Ca = 1 and N Ca = 0.02, respectively.Finally, the analytical solution of the well-known Buckley-Leveret model is compared with the numerical results here.The Buckley-Leveret problem describes immiscible displacement of one phase by another in the absence of capillary and gravity forces [50].Setting N Ca = N Gr = 0 results in Equation (A8), where viscous function, V(S), has the form of fractional flow equation in Buckley-Leveret problem.For comparison between numerical and analytical solutions, nc and nb are set to 2. Figure A6 shows the comparison between analytical and numerical results for three different mobility ratios.Results show that the numerical model reproduces the front location fairly well at both low and large mobility ratios.

Figure 1 .
Figure 1.Pathways for CO 2 migration for a possible leakage through faults (Adapted from Bachu and Celia [13]).

Figure 2 .
Figure 2. Schematic of the problem.
velocity and is defined as the summation of CO2 and brine velocities.In addition,

Figure 2 .
Figure 2. Schematic of the problem.

Figure 6 .
Figure 6.Average saturation versus dimensionless time (t D ) for: (a) Datasets 1-4 (Table1) to show the effect of varying N Ca and N Gr ; and (b) Datasets 1 and 5-7 (Table1) to show the effect of relative permeability parameters M, nc and nb.

Figure 8 .
Figure 8. Influence of (a) gravity forces and (b) capillary forces on breakthrough time with varying velocity from 1 × 10 −9 m/s to 1 × 10 −3 m/s for two mobility ratios of 1 M = (Dataset 5) and 10 M =

Figure 9 .
Figure 9. Location of large stationary CO2 sources in central Alberta, Canada and wells with core samples used in experiments [33].

Figure 9 .
Figure 9. Location of large stationary CO 2 sources in central Alberta, Canada and wells with core samples used in experiments [33].

Fluids 2018, 3 , 32 * 1 GrN = , and between capillary and viscous forces at 1 CaN
x FOR PEER REVIEW 18 of Estimated values (details are available in the text).The calculated values for M , Ca N and Gr N are presented in the last three rows of Table3.The estimated values of mobility ratio have a wide range, starting at 1.06 for Cooking Lake formation to 8.15 for Wabamun Low Perm formation.Estimated values of 10.A reference equilibrium state between gravitational and viscous forces is indicated with dotted line at = .As shown in Figure10, all rock samples fall into the lower left sector where viscous forces dominate over capillary and gravity forces.It is seen that Ca N varies three orders of magnitude, while variations of Gr N are more gradual and fall within one order of magnitude for different formations.

Figure 10 .,
Figure 10.Gravity number ( Gr N ) versus capillary number ( Ca N ) for different formations in Table3.
the total Darcy velocity leads to different capillary and gravity numbers, but the ratio of these forces remains the same.Figure11shows the estimated values for Ca N and Gr N when total Darcy velocity is varied from 10 −10 m/s to 10 −4 m/s for all the formations shown in Table

Figure 10 .
Figure 10.Gravity number (N Gr ) versus capillary number (N Ca ) for different formations in Table3.

Figure 11 .
Figure 11.Gr N versus Ca N for a varying total Darcy velocity ( v t ) ranging from 10 −10 m/s to 10 −4 m/s for different formations ofTable 3. The numbers in brackets is the ratio of Gr Ca NN for different

Figure 11 .
Figure 11.N Gr versus N Ca for a varying total Darcy velocity (v t ) ranging from 10 −10 m/s to 10 −4 m/s for different formations of Table3.The numbers in brackets is the ratio of N Gr /N Ca for different formations.

Figure 13 2 Dt
Figure13shows the normalized gas saturation and contribution of different forces on CO2 dimensionless velocity, for all rock samples at time

1 0
Sdz D .The cumulative amount of leakage is plotted versus time for all formations in Figure 15.As expected, the results are different from what has been expected based on Figure 14.This is essentially due to the different values of residual brine saturation and porosity of different formations.

Figure 13 .
Figure 13.Normalized CO2 saturation and dimensionless CO2 velocity profiles versus dimensionless length at

Figure 13 .
Figure 13.Normalized CO 2 saturation and dimensionless CO 2 velocity profiles versus dimensionless length at t D = 0.2 for rock samples in Table 3: (a) normalized saturation; (b) dimensionless CO 2 velocity, v D ; (c) capillary contribution to CO 2 velocity, v DD ; (d) gravity contribution to CO 2 velocity,v DG ; and (e) viscous contribution to CO 2 velocity, v DV .

Figure 14 .
Figure 14.Average normalized saturation versus dimensionless time ( D t ) for rock samples in Table2.

S
is the brine residual saturation and S is the average normalized saturation, The cumulative amount of leakage is plotted versus time for all formations in Figure 15.As expected, the results are different from what has been expected based on Figure 14.This is essentially due to the different values of residual brine saturation and porosity of different formations.

Figure 14 .
Figure 14.Average normalized saturation versus dimensionless time (t D ) for rock samples in Table2.
same time, the amount of leaked CO2 is lowest for this case, due to the small pore volume available for fluids flow (

Fluids 2018, 3 , 32 Figure A3 .
Figure A3.Comparison of (a) normalized saturation and (b) its average between numerical simulation (red lines) and analytical solution (black lines) for Equation (A5).Numerical simulations are performed for 200 grid blocks and with a time step of 1 × 10 −3 .

Figure A4 .Figure A3 . 32 Figure A3 .
Figure A4.Comparison of (a) normalized saturation and (b) its average between numerical simulation (red lines) and analytical solution (black lines) for Equation (A6), 1 Ca N = .Numerical simulations are performed for 500 grid blocks and with a time step of 1 × 10 −3 .Finally, the analytical solution of the well-known Buckley-Leveret model is compared with the numerical results here.The Buckley-Leveret problem describes immiscible displacement of one phase by another in the absence of capillary and gravity forces[50].Setting 0

Figure A4 .Figure A4 .
Figure A4.Comparison of (a) normalized saturation and (b) its average between numerical simulation (red lines) and analytical solution (black lines) for Equation (A6), 1 Ca N = .Numerical simulations are performed for 500 grid blocks and with a time step of 1 × 10 −3 .Finally, the analytical solution of the well-known Buckley-Leveret model is compared with the numerical results here.The Buckley-Leveret problem describes immiscible displacement of one phase by another in the absence of capillary and gravity forces[50].Setting 0

Figure A5 . 32 Figure A5 .Figure A6 .
Figure A5.Comparison of (a) normalized saturation and (b) its average between numerical simulation (red lines) and analytical solution (black lines) for Equation (A7), N Ca = 0.02.Numerical simulations are performed for 500 grid blocks and with a time step of 1 × 10 −3 .

Figure A6 .Figure A6 .
Figure A6.Comparison of analytical solution of Buckley-Leverett problem (black lines) and numerical results (red lines) for different mobility ratios of (a) M = 1, (b) M = 2.5, and (c) M = 5 at different times.Number of grid blocks is 600, dimensionless time step is 1 × 10 −4 .

Table 1 .
Simulation parameters used in qualitative analysis.

Table 1 .
Simulation parameters used in qualitative analysis.

ID * Unit Lithology Depth (m) Porosity (Fraction) Pressure (MPa) Temp. ( • C) Salinity (mg/Liters)
* For easy identification in the figures, letters have been assigned to different sample units.

Table 4 .
Time of breakthrough and cumulative leakage calculated for different formations of Table3.

bt Total Buoyancy Velocity, v t (m/s) Real time of Breakthrough, t (Years) Average Normalized Saturation, S Cumulative Leakage at Time of Breakthrough, (m 3 /m 2 )
Fluids 2018, 3, x FOR PEER REVIEW 23 of 32is delayed for these formations, correspondingly.Porosity times maximum saturation, *

Table 4 .
Time of breakthrough and cumulative leakage calculated for different formations of Table3.

m/s) Real time of Breakthrough, t (Years) Average Normalized Saturation, S Cumulative Leakage at Time of Breakthrough, (m 3 /m 2 )
Figure 15.Cumulative leakage per m 2 of leakage pathway versus time for rock samples of Table 3.Figure 15.Cumulative leakage per m 2 of leakage pathway versus time for rock samples of Table3.