Energy Fuels 2010, 24, 1431–1440 Published on Web 01/05/2010
: DOI:10.1021/ef9010498
Sharp Front Capturing Method for Carbon Dioxide Plume Propagation during Injection into a Deep Confined Aquifer Yongzhong Liu,* Le Wang, and Bo Yu Department of Chemical Engineering, Xi’an Jiaotong University, Xi’an, Shannxi 710049, China Received September 16, 2009. Revised Manuscript Received December 3, 2009
An in-depth understanding of patterns and evolutional characteristics of carbon dioxide plume has significant implications to optimize operational strategies and manage carbon dioxide during and after injection into saline aquifers for long-term geological storage. A sharp-interface mathematical model that describes the displacement process and its interfacial dynamics of carbon dioxide-brine flow in a deep confined saline aquifer was proposed for predicting the propagation of the carbon dioxide plume. On the basis of a conservative level set method (LSM), the governing equation of interfacial dynamics for capturing time-dependent fronts of the carbon dioxide plume was established. It features a unified relationship between the properties of fluids and the saturations in both the single-phase fluid regions, and the interfacial region in the computational domain was derived by the approximation of linear fraction functions. The model equations numerically solved under different operational conditions after the effects of two adjustable parameters in the interfacial dynamics equation that determine stable and convergent interfaces in the simulations were illustrated and analyzed. A comprehensive comparison to the similarity solutions to a benchmark problem was made to illustrate the reliability and accuracy of the proposed approach. The simulation results show that the maximum migration distances predicted by the proposed approach are shorter than those of the similarity solutions, and the patterns and evolutional characteristics of carbon dioxide plume are intensively dependent upon both injection durations and injection rates at the injection well. The results also show the importance of the effect of vertical mass flux in modeling the migrations of carbon dioxide plume in saline aquifers regardless of injection rates. Moreover, the influences of the underlying assumptions on the interfacial dynamics of the migration of carbon dioxide were analyzed and discussed in detail. The proposed approach and the simulation results will provide insights into the determination of optimal operational strategies and rapid identification of the consequences of carbon dioxide injection into deep confined saline aquifers.
and mineral trapping mechanisms increase, because these mechanisms represent a long-term immobilization of CO2.3 Furthermore, phase-transfer processes and dissolution of CO2 in brine, for example, will also gain importance on a large time scale.4 Accordingly, there are several major concerns arising for CO2 sequestering in deep saline aquifers to be a reliable and safe solution, such as estimation of storage capacity and efficiency, operational strategies during injection, management of injection wells, etc. For a specific saline aquifer, storage and management of CO2 is intimately related to intrinsic geological conditions of formations, operational parameters, and properties of fluids as well. Differences in density and mobility of the invading CO2 and the host brine lead to different patterns of CO2 plume spreading out in the aquifers. Therefore, the shapes and evolutional characteristics of the CO2 plumes have significant implications with regard to the screening storage volume and the pore pressure below the aquitard, optimizing operational strategies, increasing storage efficiency, and evaluating the risky possibilities of CO2 migration or leakage in the case that the remedial actions are to be taken during and after injection.
1. Introduction Storage of carbon dioxide (CO2) in deep saline aquifers, which provides the highest storage capacity, has been suggested as one of the most promising options for the geological storage of CO2 to mitigate greenhouse gas emissions into the atmosphere. It has been estimated that world aquifers could provide a storage capacity of over 1000 billion tons of CO2.1 Specifically, the total storage capacity of CO2 in the saline aquifers in mainland China reaches roughly 14.35 billion tons.2 Carbon dioxide can be stored in deep saline aquifers by three different mechanisms: hydrodynamic trapping, solubility trapping, and mineral trapping. The first mechanism is to trap carbon dioxide into deep aquifers as a supercritical fluid using hydrodynamic principles. The displacement process of multiphase fluid flow reigns during injection. On a short time scale, hydrodynamic trapping is the dominant storage mechanism. Over time, the contributions of residual, solubility, *To whom correspondence should be addressed. Telephone: þ86-2982664752. Fax: þ86-29-82668789. E-mail:
[email protected]. (1) Intergovernmental Panel on Climate Change (IPCC). Climate Change 2007: Mitigation. Contribution of Working Group III to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge, U.K., 2007. (2) Zhang, W.; Li, Y. L.; Xu, T. F.; Cheng, H. L.; Zheng, Y.; Xiong, P. Long-term variations of CO2 trapped in different mechanisms in deep saline formations: A case study of the Songliao Basin, China. Int. J. Greenhouse Gas Control 2009, 3, 161–180. r 2010 American Chemical Society
(3) Intergovernmental Panel on Climate Change (IPCC). Intergovernmental Panel on Climate Change, Special Report on Carbon Dioxide Capture and Storage, Cambridge, U.K., 2005. (4) Ide, S. T.; Jessen, K.; Orr, F. M. Storage of CO2 in saline aquifers: Effects of gravity, viscous, and capillary forces on amount and timing of trapping. Int. J. Greenhouse Gas Control 2007, 1, 481–491.
1431
pubs.acs.org/EF
Energy Fuels 2010, 24, 1431–1440
: DOI:10.1021/ef9010498
Liu et al.
Solutions to the above-mentioned issues need reliable experimental data or field data under operational conditions of injection, in which these data are often, in fact, not available before injection practice. Thus, mathematical models and numerical approaches are indispensable tools for the evaluation of the processes. Furthermore, comprehensive intercomparisons of different models and approaches, both numerical and analytical ones, can also provide an in-depth understanding of the mechanisms of the displacement process under different operating strategies.5 Approaches of capturing patterns and evolution of CO2 plumes through modeling can be generally classified into two categories. One approach is associated with identifying the shape and evolution of the CO2 plume through CO2 saturation profiles in aquifers during injection, which can be termed as a non-sharp front tracking method. In contrast, the spread of the CO2 plume can also be delineated through an apparent interfacial region by a sharp front, which can be accordingly termed as a sharp front tracking method. As for the non-sharp front tracking method, Sasaki and Fujii6 developed a simulation code and procedure to analyze the flow dynamics of CO2 injection in the sub-surface. Similarly, Liu and Smirnov7 proposed a variable saturation model to predict the storage capacity and coal-bed methane recovery. Class and co-workers8 proposed problem-oriented benchmark examples that relate to the CO2 plume evolution and leakage during injection into saline aquifers. They discussed the influences of some assumptions that were widely used in the existing models, such as constant fluid properties, immiscible interface of CO2/brine, isotopic formation, negligible capillary pressure, etc. In contrast, a number of investigations using semi-analytical mathematical approaches, which fall into the category of sharp front tracking methods, have been carried out to characterize the propagation of CO2 in overlying formations. Nordbotten and Celia9,10 presented analytical models to investigate the CO2 plume evolution and leakage through abandoned wells. The variational method was adopted by postulating that the advances of CO2 are subject to their minimal potential energy required to submerge the lighter CO2 into the denser water. The evolution function for sharp fronts was explicitly expressed. Later, they proposed similarity solutions to track the sharp fronts by reducing the partial differential governing equations to a set of ordinary equations,
in which the pressure distributions in both fluids is explicitly included in the solutions.11 Likewise, Dentz and Tartakovsky12 derived closed-form analytical solutions to describe the multidimensional interfacial dynamics for CO2 injection into homogeneous sub-surface formations. Their solutions are readily applied to CO2 geological storage with both displacement in aquifers and coal-bed methane recovery. However, the prediction on the evolution of CO2 plume diverges significantly from the semi-analytical variational solutions. In dealing with the interfacial dynamics in a multiphase fluids flow system, difficulties usually arise owing to singularities of multiphase interfaces. To resolve these difficulties, the most commonly used methods to model two-phase flow problems are the volume of fluid (VOF) method,13 the level set method (LSM),14 the front tracking method,15 and the phase field method.16 The VOF method possesses the main advantage of exactly conservative volumes of the fluids. The interface is, however, represented by a globally discontinuous function, which makes it difficult to track the moving interface accurately. Nevertheless, in LSM, the interface is defined by the zero contour of a signed distance function, viz., the level set function. Because this function is smooth across the interface, it is easier to track the interface as well as to calculate the interfacial curvature with a high order of accuracy. Osher and Sethian14 initiated the LSM to deal with the simulation of interfaces of the two-phase flow problem, which is a first-order partial differential equation in describing mass conservation at interfaces. This approach was extended to solve the problems for motion of soap bubbles and water drops in air,17,18 where large density and viscosity ratios as well as surface tension are encountered. To obtain a smooth interface and avoid discontinuities at the interface, Olsson et al.19,20 modified the conservative equation at the interface by adding a small amount of artificial viscosity to maintain the stability at the interface and further improved the modified method using a re-initialization procedure to preserve the smooth profile of the level set function that is a regularized characteristic function. The primary objectives of this work are (i) to propose a numerical approach on the basis of a conservative LSM to simulate the interfacial dynamics of the patterns and evolutional characteristics of CO2 plumes during injection into saline aquifers and (ii) to investigate the reliability and feasibility of the proposed method through intercomparison between the simulation solutions and the similarity solutions. This paper is organized as follows. In section 2, the displacement process of the CO2/brine system in a deep confined saline aquifer is depicted. In section 3, a general form of the mathematical
(5) Pruess, K.; Garcia, J.; Kovscek, T.; Oldenburg, C.; Rutqvist, J.; Steefel, C.; Xu, T. F. Code intercomparison builds confidence in numerical simulation models for geologic disposal of CO2. Energy 2004, 29, 1431–1444. (6) Sasaki, K.; Fujii, T.; Nilbori, Y.; Ito, T.; Hashida, T. Numerical simulation of supercritical CO2 injection into subsurface rock masses. Energy Convers. Manage. 2008, 49, 54–61. (7) Liu, G. X.; Smirnov, A. V. Modeling of carbon sequestration in coal beds: A variable saturated simulation. Energy Convers. Manage. 2008, 49, 2849–2858. (8) Ebigbo, A.; Class, H.; Helmig, R. CO2 leakage through an abandoned well: Problem-oriented benchmarks. Comput. Geosci. 2007, 11, 103–115. (9) Nordbotten, J. M.; Celia, M. A.; Bachu, S. Injection and storage of CO2 in deep saline aquifers: Analytical solution for CO2 plume evolution during injection. Transp. Porous Media 2005, 58, 339–360. (10) Nordbotten, J. M.; Celia, M. A.; Bachu, S.; Dahle, H. K. Semianalytical solution for CO2 leakage through an abandoned well. Environ. Sci. Technol. 2005, 39, 602–611. (11) Nordbotten, J. M.; Celia, M. A. Similarity solutions for fluid injection into confined aquifers. J. Fluid Mech. 2006, 561, 307–327. (12) Dentz, M.; Tartakovsky, D. Abrupt-interface solution for carbon dioxide injection into porous media. Transp. Porous Media 2008, 51, 1–13.
(13) Scardovelli, R.; Zaleski, S. Direct numerical simulation of freesurface and interfacial flow. Annu. Rev. Fluid Mech. 1999, 31, 567–603. (14) Osher, S.; Sethian, J. A. Fronts propagating with curvaturedependent speed: Algorithms based on Hamilton-Jacobi formulations. J. Comput. Phys. 1988, 79, 12–49. (15) Unverdi, S. O.; Tryggvason, G. A front-traking method for viscous, incompressible, multi-fluid flows. J. Comput. Phys. 1992, 100, 25–37. (16) Villanueva, W.; Amberg, G. Some generic capillary-driven flows. Int. J. Multiphase Flow 2006, 32, 1072–1086. (17) Kang, M.; Merriman, B.; Osher, S. Numerical simulations for the motion of soap bubbles using level set methods. Comput. Fluids 2008, 37, 524–535. (18) Sussman, M.; Smereka, P.; Osher, S. A level set approach for computing solutions to incompressible two-phase flow. J. Comput. Phys. 1994, 114, 146–159. (19) Olsson, E.; Kreiss, G. A conservative level set method for two phase flow. J. Comput. Phys. 2005, 210, 225–246. (20) Olsson, E.; Kreiss, G.; Zahedi, S. A conservative level set method for two phase flow II. J. Comput. Phys. 2007, 225, 785–807.
1432
Energy Fuels 2010, 24, 1431–1440
: DOI:10.1021/ef9010498
Liu et al.
3. LSM for Tracking Sharp Interfaces of CO2-Brine Flow in Confined Saline Aquifers 3.1. General Form of Equations of Two-Phase Flow in Porous Media. During CO2 injection into a confined aquifer, a typical immiscible two-phase flow problem can be described through mass balance and constitutive equations for both CO2 and brine. For each phase in porous media, the mass balance equation is written as21 DðφSi Fi Þ ð1Þ þ r 3 ðFi ui Þ ¼ qi , i ¼ n, w Dt
Figure 1. Schematic of the carbon dioxide-brine flow system in a deep confined acquifer during injection.
where, in a cylindrical coordinate system for an axi-symmetric problem, r = (∂/∂r)er þ (∂/∂z)ez and r 3 = (1/r)(∂/∂r)r þ (∂/∂z), in which er and ez stand for unit vectors at the r and z directions, respectively, φ is porosity, Si is saturation, F is density (kg m-3), t is time (s), ui is the velocity vector (m s-1), and qi stands for source or sink terms (kg m-3 s-1). The subscript n denotes the nonwetting phase (CO2 in this problem), whereas w denotes the wetting phase (brine in this problem). All variables are volumetric averages over a representative elementary volume. The volumetric velocity for each phase can be described by Darcy’s law21 kri K ui ¼ ðrpi -Fi grzÞ, i ¼ n, w ð2Þ μi
model of two-phase flow in porous media for describing the displacement process is given as fundamentals. Then, the front tracking method of the displacement process in a CO2/brine system is established on the basis of LSM, where the interfacial dynamics equation of front tracking is proposed by unifying the governing equations for single-phase fluid flow in saline aquifers. Later, the governing equations are transformed into dimensionless forms by defining a set of dimensionless variables. In section 4, comprehensive comparisons between the proposed solutions and the similarity solutions are made to validate the proposed approach. Moreover, the assumptions adopted in the proposed approach for simulations on the interfacial dynamics of the CO2 plumes during injection are clarified and discussed further.
where κ is the intrinsic permeability tensor (m2), krn and krw are the relative permeabilities for the nonwetting and wetting phases, respectively, μi represents the dynamic viscosity (Pa s), and g is the acceleration due to gravity (m s-2). To solve the above model equations, two closure relations are required. The first one is that the summation of the saturation Si satisfies X Si ¼ 1, i ¼ n, w ð3Þ
2. Problem Definition As an invading fluid, CO2 is radially injected into a saline aquifer confined by two parallel aquitars through an injection well with radius r0, as shown in Figure 1. The horizontal aquifer is homogeneous and isotropic with constant thickness H, porosity φ, and permeability κ. The injected CO2 plume progressively replaces the resident brine and spreads over the aquifer. For the geological conditions of typical reservoirs, it will be naturally at a supercritical state. The density ratio between the two fluids ranges from 0.22 to 0.75, and the viscosity ratios range from 0.026 to 0.20.9 The differences between the density and viscosity of the two fluids cause their stratification, with the invading CO2 laying above the host brine. The patterns and evolutional characteristics of the CO2 plume in the aquifer are, therefore, dependent upon the operational conditions and behaviors of the interaction between the two-phase fluids flow in the porous aquifer. To capture the evolutional characteristics of the CO2 plume in the aquifer and reduce the complexity of the mathematical model, the following assumptions and simplifications are made: (1) During injection, a time-dependent sharp front with a narrow and fixed thickness, which denotes a mushy transitional region of the two phases, separates the invading CO2 and the resident brine in the aquifer. Away from the interfacial region, the single-phase CO2 zone is fully occupied with CO2, whereas the single-phase brine zone is fully saturated with brine. (2) In comparison to the characteristic length scale of the aquifer, the magnitude of the thickness of the sharp front is small enough that the influence on the solutions could be ignored. (3) The properties of the interfacial region where CO2 and brine co-exist, such as density, viscosity, etc., are continuous functions of the saturation. (4) CO2 and brine are two separate and immiscible phases. The effect of mutual dissolution is neglected. (5) The capillary pressure is negligible. (6) Fluid properties, such as density and viscosity, are constant.
i
The other one is usually the relationship between pressures of two phases and capillary pressure, namely, pc ¼ pn -pw ð4Þ When the above equations are used to describe the displacement process in a CO2-brine system, the relative permeability kri is usually considered a nonlinear function of the saturation Si that is correlated to the capillary pressure pc.8,22,23 However, if the capillary pressure pc is to be ignored, which is usually valid for homogeneous porous media,8,24,25 it is no longer appropriate to associate the relative permeability kri with the capillary pressure pc. This implies that, to solve the equations, an additional relationship between the relative permeability kri and the saturation Si is needed to satisfy the closure condition. We will discuss this issue in detail and deduce an expression of the relationship in the next sections. 3.2. Interfacial Dynamics Equation for Tracking a Sharp Interface through the Conservative LSM. On the basis of the general form of equations of two-phase flow in porous media, an interfacial dynamics equation that is usually characterized by the saturations of fluids should be included to capture the (21) Bear, J. Dynamics of Fluids in Porous Media; Dover Publications, Inc.: New York, 1972. (22) Niessner, J.; Helmig, R. Multi-scale modeling of three-phase-three-component processes in heterogeneous porous media. Adv. Water Resour. 2007, 30, 2309–2325. (23) Suekane, T.; Ishii, T.; Tsushima, S.; Hirai, S. Migration of CO2 in porous media filled with water. J. Therm. Sci. Technol. 2006, 1, 1–11. (24) Fagerlund, F. F.; Niemi, A.; Oden, M. Comparison of relative permeability-fluid saturation-capillary pressure relations in the modelling of non-aqueous phase liquid infiltration in variably saturated, layered media. Adv. Water Resour. 2006, 29, 1705–1730. (25) Plug, W. J.; Bruining, J. Capillary pressure for the sand-CO2water system under various pressure conditions. Application to CO2 sequestration. Adv. Water Resour. 2007, 30, 2339–2353.
1433
Energy Fuels 2010, 24, 1431–1440
: DOI:10.1021/ef9010498
Liu et al.
possible to obtain a sharp interface in simulations. Olsson and Kreiss19 have proved that a constant thickness and smooth front can be achieved, provided that a set of proper parameters is chosen, and convergent solutions could reach up to secondorder accuracy and very good conservation of the area bounded by the interface. Numerical stability and mass conservation at the front could be further improved by the modification of the re-initialization step.20 Consequently, eq 7 can be used to describe the movement of the interfacial region that separates CO2 and brine in a confined aquifer. 3.3. Governing Equations for the Displacement Process of CO2/Brine in a Confined Aquifer. 3.3.1. Governing Equations for the Two-Phase Flow and the Interfacial Dynamics. For the displacement process of CO2/brine in a confined aquifer, the mass conservation equation for the two-phase flow can be obtained by the addition of eq 1 for the two phases, namely, DðφFÞ þ r 3 ðFuÞ ¼ q ð8Þ Dt
evolutional characteristics of the interfacial region. That also implies that an additional relationship that links the saturation and the relative permeability should be involved if the model equations are solvable. For the relationship between the saturation and the relative permeability, it is worthy of noting that the relative permeability kri should reach its two extremes in the region that is solely filled with a single phase, namely, kri ¼ f ðSi Þ ¼ 0 Si ¼ 0 , i ¼ n, w ð5Þ 1 Si ¼ 1 This means that the relative permeability kri should be unity in the region that is fully occupied by one fluid phase, whereas it should be zero in the counterpart domain. In the interfacial region that separates the two phases of immiscible fluids, however, it should be a certain value of zero to unity. Furthermore, according to the sharp front assumption, the interfacial region is spatially narrow. This implies that the forms of functions kri = f(Si) will not seriously affect the patterns of the interfacial region provided that they are continuous functions. To establish the interfacial dynamics equation in capturing the evolutional characteristics of the carbon dioxide plume, in this paper, we use the LSM, which has been successfully used to describe the interfacial dynamics for two-phase flow. For a twophase flow in porous media, the movement behaviors of fluids can be described by DðφSÞ þ u0 3 rS ¼ 0 ð6Þ Dt
where F = SFw þ (1 - S)Fn and q = qw þ qn = 0. Likewise, the constitutive equation for fluids flows in both single-phase regions and the interfacial region can also be obtained by the addition of eq 2 for the two phases, that is K F ð1 -SÞ K F S u ¼- 3 n ðrpn -Fn grzÞ - 3 w ðrpw -Fw grzÞ μn F μw F ð9Þ It can be seen that eqs 8 and 9 have the forms that are the same as the equations describing single-phase fluid flow in porous media. Therefore, the equations represent the flow of the CO2 phase when S = 0, whereas the equations represent the flow of the brine phase when S = 1, and they can also represent the twophase flow in the interfacial region when 0 < S < 1. On the basis of the assumption made in this paper, the capillary pressure is neglected in this problem. Thus pn ¼ pw ð10Þ
where S represents the saturation of the wetting phase, whereas the saturation of the nonwetting phase is, therefore, (1 - S). In the region that is fully occupied by a single fluid, u0 represents the substantial velocity of each phase. In the interfacial region in which two fluids co-exist, however, u0 represents the velocity of the interfacial movement, which is neither the velocity of the nonwetting phase nor the velocity of the wetting phase. In establishing the interfacial dynamics equation, one of the key problems to deal with is to depict the movement of the interface, i.e., the expression of u0 . Note that the substantial velocities of each fluid phase in the interfacial region are different. To describe the movement behaviors of the interfacial region, it is imperative to integrate the substantial velocities of each phase into the interfacial velocity u0 . That is to say that the velocity of the interfacial movement should be unified with the substantial velocities of each phase in the interfacial dynamics equation. When eq 6, which is a mass conservation equation, is used to delineate the movement of the interfacial region, it is a nonsharp front model for the interface region. The form of the equation makes it, in practice, not feasible to obtain stable numerical solutions in capturing the interfacial behaviors. To ensure stable numerical solutions and to delineate the evolution of a sharp front, an improvement can be made by adding an artificial diffusion term on the right-hand side of eq 6.26 Then, it is written as19,20 DðφSÞ rS þ u0 3 rS ¼ γr 3 εrS -Sð1 -SÞ ð7Þ Dt jrSj
Consequently, the interfacial dynamics equation for the interfacial region, eq 7, is then written as DðφSÞ rS þ u 3 rS ¼ γr 3 εrS -Sð1 -SÞ ð11Þ Dt jrSj where, in the interfacial region, the interfacial velocity u is unified with the substantial velocities of each phase by eq 9. On the basis of the sharp front assumption, the properties of fluids in the interfacial region, such as density, viscosity, and relative permeability, are usually treated to be stepwise functions, for example, Heaviside functions. However, it is difficult to implement in numerical simulations, and it may also cause front tracking to be instable in simulations. To guarantee the continuity of the properties of fluids in the interfacial region, we use a linear function (the density in eq 8) and linear fraction functions (the mobility in eq 9) to approximate the relationship between the properties of fluids and the saturation Si in the interfacial region. These approximations can effectively deal with the problems of stability and smoothness in the interfacial region, notwithstanding that the suggested linear fraction approximation may deviate from the Van Genuchten model to a certain extent. In addition, the main advantage of this approximation is that the explicit relationship between the saturation and the capillary pressure is no longer necessarily involved in the equations. 3.3.2. Initial Conditions and Boundary Conditions. To associate operational conditions of injection with boundary conditions, the flow rate of CO2 in the injection well, Qwell, is expressed as Qwell ¼ 2πr0 Hu0 ð12Þ
The introduced artificial term is to make the shape of the interface stable in computation and the thickness of the interfacial region controllable by two adjustable parameters: γ, which is a parameter to stabilize the front in simulations, and ε, which is a parameter to control the thickness of the interfacial region in simulations. In a meaningful extent of the solutions to the problem, the smaller the parameter ε, the narrower the interfacial region. Therefore, it should be a value as small as (26) Deshpande, K. B.; Zimmerman, W. B. Simulation of interfacial mass transfer by droplet dynamics using the level set method. Chem. Eng. Sci. 2006, 61, 6486–6498.
where r0 denotes the radius of the injection well (m) and u0 denotes the velocity at the wall of the injection well (m/s). 1434
Energy Fuels 2010, 24, 1431–1440
: DOI:10.1021/ef9010498
Liu et al.
variables and parameters, the mass conservative equation, eq 8, is accordingly converted into
As shown in Figure 1, the initial conditions and the boundary conditions for eqs 8-11 are as follows: Initial conditions
DF ~ þ r 3 ðF uÞ ¼ q Dt
p ¼ Fw gðH -zÞ; S ¼ 1 at r, z Boundary conditions r ¼ r0 : n 3 u ¼ u 0 ; S ¼ 0
where, in a cylindrical coordinate system for an axi-symmetric ~ = (1/r)(∂/∂r)r þ (∂/∂z), ~ = (∂/∂r)er þ (∂/∂z)ez and r problem, r 3 in which er and ez stand for unit vectors at the r and z directions, respectively. Thus, the constitutive equation, eq 9, is transformed into ~ -β rzÞ ~ u ¼ -Rðrp ð22Þ
z ¼ 0 : n 3 u ¼ 0; n 3 rS ¼ 0 r ¼ R : p ¼ Fw gðH -zÞ; S ¼ 1 z ¼ H : n 3 u ¼ 0; n 3 rS ¼ 0
where
For tracking the patterns and evolutional characteristics of the CO2 plume in the confined aquifer, one can substitute eq 9 into eqs 8 and 11, respectively, and solve the coupled equations with the corresponding initial conditions and boundary conditions to calculate the saturation S. In the simulations, the saturation S = 0.5 is considered the front of the CO2 plume in the aquifer, as usual.19 3.4. Non-dimensionalization of the Model Equations. 3.4.1. Definition of Dimensionless Parameters. It is not unusual that instable solutions are encountered in the simulations as a result of a large aspect ratio in the computational domain of interest. To adapt the computations for different geometry scales of the aquifers and different time scales as well, we conducted a non-dimensionalized transformation with respect to eqs 8-11. For a cylindrical coordinate system, we defined a group of dimensionless variables and parameters as follows: (1) Dimensionless spatial and temporal variables r ð13Þ r ¼ H z H
ð14Þ
t t ¼ t
ð15Þ
z ¼
p ¼
p -p0 Fw gH
or written as
82 3 ! 1
