Accepted Manuscript

Numerical simulation of displacement characteristics of CO2 injected in pore-scale porous media

Qianlin Zhu, Qianlong Zhou, Xiaochun Li

PII: S1674-7755(15)00107-9

DOI: 10.1016/j.jrmge.2015.08.004

Reference: JRMGE 194

To appear in: Journal of Rock Mechanics and Geotechnical Engineering

Received Date: 18 May 2015 Revised Date: 18 August 2015 Accepted Date: 24 August 2015

Rock Mechanics and

Geotechnical

Engineering

—- —

Please cite this article as: Zhu Q, Zhou Q, Li X, Numerical simulation of displacement characteristics of CO2 injected in pore-scale porous media, Journal of Rock Mechanics and Geotechnical Engineering (2015), doi: 10.1016/j.jrmge.2015.08.004.

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A A/I A \TT TO TTIT

Journal of Rock Mechanics and Geotechnical Engineering

Journal online: www.rockgeotech.org

Numerical simulation of displacement characteristics of CO2 injected in pore-scale porous media

Qianlin Zhua b, *, Qianlong Zhoua, Xiaochun Lib

a. Key Laboratory of Coal-based CO2 Capture and Geological Storage, China University of Mining and Technology, Xuzhou, 221008, China

b. State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan, 430071, China

Received 18 May 2015; received in revised form 18 August 2014; accepted 24 August 2015

Abstract: Pore structure of porous media, including pore size and topology, is rather complex. In immiscible two-phase displacement process, the capillary force affected by pore size dominates the two-phase flow in the porous media, affecting displacement results. Direct observation of the flow patterns in the porous media is difficult, and therefore knowledge about the two-phase displacement flow is insufficient. In this paper, a two-dimensional (2D) pore structure was extracted from a sandstone sample, and the flow process that CO2 displaces resident brine in the extracted pore structure was simulated using the Navier-Stokes equation combined with the conservative level set method. The simulation results reveal that the pore throat is a crucial factor for determining CO2 displacement process in the porous media. The two-phase meniscuses in each pore throat were in a self-adjusting process. In the displacement process, CO2 preferentially broke through the maximum pore throat. Before breaking through the maximum pore throat, the pressure of CO2 continually increased, and the curvature and position of two-phase interfaces in the other pore throats adjusted accordingly. Once the maximum pore throat was broken through by the CO2, the capillary force in the other pore throats released accordingly; subsequently, the interfaces withdrew under the effect of capillary fore, preparing for breaking through the next pore throat. Therefore, the two-phase displacement in CO2 injection is accompanied by the breaking through and adjusting of the two-phase interfaces. Key words: level set method; displacement; pore-scale porous media; numerical simulation

1. Introduction

The increase of greenhouse gas emissions to the atmosphere during the 20th century is considered to have contributed to global warming. A reduction in the released rate of CO2 to the atmosphere is essential for mitigating greenhouse effects. One way of achieving this is to inject CO2 into deep formations (Yang et al., 2015). Such formations include oil and gas reservoirs, saline aquifers and unminable coal seams, of which the deep saline aquifers are widely distributed with the largest CO2 storage capacity (IPCC, 2005). In China, the storage capacity of deep saline aquifers accounts for 88% of the total geological storage capacity (Li et al., 2009).

CO2 that is injected into a saline aquifer displaces the resident brine. Different to single-phase flow, the flow process that CO2 displaces brine is an immiscible two-phase flow, in which the capillary force between the CO2 and brine plays an important role. The capillary force is inversely proportional to interface curvature and directly proportional to interface tension. Pore diameter of rock is commonly smaller than millimeter-scale, and the rock is commonly hydrophilic. Therefore, the capillary force is an important factor that may control the displacement of the immiscible two-phase flow. However, tempo-spatial variation of two-phase interfaces between the immiscible two-phase fluids is difficult to be observed in porous media (Liu et al., 2012). In addition, methods based on Darcy's law are effective in investigating macroscale two-phase flow but cannot reveal the flow dynamics in the immiscible displacement (Christensen, 2006). Therefore, the dynamics of the immiscible two-phase flow process needs to be understood from a pore-

scale viewpoint, which is important but very complicated due to a large number of factors influencing the flow, such as fluid density, fluid viscosity, surface tension, flow rates, surface wettability, pore geometry and medium heterogeneity (Pan, 2003; Liu et al., 2013).

