Pore-Scale Investigation of Carbon Dioxide-Enhanced Oil Recovery

Mar 29, 2017 - Carbon dioxide (CO2) enhanced oil recovery is a green and promising way to produce oil and reduce the rapid growth of carbon dioxide ...
0 downloads 0 Views 3MB Size
Subscriber access provided by University of Newcastle, Australia

Article

Pore-scale investigation of carbon dioxide enhanced oil recovery Guangpu Zhu, Jun Yao, Aifen Li, Hai Sun, and Lei Zhang Energy Fuels, Just Accepted Manuscript • DOI: 10.1021/acs.energyfuels.7b00058 • Publication Date (Web): 29 Mar 2017 Downloaded from http://pubs.acs.org on April 4, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Energy & Fuels is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Pore-scale Investigation of Carbon Dioxide Enhanced Oil Recovery Guangpu Zhu, Jun Yao*, Aifen Li, Hai Sun, Lei Zhang* School of Petroleum Engineering, China University of Petroleum, Qingdao 266580, China

Abstract Carbon dioxide (CO2) enhanced oil recovery is a green and promising way to produce oil and reduce the rapid growth of carbon dioxide released to atmosphere. Pore-scale understanding of CO2 displacement phenomena is important to enhance oil recovery in porous media. In this work, a direct numerical simulation method is employed to investigate the drainage process of CO2 in an oil-wet porous media. The interface between oil and CO2 is tracked by the phase field method. The capacity and accuracy of models are validated using a classic benchmark: a bubble rising process. A series of numerical experiments are performed over a large range of gravity number, capillary number and viscosity ratio to investigate the flooding process of CO2 in a porous media. The results show that the pressure in the main CO2 flow path decreases dramatically after CO2 breaks through outlet. Oil begins to re-flow into large pores which have been occupied by CO2 before. This phenomenon has an important impact on the final saturation distribution of CO2. Increasing viscous force is the dominant mechanism for improving oil recovery. Selecting an appropriate depth is the primary consideration to reach the maximum recovery before injecting CO2 into subsurface. Abnormal high pressure formation is a good choice for CO2 sequestration. Gravity fingers improve the sweep area of CO2 when the viscous force is small. The oil recovery increases with the increase in contact angle. It is difficult to reach the final steady state of saturation due to the “snap-off and supplement” dynamic balance in porous media when both injection velocity and contact angle are small.

1 Introduction Anthropogenic emission of CO2 is considered as the main contributor to the global greenhouse effect.1 One promising approach to reduce the rapid growth of CO2 released to the atmosphere2 is to sequestrate it into the oil3 and gas reservoirs or deep saline aquifers.1, 4 In 1972, the first commercial-scale sequestration of CO2 for enhancing oil recovery (EOR) was initiated in West Texas.5 Since then, numerous technical improvements and operating practices have been developed to CO2 EOR in petroleum industries. A research conducted by the Department of Energy (DOE) found that 137 billion barrels of additional technically recoverable domestic oil can be provided by the Next Generation CO2 EOR technologies.6 In petroleum industries, CO2 is injected into reservoirs where pressures are normally higher than 7.4 MPa and temperatures are greater than 304.15 K.7, 8 CO2 may exist in a supercritical state under these conditions and it displaces oil from the pore space as a non-wetting phase in a variety of saturation patterns.7, 9 Fingering phenomena1, 9-11 are normally observed during flooding processes, which limit the oil recovery of reservoirs. To improve the supercritical CO2 flooding efficiency, clear understandings of pore-saturation levels, drainage stability and fluid flow patterns are of paramount importance. However, such understandings are complicated by the fluid properties (e.g., density, viscosity and interfacial tension), hydrodynamic forces and the properties of porous media12 (e.g., wettability, pore size distribution). Pore-scale visualization experiments13-15 provide a powerful tool to investigate two phase displacement phenomena. The X-ray Computed Tomography (CT) based rock core12 and pore

ACS Paragon Plus Environment

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 14