Several approaches have been applied to simulate multiphase flow at pore-scale in porous media, in which the pore-network model is a conventional method. Based on two-dimensional (2D) model that is simplified by ball and column connection, Lenormand and Touboul (1988) studied the immiscible two-phase flow using the pore-network model, and established a phase-diagram to identify the flow patterns. In this phase-diagram, the flow patterns were identified by two dimensionless parameters (i.e. capillary number and viscosity ratio) that are determined by fluid viscosity, flow velocity, contact angle and interfacial tension. The capillary number does not refer to the pore structure. Therefore, the phase-diagram is difficult to be applied but plays an instructive role in analysis of two-phase flow in many natural porous media. Blunt et al. (2002) reviewed the pore-network simulation in pore-scale flow. The pore-network model simplifies porous media by connection model of pore body and pore throat with idealized geometries, acquiring certain prediction effects. However, it is difficult to extract a pore-network that topologically and geometrically represents the complex pore geometry and physics, and therefore the pore-network model cannot reveal the fluid dynamics mechanism (Ramstad et al., 2012). In addition to pore-network model, computational fluid dynamics methods, including the lattice Boltzmann method (Liu et al., 2012), level set method (Gunde et al., 2010), volume-of-fluid method (Raeini et al., 2012) and phase field method (Anderson et al., 1998), have been developed in pore-scale simulation.

Corresponding author. Tel: +86 56 83883501; E-mail: zhuql@cumt.edu.cn

Simulating pore-scale two-phase flow by computational fluid dynamics methods can reveal flow characteristics.

In this paper, a computational fluid dynamics method with the level set method is used to simulate movement characteristics of the two-phase interfaces between the CO2 and resident brine in the displacement process at pore-scale, and the displacement mechanisms are analyzed, which may improve our knowledge on CO2 displacement flow in porous media.

2. Methodology

2.1. Geometric model

A sandstone sample is chosen from a shallow saline aquifer at a depth of about 200 m, which is located about 41 km northeast of Tongliao, Inner Mongolia, China. In this saline aquifer, a small-scale CO2 shallow injection experiment is conducted (Zhu et al., 2015). The pore structure (Fig. 1) is extracted from a microgram of sandstone slice.

Fig. 1. Pore structure extracted from a sandstone slice.

2.2. Numerical method

According to the continuum assumption of fluid, a crucial issue for numerically simulating immiscible two-phase displacement is to track the movement of interfaces in the two-phase flow. The level set method, an interface capture method, was firstly proposed by Osher and Sethian (1988). In this method, the flow field is covered by a level set function, and the interfaces are expressed with contour lines or surfaces of the level set function. The level set function moves with the flow just as the interfaces do, and the level set function remains a smooth function when it evolves. The level set function performs naturally for problems with singularities, large distortion or interface breakup, etc. Unlike the volume-of-fluid method, the level set function is a continuous function. Therefore, there are more choices for numerical discretization. In addition, some geometric parameters such as the interface normal vector, tangent vector and curvature can be derived from simple derivation.

In standard level set method, the level set function f(X, t) is defined to be a signed distance function:

f (X, t)| = d(x) = min(|x - x |) (1)

where X is the point coordinates, t is the time, i is the interface domain, f(X, t) >0 on one side of the interface, f(X, t) <0 on the other side, and f (X, t) =0 on the interface. The level set function f(X, t) evolves under advection velocity field as

ft + u-V f = 0 (2)

where u is the velocity field, and therefore is the movement velocity on the interface; the subscript t, which denotes the time, represents a partial derivative with respect to time; V is the Laplace operator. In a

simulation that uses fluid velocity u, the level sets close to the zero level set move with velocities that considerably distort the level set function. In this case, the level set function will cease to be an exact signed distance function after a few of time steps, affecting numerical calculation at the next time step (Sussman and Puckett, 2000). Therefore, a process called re-initialization is needed, which is used to stop the level set calculation at some points in time and to rebuild a level set function corresponding to the signed distance function. There are several ways to accomplish this re-initialization. The most commonly used method among them is the method introduced by Sussman et al. (1994). Its virtue is that the level set function is reinitialized without explicitly finding the zero level set. This reinitialization method is used to solve Eq. (3) to steady state (t ® ¥ ):

2 + h2