network patterns10, 11, 16 are currently two main porous media used in pore-scale experiments. Natural rock core flooding experiments based on the CT suffer from difficulties in manipulating porosity, conductivity and wetting properties independently.1 These limitations can be easily overcome by two-dimensional (2D) pore network patterns, which can be etched into materials, such as glass, silicon and polydimethylsiloxane (PDMS).11 Pore network patterns allow for direct visualization using cameras with or without fluorescent microscopy, and subsequent quantitative evaluation may lead to better understandings of physical displacement processes. For example, Lenormand et al.16 conducted a series of displacement experiments and established a 2D phase diagram depicting regions for capillary fingering, viscous fingering and stable displacement. Zhang et al.10 complemented this work using water-wet micromodels over a large range of capillary number and viscous ratio. Nevertheless, their studies limited to capillary and viscous force dominated flows where gravity is neglected. Ewing and Berkowitz17 extended Lenormand et al.’s work and established a three-dimensional (3D) phase diagram considering the effect of gravity. However, this extension was limited to a qualitative description. Pore-scale numerical simulations can complement experimental studies, providing an efficient and economical way to explore the influence of flow parameters in various porous media.9, 18-23 Despite a large number of numerical pore-scale investigations, only a limited number of them of CO2 sequestration have been published to date. Most of these researches use the smoothed particle hydrodynamics (SPH) model7 and the Lattice-Boltzmann method1, 9, 23, 24 and only some of gravitational, viscous, or capillary forces are included in the simulations.7 Moreover, almost none of them analyzes the mechanical mechanisms for flow phenomena according to our current level of knowledge. Here a direct numerical simulation (DNS) method is used to investigate pore-scale flooding processes. The main advantage22 of DNS method is the natural and consistent way that flow models can be extended to model processes other than just the flow. It can simulate fluid flow in a complex geometry with heterogeneous pore space in conjunction with well-established numerical methods (e.g., finite volume method, finite element method) used in computational fluid dynamics (CFD) and high-performance computation.19, 21, 25 However, the physics of interface movement are not included in Navier-Stokes (N-S) equation. Consequently, for multiphase flow, the N-S equation should be coupled with an interface tracking method.26-28 The phase field method26, 28 is employed to track the interface between oil and supercritical CO2 in this study and details of this method are elaborated in section 2.2. In this study, a direct numerical simulation method is employed to simulate fluid flow in a porous media. The interface between oil and CO2 is tracked by the phase field method. A series of numerical experiments are performed to investigate flow behaviors of CO2 and oil recovery over a large range of gravity number, capillary number and viscosity ratio. Some factors (e.g., wettability) are also analyzed to explore their potential impacts on drainage processes.

2 Pore-scale Numerical simulation 2.1 Governing equations The N-S equation represents the conservation of momentum for fluid flow. For an incompressible fluid, it can be written as,18, 19, 27 T  ∂u  ρ  + ( u ⋅ ∇ ) u  = ∇ ⋅ − pI + µ ∇u + ( ∇u )  + ρ g + Fst (1)    ∂t  The conservation of mass equation for an incompressible flow is given by,18, 19, 27

(

)

ACS Paragon Plus Environment

Page 3 of 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

∇⋅u = 0

(2) Here u is the flow velocity of the fluid, m/s; ρ denotes the density of the fluid, kg/m ; I is the unit vector; t is the time, s; µ is the viscosity of fluid, Pa·s. g is the gravitational acceleration, m/s2; Fst is the term representing interfacial tension force, Pa/m. 3

2.2 Phase field method The phase field method treats the interface between two immiscible phases as a physically diffuse thin layer. The dynamics of the thin layer are affected by convections and diffusions. The phase field variable ϕ is introduced to distinguish phases and interface. For example, phase A assumes the value ϕA is -1 while phase B takes the value ϕB as -1 and interface region can be described by the transition from -1 to 1. According to the theory of Ginzburg-Landau, the free energy of interface region for an isothermal two-phase system can be expressed as, 29 1 2  F = ∫  f (φ ) + λ ∇ φ  d V V 2  

(3)

The mixing free energy of interface region reaches the minimum in the equilibrium state. Where, f (ϕ) is the bulk energy, which is a typical double-well positive function and it can be written as,28, 29 f (φ ) =

λ 4ε 2

(φ + 1) (1 − φ ) 2

2

(4)

λ is the magnitude of the mixing energy, N; ε is the interface thickness, m; 1/2λ|▽ϕ|2 is the gradient energy. The chemical potential of interface is the rate of change of F with respect to the phase field variable,18  φ (φ − 1)  δF  = λ  −∇ 2φ + δφ ε2    2

G=

(5)

Cahn and Hilliard generalized the theory to time dependent situations by approximating interfacial diffusion fluxes as being proportional to chemical potential gradients. The convective Cahn-Hilliard equation can be written as,18, 29 ∂φ  γλ   ∂t + u ⋅∇φ = ∇ ⋅  ε 2 ∇ψ     ψ = −∇ ⋅ ε 2∇φ + (φ 2 − 1)φ 

(6)

Where γ is the mobility, which denotes the moving velocity of the interface under a unit driven force, m3·s/kg. The surface tension acting on the fluid-fluid interface can be represented as,18, 29 (7) Fst = G∇φ Note that the density and viscosity of fluid in above equations are defined as functions of the phase field variable, as follows, ρ=

1+φ 1−φ ρn + ρw 2 2

(8)

µ=

1+φ 1−φ µn + µw 2 2

(9)

where the subscripts (w) and (n) denote the wetting phase and non-wetting phase, respectively.

ACS Paragon Plus Environment

Energy & Fuels

3 Results and Discussion 3.1 Validation of numerical model The N-S equation and phase field equation are highly coupled together and COMSOL MULTIPHYSICS19, 30 software is employed to get numerical solutions. A benchmark case is provided here to validate our numerical models. The diameter of the circular bubble (fluid 1) is D=0.5 m, as illustrated in Figure 1. The height of the rectangular domain is H=4 D and the width is W=2 D. The bubble is surrounded by fluid 2 and its center locates at [0.5 m, 0.5 m] at the initial stage. The densities of fluid 1 and fluid 2 are 100 kg/m3 and 1000 kg/m3, respectively. The viscosities of fluid 1 and fluid 2 are 1 Pa·s and 10 Pa·s, respectively. The acceleration is 0.98 m/s2 (body force). The interface tension is 24.5 N/m. The bubble will rise due to the buoyancy force. No-slip boundary conditions are set at the top and bottom boundaries, whereas free slip conditions are imposed on the vertical walls. The mass center and rising velocity of the bubble at different times are calculated respectively, as shown in Figure 2. The simulation results agree quite well with the results of Hysing et al.31, which validates the correctness of numerical models.

Figure 1 Initial configuration and boundary conditions of the bubble rising case 0.25

1.1 The result of S. Hysing et al Our result

0.20

rise velocity (m/s)

1.0

center of mass (m)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 14

0.9 0.8 0.7 0.6

0.15

The result of S. Hysing et al Our result

0.10

0.05

0.5 0.0

0.5

1.0

1.5

2.0

2.5

3.0

0.00 0.0

0.5

t (s)

1.0

1.5

2.0

2.5

3.0

t (s)

Figure 2 Time evolution of the mass center (left) and rising velocity (right) of the bubble 3.2 Numerical investigation of carbon dioxide EOR Figure 3 gives the schematic of the porous media. The dimensions of the porous media are 8100 µm × 4080 µm. It contains two different permeability zones. The top of the domain (dense zone) is created with closely packed layers of soil grains. The thickness of dense zone is 650 µm

ACS Paragon Plus Environment

Page 5 of 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

and the diameter of each grain in this region is 300 µm. This low permeability zone exists as a cap-rock, which can curtail the upward migration of CO2. The rest high permeability zone consists of a staggered periodic array of circular grains, 260 µm in diameter, with 220 µm pore bodies and 72 µm pore throats.

Figure 3 Schematic of the porous media. Some dimensionless numbers are defined to represent different flow conditions in reservoirs. The capillary number7 is defined as Ca=µcUin/σ, where µc is the viscosity of CO2, Pa·s; Uin is the average injection velocity at the inlet, m/s; σ is the interface tension coefficient, N/m. Ca represents the effect of viscous force versus capillary force. The gravity number7 represents the ratio of gravity force versus viscous force and is defined as Gr=∆ρgφδ2/µcUin, where ∆ρ is the density difference between two phases, kg/m3; φ is the porosity of porous media, dimensionless; δ is the size of pore throat, m. The viscosity ratio is defined as M=µc/µo, where µo is the viscosity of oil, Pa·s. Initially, the oil-wet porous media is saturated with oil. To observe the effect of gravity force clearly, CO2 is injected form the bottom part (0-1300 µm) of left hand side at a constant flow rate. A pressure of 0 Pa is set at the outlet. The rest boundaries are assumed to be solid walls. Detailed simulation parameters of each section are provided in a supplementary file. A series of numerical experiments are performed over a large range of gravity number, capillary number and viscosity ratio (M=0.001, M=0.01 and M=0.1) to investigate systematically flooding processes of supercritical CO2 in an oil-wet porous media. For each case, the simulation is run until the saturation of CO2 reaches the steady state. 3.2.1 Flow patterns and mechanism analysis Figure 4 demonstrates the final distribution of CO2 in porous media at M=0.001. The non-wetting CO2 has completely overcome the entry pressure and formed continuous flow paths across the whole porous media. Fingering phenomena can be clearly observed during the flooding process due to the unfavorable viscosity ratio M. These fingers are associated with capillary, gravity and viscous force under different flow conditions in reservoirs.