where f0 is the level set function distribution calculated by Eq. (2) before the re-initialization; h is a parameter of amplitude that is a function of grid size; the subscript t denotes a partial derivative with respect to the re-initialization time t that is different from the time t. Calculating Eq. (3) to steady state will provide a new value for f with a property that |Vf| = 1 . In addition, this re-initialization calculation ensures that the zero level set does not move (Sussman et al., 1994). However, a main defect of this level set method is that the advection of f , including the re-initialization steps, is not done in a conservative way, not even for divergence-free velocity fields. This implies that the area bounded by the zero level set is not conserved (Sussman and Puckett, 2000; Olsson and Kreiss, 2005). In order to ensure mass conserved, Olsson and Kreiss (2005) and Olsson et al. (2007) modified this level set function by replacing the signed distance function with a smeared-out Heaviside function f (X, t):

0 [f(X, t) <-e]

f(x, t) =

1 + f + ^sin f [-£<f(x, t) <e]

2 2e 2n e L ' J

[f(x, t) > e]

where e corresponds to half of the interface thickness, which can be taken from half of the characteristic dimension of the smallest cell in the mesh of computational geometry (Gunde et al., 2010). After the modification, the interface location is represented by the contour f (X, t) = 0.5. In this paper, this modified level set method is chosen to be the interface tracking method in the CO2 displacement simulation. The normal vector n and curvature k of the interface can easily be obtained from the level set function field as (the subscript of modified function f will be omitted in the following part): Vf |

f=0 5 I

||f=0.5

k = -V-1

|f=05 I

The momentum equation of fluid is described by the incompressible Navier-Stokes equation as

put + p(u ■ V)u = V ■ {-pI+h[Vu+(Vu)T]} + Fsv (6)

and the continuity equation as

V ■ u = 0 (7)

where p is the fluid density, n is the dynamic viscosity, I is the identity tensor, p is the pressure field, and Fsv is the term representing interfacial tension. The level set function evolving under the effect of fluid velocity can be described as

f + u ^0 =—V-[eVf-f(1 -f)n]

where m denotes the amount of re-initialization or stabilization of the level set function f , and e determines the interface thickness between the two fluids. The density p and viscosity n are related to the level set f in the region of the two-phase flow: P = Pl + (Pg ~Pl)f h = h + (hg -hl)f

where the subscripts l and g refer to the liquid and gas phases, respectively. The interfacial tension Fsv is calculated according to the continuum surface tension model proposed by Brackbill et al. (1992) as follows:

F„(x) = a8(-V- n)n (10)

where a is the surface tension coefficient; and 8 is the Dirac delta function, which can be written as

8 = 6 |Vf||f(1 -f)\ (11)

This function ensures the following equation:

J0 Fsvdf = an (12)

In addition, the continuum surface tension model is independent of grid size.

3. Numerical solution and discussion

The left and right sides of the porous media shown in Fig. 1 are taken as inlet and outlet boundary conditions, respectively; the lower and upper boundaries are taken as a symmetric boundary condition. Pores are considered as water wettable with a contact angle of 180°. The CO2 injected into porous media will automatically choose flow channels. Therefore, a buffer region, where the CO2 can automatically choose flow channels, is connected to the inlet boundary of the geometric model. The modified geometric model is shown in Fig. 2.

Fig. 2. Porous media with a buffer region.

The inlet velocity is set to 0.5 cm/s. According to conditions of the shallow saline aquifer, the density and viscosity of the water are set to 1 x 103 kg/m3 and 1.34 x10 3 Pa s, respectively. The density and viscosity of the CO2 are set to 38.5 kg/m3 and 1.47x10-5 Pa s, respectively. The interfacial tension between water and CO2 is set to 44.57x 10-3 N/m (Chalbaud et al., 2010).

The finite element method is used to solve the coupled equation system. The geometric model is discretized with 2D unstructured triangle meshes, and boundaries with large curvature are discretized with finer meshes as shown in Fig. 3. The number of the mesh is 251,327. A linear system solver PARDISO (parallel direct solver) is used. A relative tolerance of 0.01 and an absolute tolerance of 0.1 are assumed for the time discretization scheme. The time-stepping method uses the generalized alpha method to determine each time step.

An initial condition given by an assumed phase distribution may not be reasonable. Therefore, before the displacement simulation, the flow field is calculated with zero inlet velocity within a short time to make the assumed phase distribution to reach a reasonable pattern. The calculated initial phase distribution is shown in Fig. 4.