(a) logCa=-4.62

Gr=33.53

(b) logCa=-4.25

ACS Paragon Plus Environment

Gr=14.37

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(c) logCa=-4.09

Gr=10.06

Page 6 of 14

(d) logCa=-3.92

Gr=6.71

(e) logCa=-3.80 Gr=5.03 (f) logCa=-3.62 Gr=3.35 Figure 4 Final CO2 distributions in the porous media (M=0.001). The white circles are solid grains. The blue and red represent oil and CO2, respectively. At the low capillary number (logCa=-4.62) and high gravity number (Gr=33.53), the effect of viscous force is negligible and flow behaviors of liquid CO2 in porous media are mainly controlled by the capillary force and gravity force. However, a long loosely connected flow path is found in Figure 4 (a), which is the typical characteristic of viscous flow. Viscous finer should not appear under this situation. To investigated fundamental reasons for this “viscous finger”, it is necessary to analyze the whole flooding process.

(a) t=0.08 s (b) t=0.109 s (c) t=0.2 s Figure 5 CO2 distributions during the flooding process (lgCa=-4.62, Gr=33.53, M=0.001) Flow behaviors in Figure 5 show the typical characteristics of capillary fingering. The injected CO2 prefers to advance from a large pore in any directions due to the low entry pressure. The CO2 has almost occupied a pore body completely before it overcomes capillary pressure and flows into nearby pores. Nevertheless, the phenomenon that CO2 finger processes into new pores in the backward direction is not observed during displacement processes and it may be attributed to the homogeneous structure of porous media. Before the CO2 breaks through the outlet, only capillary fingers or gravity fingers are observed and no viscous finger appears. This results are consistent with our expectation at the low logCa and high Gr. Additionally, no CO2 plume break-up phenomenon (snap-off) is observed when it progresses in the porous media, which indicates that the capillary force is not large enough although flow behaviors are dominated by it. Once CO2 column flows away from the outlet, some changes in the distribution of CO2 will lead to the so called “viscous flow”.

ACS Paragon Plus Environment

Page 7 of 14

1200

1000

800

p/Pa

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

600

400

200 0.100

0.102

0.104

0.106

0.108

0.110

0.112

0.114

t/s

Figure 6 The pressure change at point R near the outlet Figure 6 shows the change in pressure at the point R near the outlet. Before the CO2 front reaches R, the fluid pressure in R has been steadily increasing. When the CO2 plume flows through R (0.107 s), the pressure in R sees a significant increase and the figure peak occurs at 0.109 s. We denote pressures in CO2 plume and oil phase as pCO2 and po, respectively. The capillary pressure is written as pc. When the CO2 displaces oil from the pore space (0.1-0.109 s), pCO2 is larger than the sum of po and pc. After the CO2 front breaks through the outlet (0.109-0.113 s), the pressure (pCO2) in main flow path decreases dramatically, as illustrated in Figure 6. At this time, pCO2 is significantly lower than the sum of po and pc and the imbibition of oil occurs. The oil begins to re-flow into large pores which have been occupied by CO2 before. The imbibition process leads to the shrinkage of main flow path (CO2), as illustrated in Figure 5 (c). Eventually, when pressure comes to balance (pCO2=po + pc), the flow reaches a steady state and the flow path shows the characteristics of viscous flow in Figure 4 (a). Based on the above analysis, the so called “viscous finger” in Figure 4 (a) is caused by the imbibition of oil and outlet effect is responsible for the imbibition process. Some vertical flow paths can be found in Figure 4 (a) and these paths are normally treated as gravity fingers because of a large Gr. Although gravity force plays an important role under this condition, the flow patterns of CO2 are also dominated by capillary force at a low capillary number (-4.62). The injected CO2 prefers to advance into a large pore in any directions and it may flow along vertical flow paths. Hence, we can not sure whether the two flow paths belong to gravity fingers. Details of this issue will be discussed in section 3.2.2. To sum up, at the low Ca and high Gr, flow behaviors of CO2 in Figure 4 (a) are determined by capillary force, gravity force and outlet effect together. Figure 4 (b) shows the same flow patterns with Figure 4 (a). With the increase in injection flow rate, the effect of viscous force strengthens. In Figure 4 (c) (logCa=-4.09, Gr=10.06), capillary fingers and “viscous fingers” are observed. Nevertheless, it is difficult to judge whether these “viscous fingers” are cause by viscous force or outlet effect. Additionally, obvious vertical flow paths occur in Figure 4 (c). These paths may be caused by gravity effect or some other reasons. To find out fundamental mechanisms for these phenomena, studying the drainage process is necessary.

ACS Paragon Plus Environment

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 14

(a) t=0.016 (b) t=0.017 (c) t=0.018 (d) t=0.019 Figure 7 Breakups of CO2 plume during the displacement process (logCa=-3.92 and Gr=6.71) Figure 7 demonstrates breakups of CO2 during the flooding process at logCa=-3.92 and Gr=6.71. Breakups of non-wetting fluid is normally caused by two different mechanisms: snap-off and viscous breakup. Snap-off is an interfacial instability and it occurs when capillary force dominates flow behaviors. As mentioned before, snap-off phenomenon is not observed in a capillary force dominated case (Figure 4(a), logCa=-4.62). Snap-off phenomenon is impossible for the case with relatively weak effect of capillary force (logCa=-3.92). Therefore, the breakups in Figure 7 belong to typical viscous breakup. Viscous breakups indicate that viscous force has played a significant role in the displacement process.

(a) t=0.051 s

(b) t=0.061 s

(c) t=0.13 s (d) t=0.16 s Figure 8 The formation process of vertical flow paths Figure 8 demonstrates the formation process of vertical flow paths clearly. We can see that these vertical flow paths are not gravity fingers, and they are formed by downward flow. The effect of gravity force is weak in this situation (Gr=6.71). The main reason for the formation of vertical flow paths is the pressure difference between different flow plumes. After the main flow plume breaks through the outlet, the pressure in the lower main plume reduces significantly. The pressure difference between upper and lower flow paths is large enough to overcome capillary force and buoyancy force. Besides, some typical characteristics of capillary fingering are observed near the inlet. Based on the above analysis, it can be concluded that the capillary force, the viscous force and outlet effect dominate the flow behaviors in Figure 7 (c), and the effect of gravity force is not obvious. The effect of viscous force increases with the increase in flow rate at the inlet. In the situation of Figure 7 (f), the dominant viscous force causes some fingers to break up. Outlet effect may have some impacts on final CO2 distributions and other effects can almost be neglected (e.g., gravity force). As the capillary number increases, the number of CO2 flow paths connected to the outlet increases, which is consistent with previous experimental and numerical results.9, 11

ACS Paragon Plus Environment

Page 9 of 14

90 ρC=150 kg/m3, M=0.001 ρC=750 kg/m3, M=0.001

80

ρC=150 kg/m3, M=0.01 ρC=750 kg/m3, M=0.01

70

ρC=150 kg/m3, M=0.1

SC/%

60

3

ρC=750 kg/m , M=0.1

50 40 30 20 10 -4.6

-4.4

-4.2

-4.0

-3.8

-3.6

-3.4

-3.2

-3.0

-2.8

logCa

Figure 9 Final saturations of CO2 in the porous media in different conditions To investigate EOR with CO2 sequestration, a series of scenarios are performed to simulate different underground conditions. Figure 9 gives the final CO2 saturations, which reflect oil recoveries in porous media. As the capillary number logCa increases, oil recovery sees a significant increase. According to the definition of Ca=µcUin/σ, we can conclude that increases in viscosity and injection velocity (inlet flow rate) of CO2 are keys to improve oil recovery. For the same viscosity ratio and capillary number, the flooding efficiency of weight CO2 (corresponds to small Gr) is higher than that of light CO2 (large Gr). In fact, both large logCa and small Gr indicate a strong effect of viscous force. Hence, increasing viscous force is the dominant mechanism for improving oil recovery. In field practices, the properties of CO2 are strongly affected by reservoir conditions. As the formation temperature decreases or pressure increases, both viscosity and density of CO2 increase, which are conducive to enhance oil recovery. Low-temperature or high-pressure formation is more suitable for CO2 sequestration. However, with the increase in depth, both temperature and pressure increase. It is difficult to find a formation with low temperature but high pressure. Only some abnormal high pressure formations can meet this requirement. Therefore, selecting an appropriate depth is the primary consideration to reach maximum recovery before injecting CO2 into subsurface. M=0.1 M=0.01 M=0.001

0.03

average inlet velocity Uin (m/s)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Uin=0.03 m/s, M=0.1