At the beginning of the injection, the injected CO2 displaces the resident water in the porous media through a relatively large channel, as shown by the arrow in Fig. 5. At this time, the CO2 in relatively narrow pores does not show obvious flow. When the CO2 in the relatively large channel reaches a narrow pore throat as shown by the arrow in Fig. 6a, the CO2 will stop displacing through this channel, and the CO2 in relatively narrow pores begins to displace the water as shown by the arrows in Fig. 6b.

Fig. 3. Pore geometry discretized with unstructured 2D meshes.

Fig. 4. Initial CO2-water distribution.

Fig. 5. Injected CO2 flows through the relatively large pore.

Comparing CO2 distribution in Fig. 6b and c, we may see that although CO2 firstly flows through a relatively large pore, a relatively narrow throat restricts the CO2 to flow through this channel. Therefore, pore throat, rather than pore body, is determinant that determines channel of the CO2 flow. This result indicates that capillary force is determinant in CO2 flow in porous media.

Distribution of the two-phase interface in Fig. 6c shows that, when reaching a pore throat, CO2 will displace through another channel. Consequently, there will be a moment when all two-phase interfaces reach respective pore throats of potential flow channels. At this time, once the maximum pore throat is broken through by the CO2, meniscuses at the other pore throats will retreat quickly; this process can be shown by comparing Fig. 6c and d. Therefore, the two-phase interfaces are in a self-adjusting process in the displacement process. This result also indicates that CO2 may have ever reached a position in porous media where although the CO2 is not seen at a certain time, such as the positions shown with pink arrows in Fig. 6d. This push-pull movement of interfaces may improve the dissolution and reaction of CO2 in porous media (Li, 2011).

Time 0.059 s

Fig. 6. Variations of two-phase

Fig. 7 shows the variation of injection pressure during the displacement process. The pressure variation indicates that injection pressure fluctuates throughout the displacement process. The variations of two-phase distribution and injection pressure reveal that, when CO2 reaches a pore throat and needs to break through it, the injection pressure needs to increase high enough. Once the CO2 breaks through the pore throat under effect of the increased injection pressure, the injection pressure will release quickly, and the stored pressure by the capillary pressure of meniscus will release with the meniscuses retreating. Subsequently, the injection pressure gradually increases to prepare for breaking through the next pore throat.

Time (s)

Fig. 7. Variation of injection pressure in the injection process.

Time 0.063 s

interface during displacement.

Fig. 8 shows the finally simulated results of two-phase distribution. This result shows that 47.71% of the water is displaced. Therefore, more than half of the water remains in the porous media.

Fig. 8. Simulated results of CO2 distribution.

4. Conclusions

A 2D pore structure is extracted from a sandstone sample, and CO2 displacement in the sandstone is simulated using the Navier-Stokes equation combined with the conservative level set method. The simulation results show that pore throat, rather than pore body, is a crucial factor that determines flow channel of the CO2 displacement in porous media. The injected CO2 firstly flows through the largest pore. When reaching a pore throat, the pore throat will restrict the CO2 from flowing in that channel, and the CO2 will choose to flow through the second-largest pore, and so on. In this way, CO2 displacement needs to break through one or more pore throats. At this time, the CO2 will choose the relatively large pore throat to break through, and therefore this pore throat determines flow channel of the CO2.

The two-phase meniscuses in each pore throat are in a self-adjusting

process during the displacement. Before breaking through the relatively large pore throat, the injection pressure of the CO2 increases gradually, and the curvature and the position of interfaces in other pore throats adjust accordingly. Once the maximum pore throat is broken through by the CO2, the capillary force in other pore throats releases subsequently, and the interfaces retreat accordingly, preparing for breaking through the next pore throat. The two-phase displacement in CO2 injection is accompanied by the breaking through and adjusting of the two-phase interfaces.

Conflict of interest

The authors wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgements

This work is funded by Key Laboratory of Coal-based CO2 Capture and Geological Storage, Jiangsu Province, China and US Advanced Coal Technology Consortium (No. 2013 DFB60140-08). We also thank Prof. Qi Li for useful discussions.

References

Anderson DM, McFadden GB, Wheeler AA. Diffuse-interface models in fluid