0.05 Uin=0.03 m/s, M=0.01

0.07 Uin=0.03 m/s, M=0.001

0.1

0

5

10

15

20

25

30

35

40

45

50

55

SC (%)

Figure 10 Saturations of CO2 at different viscosity ratios and average inlet velocity

ACS Paragon Plus Environment

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 14

A favorable viscosity ratio M means a high oil recovery according to our previous analysis. Nevertheless, this conclusion maybe not true when the injection flow rate is small. Figure 10 gives the final CO2 saturations with various average inlet velocities at different viscosity ratios. An interesting phenomenon can be observed that SC (M=0.01)> SC (M=0.001)>SC (M=0.1) when the average inlet velocity Uin is 0.03 m/s. The result that SC(M=0.01)> SC(M=0.001) is consistent with the previous analysis. However, the interesting phenomenon that SC(M=0.001)> SC(M=0.1) is difficult to understand. Comparing CO2 distributions at different viscosity ratios in Figure 10, we can find that the pressure drop is not large enough to form a continuous flow across the whole zone even it is able to overcome the entry pressure at an average inlet velocity Uin of 0.03 m/s, and this may be the main reason for SC (M=0.001)> SC (M=0.1). The conclusion that a favorable viscosity ratio M means a large oil recovery can be true when the CO2 breaks through the outlet. 3.2.2 The effect of gravity force To investigate the effect of gravity force, six different cases at the same viscosity ratio (M=0.001) are investigated. The results are shown in Figure 11.

(a) logCa=-4.62 Gr=33.53

(b) logCa=-4.25 Gr=14.27

(c) logCa=-3.80 Gr=5.03

(d) logCa=-4.62 (e) logCa=-4.25 (f) logCa=-3.80 Figure 11 Final CO2 distributions in porous media at a viscosity ratio M of 0.001. (a)-(c): gravity is considered in simulations; (d)-(f): no gravity is considered. Figure 11 (a)-(c) show CO2 flow behaviors as gravity is incorporated into simulations. Figure 11 (d)-(f) give the results without the effect of gravity force. Comparing the vertical flow paths in Figure 11 (a) and Figure 11 (d), it can be found that the first vertical path in Figure 11 (a) is a typical capillary finger and the second finger is mainly controlled by gravity force at the high Gr and low logCa. With the increase in flow rate, the effect of capillary and gravity force decrease. Only a small gravity finger occurs in Figure 11 (b). No gravity finger shows in Figure 11 (c) and gravity has a limited effect on CO2 flow behaviors.

ACS Paragon Plus Environment

Page 11 of 14

55 50 45 40

SC (%)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

35 30

No gravity M=0.001 M=0.001 No gravity M=0.1 M=0.1

25 20 15 0.04

0.06

0.08

0.10

0.12

0.14

0.16

0.18

0.20

Uin (m/s)

Figure 12 Saturations of CO2 with different injection velocities at various viscosity ratios Figure 12 gives saturations at different viscosity ratios. For comparison purpose, the x-coordinate is set as injection velocity. When the injection velocity is low (small logCa and large Gr), flow behaviors of CO2 are dominated by capillary force and gravity force. Gravity fingers improve the sweep area of CO2 and gravity force has a positive impact on oil recovery. As the viscous force increases, the case with no gravity effect has better flooding efficiency than that with gravity effect at the same viscosity ratio. Gravity effect plays a resistance role. 3.2.3 The wettability of porous media To investigate the impact of wettability on oil recovery, six different cases at the same viscosity ratio (M=0.01) are performed to simulate various injection and formation conditions.

(a) θ=75°

(b) θ=45°

(c) θ=30°

(d) θ=75° (e) θ=45° (f) θ=30° Figure 13 Final distributions of CO2 in different wetting conditions at M=0.01. (a)-(c): logCa=-3.40, Gr=2.01; (d)-(f): logCa=-3.85, Gr=5.75. Flow behaviors at logCa=-3.40 and logCa=-3.85 are similar, as illustrated in Figure 13. The flow patterns at logCa=-3.85 (Uin=0.07 m/s) are taken as examples to analyze the impact of wettability. When the contact angle θ is 75°, it is easy for non-wetting CO2 to enter and occupy pore bodies completely due to the low entry pressure (pc=2σcosθ/r), which contributes to the large sweep area of CO2, as shown in Figure 13 (d). Capillary fingers can be observed during flooding processes. Also, the outlet effect has a limited impact on the final distribution of CO2 because of the low capillary force and a small change in pressure of main flow path after the outlet is broken through. The capillary pressure increases as the contact angle decreases. A part of the domain is selected to elaborate the flow patterns in detail in Figure 13 (e).

ACS Paragon Plus Environment

Energy & Fuels

Figure 14 Flow behaviors of CO2 in a selected domain (θ=45°) Figure 14 gives the flow behaviors in detail at contact angle of 45°. It can be clearly observed that CO2 must occupy a pore body completely before it processes into neighboring pores, which exhibits the displacement characteristics associated with capillary fingerings. When the CO2 front advances toward the outlet, the pressure in the main flow path decreases, as shown in Figure 15, and oil reflows into pore bodies. Hence, a part of the loosely connected flow path shows the characteristics of viscous fingering but it is not real viscous fingers. The drained CO2 flows along the main path and services as an important supplement of front. The outlet effect begins to play an important role in determining final flow behaviors of CO2 when the CO2 front breaks through the outlet. 2400 2200 2000

p/Pa

1800 1600 1400 1200 1000

CO2 breaks through the outlet 800 0.055

0.060

0.065

0.070

0.075

0.080

t (s)

Figure 15 The pressure in point M before CO2 break through the outlet (θ=45°) 60 55 Uin=0.2 m/s, M=0.01 50

Uin=0.07 m/s, M=0.01

45

SC (%)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 14

40 35 30 25 20 15 20

30

40

50

60

70

contact angle (°)

Figure 16 The oil recovery at different contact angles The flow behaviors of CO2 at contact angle of 30° are similar to that at 45° and they are

ACS Paragon Plus Environment

Page 13 of 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

determined by the capillary force and outlet effect. The capillary force at contact angle of 30° is almost 2.4 times than that at contact angle of 75°. Due to the large capillary pressure and low viscous force (low injection velocity), snap-off occurs32 during the flooding, as shown in Figure 13 (f). When the CO2 front breaks up from the main flow path and flows away from the outlet, injection fluid serves as an important supplement for the main flow path and then snap-off phenomenon occurs again. It is difficult to reach the final steady state due to the “snap-off and supplement” dynamic balance in porous media. The oil recovery decreases with the decrease in contact angle, as shown in Figure 16, which is consistent with the above analysis.

Conclusion In this study, a direct numerical method together with the phase-field method is employed to investigate the flow behavior of CO2 and oil recovery during the drainage process in an oil-wet porous media. The specific conclusions of this study are summarized below: (1) The pressure in the main CO2 flow path decreases dramatically after CO2 breaks through the outlet. Oil begins to re-flows into the large pores which have been occupied by CO2 before. This phenomenon has an important impact on final saturation distribution of CO2. Increasing viscous force is the dominant mechanism for improving oil recovery. (2) Selecting an appropriate depth is the primary consideration to reach maximum recovery before injecting CO2 into the subsurface. Abnormal high pressure formation is a good choice for CO2 sequestration. (3) Favorable viscosity ratio M means large oil recovery can be true when CO2 breaks through outlet. (4) Gravity fingers improve the sweep area of CO2 when viscous force is small. The oil recovery decreases with the decrease in contact angle. It is difficult to reach final steady state of saturation due to the “snap-off and supplement” dynamic balance in the porous media when the injection velocity and contact angle are small. Supporting Information Information about simulation parameters Author Information Corresponding Author *E-mail: [email protected] (J.Y.) *E-mail: [email protected] (L.Z.) Author Contributions *J.Y. and L.Z. contributed equally to this work. Notes The authors declare no competing financial interest.

Acknowledgement This work was financially supported by the National Science and Technology Major Project (2016ZX05011-001, 2016ZX05061), Shandong Provincial Natural Science Foundation, China (ZR2016EL09), National Natural Science Foundation of China (51490654, 51504276 and 51304232) and National Science Foundation for Young Scientists of China (51404291).

ACS Paragon Plus Environment