mechanics. Annual Review of Fluid Mechanics 1998; 30(1): 139-65. Blunt MJ, Jackson MD, Piri M, Valvatne PH. Detailed physics, predictive capabilities and macroscopic consequences for pore-network models of multiphase flow. Advances in Water Resources 2002; 25(8-12): 1069-89. Brackbill JU, Kothe DB, Zemach C. A continuum method for modeling surface tension.

Journal of Computational Physics 1992; 100(2): 335-54. Chalbaud C, Robin M, Lombard JM, Bertin H, Egermann P. Brine/CO2 interfacial properties and effects on CO2 storage in deep saline aquifers. Oil & Gas Science and Technology 2010; 65(4): 541-55. Christensen BSB. Using X-ray tomography and lattice boltzmann modeling to evaluate pore-scale processes in porous media. PhD Thesis. Lyngby, Denmark: Institute of Environment & Resources, Technical University of Denmark, 2006. Gunde A, Bera B, Mitra SK. Investigation of water and CO2 (carbon dioxide) flooding using micro-CT (micro-computed tomography) images of Berea sandstone core using finite element simulations. Energy 2010; 35(12): 5209-16.

Intergovernmental Panel on Climate Change (IPCC). IPCC special report on carbon dioxide capture and storage. Cambridge, UK: Cambridge University Press, 2005.

Lenormand BR, Touboul E. Numerical model and experiments on immiscible displacements in porous media. Journal of Fluid Mechanics 1988; 189: 165-87.

Li Q. Coupled reactive transport model for heat and density driven flow in CO2 storage in saline aquifers. Journal of Hazardous, Toxic, and Radioactive Waste 2011; 15(4): 251-8.

Li XC, Fang ZM, Wei N, Bai B. Discussion on technical roadmap of CO2 capture and storage in China. Rock and Soil Mechanics 2009; 30(9): 1674-8 (in Chinese).

Liu HH, Valocchi AJ, Kang QJ. Pore-scale simulation of high-density-ration multiphase flows in porous media using lattice Boltzmann method. In: Proceedings of the 19th International Conference on Computational Methods in Water Resources. University of Illinois at Urbana-Champaign, 2012.

Liu HH, Valocchi AJ, Kang Qj, Werth C. Pore-scale simulations of gas displacing liquid in a homogeneous pore network using the lattice Boltzmann method. Transport in Porous Media 2013; 99(3): 555-80.

Olsson E, Kreiss G. A conservative level set method for two phase flow. Journal of Computational Physics 2005; 210(1): 225-46.

Olsson E, Kreiss G, Zahedi S. A conservative level set method for two phase flow II. Journal of Computational Physics 2007; 225(1): 785-807.

Osher S, Sethian JA. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics 1988; 79(1): 12-49.

Pan CX. Use of pore-scale modeling to understand flow and transport in porous media. PhD thesis. Chapel Hill, USA: University of North Carolina, 2003.

Raeini AQ, Blunt MJ, Bijeljic B. Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method. Journal of Computational Physics 2012; 231(17): 5653-68.

Ramstad T, Idowu N, Nardi C, 0ren PE. Relative permeability calculations from two-phase flow simulations directly on digital images of porous rocks. Transport in Porous Media 2012; 94(2): 487-504.

Sussman M, Puckett EG. A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows. Journal of Computational Physics 2000; 162(2): 301-37.

Sussman M, Smereka P, Osher S. A level set approach for computing solutions to incompressible two-phase flow. Journal of Computational Physics 1994; 114(1): 146-59.

Yang DX, Li Q, Zhang LZ. Propagation of pore pressure diffusion waves in saturated porous media. Journal of Applied Physics 2015; 117(13): 134902. Doi: 10.1063/1.4916805.

Zhu QL, Li XC, Jiang ZB, Wei N. Impacts of CO2 leakage into shallow formations on groundwater chemistry. Fuel Processing Technology 2015; 135: 162-7.

Qianlin Zhu obtained his M.Sc. and Ph.D. degrees in geotechnical engineering from Hohai University in 2008 and Institute of Rock and Soil Mechanics, Chinese Academy of Sciences in 2011, respectively. He is a research assistant at Institute of Low Carbon Energy, China University of Mining and Technology (Xuzhou). His research interests include CO2 geological storage, permeability theory in geotechnical engineering and isotopic hydrology. In recent years, he is striving to study two-phase displacement flow in porous media and coupled fluid flow-geomechanics characteristics in CO2 saline aquifer storage.

ACCEPTED MANUSCRIPT

The authors wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.