Energy & Fuels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Reference (1) Liu, H.; Zhang, Y.; Valocchi, A. J. Phys. Fluids 2015, 27, 052103. (2) Aycaguer, A.-C.; Lev-On, M.; Winer, A. M. Energy Fuels 2001, 15, 303-308. (3) Rudyk, S.; Spirov, P. Energy Fuels 2013, 27, 5996-6001. (4) Cao, S. C.; Dai, S.; Jung, J. Int. J. Greenhouse Gas Control 2016, 44, 104-114. (5) Parker, M. E.; Meyer, J. P.; Meadows, S. R. Energy Procedia 2009, 1, 3141-3148. (6) Li, L.; Yao, J.; Li, Y.; Wu, M.; Zhang, L. J. Nat. Gas Sci. Eng. 2016, 33, 30-36. (7) Bandara, U. C.; Tartakovsky, A. M.; Palmer, B. J. Int. J. Greenhouse Gas Control 2011, 5, 1566-1577. (8) Jung, J.-W.; Wan, J. Energy Fuels 2012, 26, 6053-6059. (9) Liu, H.; Valocchi, A. J.; Werth, C.; Kang, Q.; Oostrom, M. Adv. Water Res. 2014, 73, 144-158. (10) Zhang, C.; Oostrom, M.; Wietsma, T. W.; Grate, J. W.; Warner, M. G. Energy Fuels 2011, 25, 3493-3505. (11) Zhang, C.; Oostrom, M.; Grate, J. W.; Wietsma, T. W.; Warner, M. G. Environ. Sci. Technol 2011, 45, 7581-7588. (12) Chaudhary, K.; Bayani Cardenas, M.; Wolfe, W. W.; Maisano, J. A.; Ketcham, R. A.; Bennett, P. C. Geophys. Res. Lett. 2013, 40, 3878-3882. (13) Kim, Y.; Wan, J.; Kneafsey, T. J.; Tokunaga, T. K. Environ. Sci. Technol 2012, 46, 4228-4235. (14) Li, R.; Jiang, P.-X.; Gao, C.; Huang, F.; Xu, R.; Chen, X. Energy Fuels 2016. (15) Ogden, S.; Bodén, R.; Do-Quang, M.; Wu, Z.; Amberg, G.; Hjort, K. Microfluid. Nanofluid. 2014, 17, 1105-1112. (16) Lenormand, R.; Touboul, E.; Zarcone, C. J. Fluid Mech. 1988, 189, 165-187. (17) Ewing, R. P.; Berkowitz, B. Water Resour. Res. 1998, 34, 611-622. (18) Zhu, G.; Yao, J.; Zhang, L.; Sun, H.; Li, A.; Shams, B. Langmuir 2016, 32, 11736-11744. (19) Gunde, A. C.; Bera, B.; Mitra, S. K. Energy 2010, 35, 5209-5216. (20) Chen, L.; Luan, H.-B.; He, Y.-L.; Tao, W.-Q. Int. J. Therm. Sci. 2012, 51, 132-144. (21) Qiao, Z.; Sun, S. SIAM. J. Sci. Compyt. 2014, 36, B708-B728. (22) Peszynska, M.; Trykozko, A. Computat. Geosci. 2013, 17, 623-645. (23) Zhang, L.; Kang, Q.; Yao, J.; Gao, Y.; Sun, Z.; Liu, H.; Valocchi, A. J. Sci. China: Technol. Sci. 2015, 58, 1375-1384. (24) Fei, K.; Hong, C. Microfluid. Nanofluid. 2007, 3, 77-88. (25) Yue, P.; Feng, J. J.; Liu, C.; Shen, J. J. Fluid Mech. 2004, 515, 293-317. (26) Boyer, F.; Lapuerta, C.; Minjeaud, S.; Piar, B.; Quintard, M. Transp. Porous Media 2010, 82, 463-483. (27) Tryggvason, G.; Bunner, B.; Esmaeeli, A.; Juric, D.; Al-Rawahi, N.; Tauber, W.; Han, J.; Nas, S.; Jan, Y.-J. J. Comput. Phys. 2001, 169, 708-759. (28) Yue, P.; Zhou, C.; Feng, J. J.; Ollivier-Gooch, C. F.; Hu, H. H. J. Comput. Phys. 2006, 219, 47-67. (29) Villanueva, W.; Amberg, G. Int. J. Multiphase Flow 2006, 32, 1072-1086. (30) Zhu, G.-p.; Yao, J.; Sun, H.; Zhang, M.; Xie, M.-j.; Sun, Z.-x.; Lu, T. J. Nat. Gas Sci. Eng. 2016, 28, 305-316. (31) Hysing, S.-R.; Turek, S.; Kuzmin, D.; Parolini, N.; Burman, E.; Ganesan, S.; Tobiska, L. Int. J. Numer. Methods Fluids 2009, 60, 1259-1288. (32) Deng, W.; Cardenas, M. B.; Bennett, P. C. Adv. Water Res. 2014, 64, 34-46.

ACS Paragon Plus Environment

Page 14 of 14