Computer Modeling of the Displacement Behavior of Carbon Dioxide in Undersaturated Oil Reservoirs

The injection of CO2 into oil reservoirs is performed not only to improve oil recovery but also to store CO2 captured from fuel combustion. The objective of this work is to develop a numerical simulator to predict quantitatively supercritical CO2 flooding behaviors for Enhanced Oil Recovery (EOR). A non-isothermal compositional flow mathematical model is developed. The phase transition diagram is designed according to the Minimum Miscibility Pressure (MMP) and CO2 maximum solubility in oil phase. The convection and diffusion of CO2 mixtures in multiphase fluids in reservoirs, mass transfer between CO2 and crude and phase partitioning are considered. The governing equations are discretized by applying a fully implicit finite difference technique. Newton-Raphson iterative technique was used to solve the nonlinear equation systems and a simulator was developed. The performances of CO2 immiscible and miscible flooding in oil reservoirs are predicted by the new simulator. The distribution of pressure and temperature, phase saturations, mole fraction of each component in each phase, formation damage caused by asphaltene precipitation and the improved oil recovery are predicted by the simulator. Experimental data validate the developed simulator by comparison with simulation results. The applications of the simulator in prediction of CO2 flooding in oil reservoirs indicate that the simulator is robust for predicting CO2 flooding performance.


INTRODUCTION
During the period since 1750, CO 2 concentration in the atmosphere has increased from 280 ppm to 368 ppm in 2000 and 388 ppm in 2010 (Rackley, 2010).Annual CO 2 emission in the 11 southeastern states (Alabama, Arkansas, Florida, Georgia, Louisiana, Mississippi, North Carolina, South Carolina, Tennessee, Virginia and east Texas), of USA is up to about 1045 million metric tons (Petrusak et al., 2009).
Chinese existing and planned projects (nearly 400) for making ammonia, methanol and other fuels from coal will emit some 270 million tonnes of CO 2 per annum (Zheng et al., 2010).
The CO 2 concentration in the atmosphere is increasing at the rate of 2 ppm yearly.Increase in atmospheric emission of CO 2 with the attendant global warming and environmental degradation are driven by global energy demand (Lal, 2008).
Depleted or mature oil fields may provide feasible sites for storing CO 2 in porous and permeable reservoirs (Petrusak et al., 2009;Shtepani, 2006).CO 2 is injected into an oil reservoir to reduce oil viscosity and interfacial tension, and cause oil swelling, which improves oil recovery (Jaramillo et al., 2009).In most oil reservoirs, CO 2 is a supercritical fluid as reservoir conditions are above its critical temperature (31.4°C) and pressure (7.4 MPa).Supercritical CO 2 is a favorable flooding agent for EOR.Many oil fields in main oil production countries offer good opportunities for CO 2 injection in oil formations for EOR purposes.The history of CO 2 flooding for EOR purpose in oil fields and successful projects verifies that it can improve oil recovery to a larger extent.However, one adverse factor of CO 2 flooding for EOR is the problem induced by asphaltene precipitation and deposition.It may not only lead to the formation damages by reducing porosity and permeability (Shedid and Zekri, 2006) but also have some adverse influences on production facilities such as tubing and pumps (Thawer et al., 1990;Rashid and Sultan, 2003).
During CO 2 flow in oil reservoirs, the main physical and chemical interactions among the water, crude oil and rocks include: -the immiscible and miscible displacement of oil, water and CO 2 ; -convection and diffusion; -phase change and behavior of in-situ fluids on flow and mass transfer; -asphaltene precipitation and its effects on porosity and permeability.
TOUGH2-ECO2N (one member of TOUGH2 family) provides a flexible grid treatment and robust numerical solver to solve the multiphase flow equations in porous media.PSU-COALCOMP (Sams et al., 2002) and SIMED (Pinczewski and Stevenson, 2007) were reported to simulate the flow of CO 2 in coal beds.However, the above simulators provide the approaches for simulating CO 2 sequestration in aquifer formation or coal beds.TOUGH2-TMGAS (Battistelli and Marcolini, 2009) can be used to model the two-phase equilibrium between NaCl brine and non-aqueous phase mixtures made up by both organic and inorganic components.The main limitation of the approach is that only two-phase conditions are allowed.Especially, the physical and chemical mechanisms of CO 2 flooding for EOR differ from the CO 2 sequestration in aquifer and coal beds.More recently, ECLIPSE 300 and CMG-GEM modules has been developed for predicting the CO 2 flooding performances for EOR.Peng-Robinson's (PR) and Soave-Redlich-Kwong's (SRK) equations are selected to predict phase equilibrium and properties in the above simulators.However, high CO 2 concentration dissolved in crude oil at reservoir conditions will lead to the severe inaccurate prediction of CO 2 -oil mixture properties such as densities and compressibility factors.
The authors developed a non-isothermal compositional mathematical model as well as a reservoir simulator for quantitative description of processes in CO 2 flooding and subsurface storage in reservoirs.The formation damage caused by asphaltene precipitation and its effect on the performances of oil production are also considered in the model.The MMP and CO 2 solubility in crude oil were chosen as the criterions for phase transfer.Four simulation examples were designed to demonstrate the various functions of the simulator developed in this work.It is found that the predicted recovery flooded by CO 2 injection at a miscible condition has a good match with the experimental data.

MATHEMATICAL MODEL AND GOVERNING EQUATIONS 1.Assumptions and Considerations
Physical and chemical processes of CO 2 flow and transport are simulated by using a conceptual model of multiphase fluid flow and multi-component transport in porous media.
The CO 2 -EOR mathematical model is based on the following assumptions and considerations: -heat transfer in porous media is considered; -flow and transport mechanisms include multiphase Darcy's flow and multi-component diffusion, subject to adsorption and precipitation; -incorporation of a number of mass components: CO 2 , H 2 O, oil with some pseudo-oil components; -description of CO 2 adsorption on pore walls of reservoir formations; -dissolution of CO 2 in oil; -phase change and behavior of multiphase and multicomponent systems of CO 2 , oil and water on flow, displacement in reservoirs; -porosity and permeability changes induced by asphaltene deposition on pore walls are considered.

Phase Behavior and Fluid Types
The base case scenarios of fluid phases in the model are: -two phases of water and crude oil (pressure is lower than MMP without free CO 2 phase); -three phases of water, gas and oil under reservoir pressure below MMP; -two phases of water and CO 2 -oil mixture above MMP.

Phase Transfer and Condition
According to thermodynamic conditions and relative abundance of the different components, the fluids may exist in three scenarios, as shown in Figure 1.Arrows denote routes for appearance or disappearance of phases that are checked after the updates of primary variables.There are two key parameters (MMP and CO 2 equilibrium solubility in oil) determining phase switch in Figure 1.
Before CO 2 injection into undersaturated reservoirs, only oil and water phases coexist (see scenario 1).If the reservoir pressure is less than MMP or equal to MMP, scenario (1) will switch to scenario (2) if the CO 2 solubility exceeds its saturation state.If the pressure becomes larger than MMP, no free CO 2 releases and CO 2 -hydrocarbons will arrive at miscible state.In this case, scenario (1) will switch to scenario (3).
For scenario (2), when the pressure becomes lower than MMP or equal to MMP and if free gas disappears, then scenario (2) will switch to scenario (1).If the pressure exceeds MMP, the gas dissolves into oil.In this case, scenario (2) will switch to scenario (3).
For scenario (3), when the pressure becomes lower than MMP and if the CO 2 solubility does not reach its saturation state, then scenario (3) will switch to scenario (1).If the pressure becomes less than MMP, the free gas releases from oil.In this case, scenario (3) will switch to scenario (2).

MMP Calculation
A feasible approach for predicting MMP is to combine experimental data and consider sensitive hydrocarbon components affecting MMP.Glaso (1985) presented a MMP correlation, which not only requires less CPU time but considers the sensitive factors such as reservoir temperature, the heavy components (C 7+ ) and intermediate components (C 2 -C 6 ).If the international system of units is used in Glaso's correlation, then: see Equation ( 1 Phase scenarios and transfer path diagram (Note: gas phase is a special component mixture, including CO 2 , methane, ethane and other light hydrocarbons.At reservoir conditions, the phase may be at a supercritical state.Therefore, gas phase specified in this work might not be a real gas state condition).
O: Oil phase or Oil-CO 2 mixture phase; W: Water phase; G: Gas phase; P: Pressure; MMP: Minimum miscibility pressure; Sg: Gas (CO 2 ) saturation.X CO 2 : CO 2 mole fraction dissolved in oil phase; X CO 2 eq : Maximum CO 2 mole fraction dissolved in oil phase when equilibrium reaches at given pressure and temperature.

CO 2 Equilibrium Solubility in Oil
CO 2 equilibrium solubility in a crude oil depends on oil components, pressure and temperature.It is difficult to accurately calculate the solubility by using a correlation at any pressure and temperature.An experimental approach can provide CO 2 equilibrium solubility in a specific crude oil.In this work, the experimental solubility in the JL oil sample (CNPC, China) is shown in Figure 2. Interpolation approach is used to calculate the solubility in the oil at a specified pressure and temperature (P, T).

Governing Equations
The physical processes associated with CO 2 , oil and water flow in porous media are governed by the laws such as conservation of mass, momentum and energy.These physical laws can be represented mathematically at the macroscopic level by a set of partial differential or integral equations.To derive a set of generalized governing equations for multiphase fluid flow, multicomponent transport and heat transfer, we assume that these processes can be described using a continuum approach within a Representative Elementary Volume (REV) in a porous medium.In addition, a condition of local thermodynamic equilibrium and partitioning is assumed so that temperatures, phase pressures, densities, viscosities, enthalpies, internal energies and component concentrations (or mole fractions) at any time are the same locally, respectively, at each REV of the porous medium.

Mass Conservation Equations
According to mass conservation principles, a conservation equation of mass components in a porous medium can be written as follows: (2) where superscript k is the index for the components, k = 1, 2, 3, …, N c , with N c being the total number of mass components; A k is the accumulation term of component k; q k is an external source/sink term for mass component k; and F k is the "flow" term of mass or diffusive and dispersive mass transport.
The accumulation terms in Equation ( 2) for component k is written as: (3) where subscript β is an index for fluid phase (β = g for gas phase, = w for aqueous phase, = o for oil phase); φ is the porosity of porous media; ρ β is the density of phase β; and S β is the saturation of phase β; X k β is the mole fraction of component k in fluid β; δ k is the volume of component k deposited on the pore surfaces per unit bulk volume.Its adsorption rate may be described as: (4) where α k d is a capture rate coefficient of component k at pore walls.v o is the flow velocity of oil phase.R k s is the adsorption rate of component k on rock solids.Some heavy components such as asphalt and asphaltene have more opportunity to attach on pore walls, which leads to formation damage.
The mass component transport is governed by processes of advection, diffusion and dispersion.Based on T2R3D code of the TOUGH2 family, the total mass flow term (Scheidegger, 1961;Wu and Pruess, 2000) for a component k is written as: (5) where v β is a vector of volumetric flow rate for fluid flow, defined by Darcy's law to describe the multiphase flow as: where, P β , μ β and g are pressure, viscosity of fluid phase β and gravitational constant, respectively; z is the vertical coordinate.k ' is absolute permeability; and k ' rβ is the relative permeability to phase β.In Equation ( 5), D k β is the hydrodynamic dispersion tensor accounting for both molecular diffusion and mechanical dispersion for component k in phase β, defined by an extended dispersion model (Scheidegger, 1961) to include multiphase effects (Wu and Pruess, 2000) as: (7) where α β T and α β L are transverse and longitudinal dispersivities, respectively, in fluid β; τ is the tortuosity of a porous medium; d k β is the molecular diffusion coefficient of component k within fluid β; and δ ij is the Kronecker delta function (δ ij = 1 for i = j; and δ ij = 0 for i ≠ j), with i and j being coordinated indices.

Energy Conservation Equations
Heat transfer in porous media is in general a result of both convective and conductive heat transfer processes, while radiation may be ignored in most cases.These processes are complicated by mass transfer between multiphase fluids, multicomponents, associated changes in phases, internal energy and enthalpy.The energy conservation equations in a porous medium system (rock) may be described as: where ρ s is the density of rock grains; U β and U s are the internal energies of fluid β and rock grains, respectively.q T is the external source/sink term for energy; h β and h k β are the specific enthalpies of fluid phase β and of component k in fluid β, respectively; K T is the overall thermal conductivity; and T is temperature.

Constitutive Relationships and Fluid Properties
To complete the mathematical description of multiphase flow, multicomponent transport and heat transfer in porous media, the mass and energy-balance equations need to be supplemented with a number of constitutive equations.These correlations express interrelationships and constraints of physical processes, variables and parameters and allow the evaluation of secondary variables and parameters as functions of a set of primary unknowns or variables selected to make the governing equations solvable.The constitutive correlations are expressed as the following.
Constraint on the summation of total fluid saturation in porous media satisfies as: Constraint on mole fractions within phase β is written as: (10) In a multiphase system, the capillary pressure between two phases is defined as functions of fluid saturation: The relative permeability of a fluid phase in a two-phase system is normally assumed to be functions of its saturation: Although the relative permeabilities and capillary pressures functions can be described by previous correlations (Pruess et al., 1999), experimental work generally provides experimental data.In this work, relative permeabilities and capillary pressures are show in Figure 3.
The density of a fluid phase is treated as a function of pressure and temperature, as well as mole fraction (k = 1, 2, 3, ..., N c ): (13) The density of each phase can be regarded as a function of pressure, temperature and compositions.For the aqueous phase, the effects of hydrocarbons dissolved in the phase on its density can be neglected.A modification of the model (Garcia, 2001;Pruess, 2005) is: For gas phase, a version of the SRK equation with modified coefficients is used to gas density.The multicomponent versions of two parameters in SRK are written in terms of parameters a k , b k and gas phase mole fractions X g k of the individual components using the mixing rules recommended by Reid et al. (1987).
The density of CO 2 -oil mixture calculated by PR or SRK equations is inaccurate.Based on PR equation, we developed a New State Equation (NSE) to calculate the density: (15) where: = ( ) Table 1 gives the densities of CO 2 -oil mixture predicted by the equations (PR, SRK and NSE).The Comparison of experimental data and calculated results validates NSE.
The empirical expression of viscosity of a fluid is treated as a function of pressure, temperature and composition: (20) The mixture viscosity of each phase is obtained in terms of the pure component viscosities as: The viscosity of each component is calculated as a function of temperature from an equation due to Yaws et al. (1976): (22) Equilibrium partitioning X α h and X β h are the mole fraction of component k in phase α and β, respectively; and K h α:β is the equilibrium partitioning coefficient of component k between phases α and β.Their relation is: The relative permeabilities and capillary pressures.a) The relative permeability curve of oil and water; b) the relative permeability curve of gas and oil; c) the capillary pressure curve of oil and water; d) the capillary pressure curve of gas and oil.Partitioning coefficient depends on chemical properties of the component and is a function of temperature, pressure and composition: (24) Considering the choice of primary variables, we treat oil as reference phase in this work.The details of calculating equilibrium partitioning coefficient are discussed by Ahmed (2007).
Specific enthalpy of a liquid is expressed as: (25) Internal energy, U β , of liquid phase β is a function of composition, pressure and temperature.Specific enthalpy of a gas is expressed as: g is the specific internal energy of component k in the gas phase; C k g is the partial density of component k in gas phase.
The thermal conductivity of a porous medium is treated as a function of fluid saturation: (27) Porosity, φ 0 , is the effective porosity at a reference pressure, P 0 and a reference temperature, T 0 ; and C r and C T are the compressibility and thermal expansion coefficient of the medium, respectively: (28) The adsorption of heavy components on pore surfaces may lead to the reduction in a local porosity and permeability of the medium.The instantaneous porosity is expressed by: (29) where Δφ denotes the variation of the porosity caused by an adsorption of heavy components on pore walls and it is expressed by: (30) δ k can be calculated by the integration of Equation (4).That is: (31) To describe mathematically the formation damage caused by asphaltene precipitation, the heaviest pseudo-oil component (k = asphaltene) is approximately treated as asphaltene.δ k is generally equal to zero at initial time.The theory of asphaltene precipitation and deposition and its mathematical modeling was discussed in details by Ju et al. (2010).The change in permeability addressed in this study is caused by asphaltene.According to the reference (Liu and Faruk, 1993), the expression for instantaneous permeability changed by asphaltene precipitation is given by: (32) in which k ' 0 and φ 0 are the initial permeability and porosity; k ' and φ are instantaneous local permeability and porosity of a porous medium; and λ f is a constant for fluid seepage allowed by the plugged pores.The range of index n is from 2.5 to 3.5.f is a flow efficiency factor (Ju et al., 2002).

NUMERICAL SOLUTION SCHEME
The total number of mass components and energy equation is N c + 1.There is a total of N c mass conservation equations to solve for an isothermal reservoir; and N c + 1 equations to solve for a non-isothermal system.The selections of primary variables are listed in Table 2.
In Table 2, phase g denotes CO 2 phase of gaseous, supercritical, or liquid state; m is phase index for oil-CO 2 mixture; and X h m is a mole fraction of component k in oil-CO 2 mixture phase.
In this work, the integral finite-difference method is used for spatial discretization.Time discretization is carried out using a backward, first-order, fully implicit finite-difference scheme.The details of the discretization and numerical solution procedure are discussed in the TOUGH2 literature (Pruess et al., 1999).

NUMERICAL RESULTS AND DISCUSSIONS
This section demonstrates the functions of the developed simulator and validates the mathematical model describing the multiphase and multicomponent flow involving CO 2 displacement in oil formations.Although the model developed in this work can predict water solubility in oil phase, we neglect this effect because of very low solubility of water in the oil phase.In addition, CO 2 dissolving in the aqueous phase is also neglected for its slight effects on water properties.

The Comparison between Simulation and Slim Tube Data
In this section, the isothermal simulation result of one-dimensional miscible flooding is compared with the experimental data.The main parameters used in the numerical simulation are shown in The model grids and locations of the CO 2 injection and production sink are demonstrated in Figure 4.
Figure 5 shows that the relation between the oil recovery and injection PV of CO 2 is almost linear before CO 2 breakthrough (0.88 PV at horizontal axis) at the outlet of the model.After CO 2 breakthrough, the improvement in oil recovery slows down with continuous injection.It also shows that the oil recovery obtained by numerical simulation has a good match with experimental data.Only small errors (maximum relative error is less than 2%) between the predicted results and experimental data.The oil recovery after 1.2 PV injection of CO 2 is also very closed to the result of Dong et al. (2001).The comparison with our experimental data and the published work (Dong et al., 2001) validates the reliability of the mathematical model in case of a miscible displacement.The high percentage of oil recovery at the miscible flooding indicates that the flooding in the 1D model can be regarded as piston-like displacement.
To demonstrate the temperature change during flow in 1D simulation, the temperature distribution of a non-isothermal numerical run is shown in Figure 6.The injection temperature of CO 2 is 51.2°C.The temperatures of the elements closed to the injection grid element (the first grid shown in Fig. 4) become much lower than that of other elements.The number of the elements with the temperatures lower than original temperatures (99°C) increases with the elapsed time of CO 2 injection.The temperatures of all elements become lower than original temperature at the 69th day.The comparison between the numerical results and experimental data (PV: pore volume at experimental pressure and temperature).
Figure 6 The temperature distribution of 1D simulation.

One Dimensional Immiscible Flooding Sample
The sample is a numerical test with a homogeneous geological model.Three sources of CO 2 at a same rate and four production sinks at a same specified wellbore flow pressure are considered (Tab.4).The geological model and the locations of these sources and sinks are shown in Figure 7.The distance between injection sources and production sinks is constant.
The oil-phase pressures along the one-dimensional model at the 10th, 230th and 415th day are shown in Figure 8.The pressure reaches a maximum value at each injection element (a source location) and it reaches a minimum value at each production element (a sink location).It indicates that the pressure gradients between the injection and production grids reduce with the extension of injection time.The reduction in the pressure gradients may be explained as the decrease in flow resistance, which is caused by the decline of the viscosity of the oil-CO 2 mixture for more CO 2 dissolved in oil.
Figure 9 shows the saturations of three phases along the one-dimensional model after 415 days of CO 2 injection.At the locations (Element 11, 31 and 51) of CO 2 injection, gas saturations reach a maximum value and oil saturations reach a minimum value due to CO 2 injection into these grids at an immiscible condition.The symmetrical distributions of the three phases predicted by the simulator validate the stability of the numerical solution in view of symmetrical check for numerical iteration.Figure 10 shows the gas saturation distributions at three different times.At the 10th day, gas saturation in each grid is zero as the injected CO 2 is completely  The saturations of oil, water and gas at the 415th day.dissolved into the oil phase.Gas (the mixture of CO 2 and light hydrocarbons such as methane, ethane) doesn't form from oil phase until the equilibrium solubility of CO 2 in oil is reached.Continuous CO 2 injection leads to CO 2 and light hydrocarbon's formation from oil and the increase in the numbers of the grids containing free gas.The gas phase almost exists in all grids at the 415th day.
Figure 11 shows the mole fraction of each component in oil phase in each grid at the 132nd day.In the simulation, the hydrocarbons in the crude are grouped into five pseudocomponents ...,.It is assumed that CO 2 dissolves in oil, and water does not dissolve in oil.In the vicinity of the injection elements, CO 2 mole fraction reaches about 0.46 and it does not continuously rise even if more CO 2 injected.The reason is that the maximum mole solubility of CO 2 in the oil at the immiscible conditions is 0.46 and successive CO 2 injected in those grids will form free gas phase.

Two Dimensional Vertical Profile Flooding Sample
In this example, the CO 2 flooding performance in a twodimensional cross sectional model is simulated.The generated grid: is shown in Figure 12.The thickness of the model is 40 m, and its length is 300 m.The horizontal and vertical permeability are 0.16 × 10 -12 m 2 and 0.02 × 10 -12 m 2 respectively.The initial pressure is 20.18 MPa, and the other parameters are listed in Table 4. Figure 13 shows the CO 2 flooding profile at the 500th day.Three zones (regions) are shown in the profiles.Miscible zone locates in the upper region close to the perforated interval.In this zone, no gas phase exists due to the pressure greater than MMP.In immiscible zone, CO 2 does not completely dissolve in the oil.Oil-CO 2 mixture, water and gas phases coexist.In the unswept zone ahead of the immiscible region, only water and oil phases coexist.It implies that the spatial distribution of the three zones and flow characteristics is affected by gravitational difference between two phases.The sample demonstrates that the simulator is capable of simultaneously simulating miscible and immiscible flooding in different regions of an oil reservoir.Phase transfer will carry out automatically at the flooding front.The CO 2 flooding profile at the 500th day since CO 2 injection.

Two Dimensional Horizontal Flooding Sample
In this example, a simulation with nine wells has been performed to study the CO 2 immiscible flooding performances and the effects of asphaltene precipitation on formation damages.The main parameters for this example are shown in Table 5.A CO 2 injection well is located in the center of the square domain.Four production wells are located at four corners and the other four production wells are located at the mid-point of the four outer boundaries.All grids in the x-y plane are squares with same size.The wellbore flow pressure of each production well is assumed constant during the simulation.Figure 14 gives the pressure contours of oil phase in oil formation.The symmetrical distribution of the pressure predicted by the simulator in the plane accords with the theoretical pressure distribution.It indicates that the pressure gradients near the injection well are much less than that of the production wells.The characteristics of pressure gradients imply that the flow resistance of the region near the source (injection well) is much lower than that of the sink locations.The production wells will benefit from the reduction in the flow resistance induced by CO 2 flooding.
The contours of gas saturation are shown in Figure 15.Gas breakthrough occurs at the edge wells (OILwell5 to OILwell8) on the 400th day since CO 2 injection.It is estimated that the area swept by injected CO 2 accounts for less than a half area of the domain after the gas breakthrough at the edge wells.This viscous fingering phenomenon prevents the expansion of CO 2 swept area, which is unfavorable for further EOR in view of corner production wells (OILwell1 to OILwell4).In addition, rectangular grids might sharpen the displacement fronts along the grid orientations (x or y direction), which might also affect the displacement performances.These effects could be reduced by using a nine or seven-point finite difference discretization scheme (Yanosik and McCracken, 1979;Pruess and Bodvarsson, 1983).Hybrid or polygonal grids also could be used to reduce the grid orientation effects.
Figure 16 shows the contours of the CO 2 mole fraction in oil phase after 400 days of CO 2 injection.The CO 2 swept area can be divided into two regions by the dot line (0.46 CO 2 mole fractions).The inner region, (the close region surrounded by the dot line) is saturated by water, oil and CO 2 rich gas.In addition, part of CO 2 dissolves in oil phase.The outer region of the dot line swept by CO 2 is saturated by water and oil  phases.Although the oil containing CO 2 in this region, the CO 2 mole fraction in oil is less than the equilibrium solubility.Therefore, no free gas phase exists in the region.Asphaltene will start to flocculate from oil phase when the mole fraction of CO 2 in oil increases high enough.It may build up on pore walls or at pore throats of oil formation, leading to formation damage.The term permeability change is introduced to evaluate formation damage caused by asphaltene precipitation.It is defined as: Figure 17 shows that the formation damage region expands for the increase in sweep area by injected CO 2 .However, the permeabilities in the location closed to CO 2 injection well are higher than that of the region away from the injection well.One reason is that the high-flow velocity in the region around the well may flush away part of the asphaltene adsorbed on pore surfaces, which leads to a reduction in the capture rates of asphaltene at the sites.The other reason is that the pores in the region closed to injection well is almost saturated by CO 2 .Therefore, only a small amount of solid asphaltene formed because of low oil saturation in the pores of the region.

CONCLUSION
This study was performed to develop a non-isothermal compositional mathematical model for predicting CO 2 miscible and immiscible flooding performances in oil reservoirs for EOR.Physical and chemical processes such as heat transfer, convection-diffusion-adsorption, porosity and permeability changes induced by asphaltene deposition are fully considered during CO 2 flow in porous media.Three scenarios of fluid phases and phase transfer are implemented by checking the MMP and CO 2 equilibrium solubility after the updates of primary variables.Automatic switches of the variables are assigned for the numerical solution according to corresponding phase scenarios.
The governing equations for CO 2 flooding were discretized in space using the integral implicit finite-difference method developed for the TOUGH2 family of reservoir simulators.Newton-Raphson iterative technique was used to solve the nonlinear equation system and a numerical reservoir simulator was developed in FORTRAN code based on the T2R3D simulator belonging to the TOUGH2 family of computer codes.
Four examples demonstrate the various functions of the simulator developed in this work.The comparison of onedimensional miscible simulation results with the experimental data validates the reliability of the mathematical model and the simulator.Symmetrical distributions of the predicted parameters such as pressure, saturations and mole fraction in the second and fourth examples imply that the numerical solution passes the symmetrical test.The displacement front and the interface between miscible and immiscible zones can be traced during the simulation of the third sample.The distributions of temperature and permeability change caused by temperature-independent asphaltene precipitation were also analyzed in the first and fourth examples.
The four examples not only imply that the simulator developed for CO 2 flooding for EOR is robust but also demonstrate the functional variety of phase scenarios.However, further research should focus on a precise correlation (a function of pressure, temperature and components) for predicting CO 2 solubility in oil phases and the effects of geological heterogeneity and grid orientations on flooding performances in actual oil reservoirs.

βββ
Accumulation term of component k per unit bulk volume (kg/m 3 ) A' Parameter in the correlation for calculating viscosity d k Molecular diffusion coefficient of component k within fluid β (m 2 /s) b "Repulsion" parameter in state equation B' Parameter in the correlation for calculating viscosity c Parameter in state equation C k s Concentration of component k in gas phase (kg/m 3 ) C perm Permeability change C' Parameter in the correlation for calculating viscosity D k Hydrodynamic dispersion tensor accounting for both molecular diffusion and mechanical dispersion for component k in phase β (m 2 /s) D' Parameter in the correlation for calculating viscosity f Flow efficiency factor f RF Mole fraction of components, C 2 -C 4 F k "Flow" term of mass or diffusive and dispersive mass transport per unit bulk volume (kg/(s.m 3 )) g Gravitational constant (m/s 2 ) h β Specific enthalpy of fluid phase β (J/kg) h k Specific enthalpy of component k in fluid β (J/kg) Equilibrium partitioning coefficient of component k between phases α and β MMP Minimum Miscibility Pressure (Pa) Adsorption term of component k on rock solids P C Capillary pressure (Pa) P β Pressure of phase β (Pa) q k Figure2CO 2 equilibrium solubility in the crude oil (experimental data).

Figure 4
Figure 4Model grid of CO 2 flooding for the numerical simulation example.I: CO 2 injection; P: oil production sink.

Figure 8
Figure 8The oil phase pressure distribution at different injection time.
Figure 10Gas saturations at the different time.

Figure 11 Mole
Figure 11Mole fraction distribution in crude oil at the 132nd day.
Figure 12Model grid for the two-dimensional vertical section model.
The distribution of pressure and temperature, phase saturations, mole fraction of each component in each phase, formation damage caused by asphaltene precipitation and the improved oil recovery are predicted by the simulator.Experimental data validate the developed simulator by comparison with simulation results.The applications of the simulator in prediction of CO 2 flooding in oil reservoirs indicate that the simulator is robust for predicting CO flooding performance.
β Mole fraction of component k in phase β y Coordinate at Y direction (m) z Coordinate at Z direction (m)

Table 3 .
The MMP of the CO 2 -oil system is 27.45 MPa, which is obtained by slim tube experiments.The crude oil used for the slim tube and CO 2 flooding experiments sampled from H75-27-7 well, H.G. Oil field, CNPC.Its viscosity and density at the reservoir condition are 1.86 mPa.s and 0.7542 g/cm 3 .The initial pressure of each grid element and flow pressure at the outlet (the 10th element) is set at 32.27 MPa and 32.26 MPa, respectively, which results in miscible flooding during the flooding.The conditions such as pressure and temperature are the same in the numerical runs and experiments.

TABLE 4
Main parameters of one-dimensional immiscible flooding simulationModel grid of CO 2 flooding for one dimensional simulation example.I: CO 2 injection source; P: oil production sink; the bulk of each block is 1 m 3 .

TABLE 5 Main
Binshan Ju et al. / Computer Modeling of the Displacement Behavior of Carbon Dioxide in Undersaturated Oil Reservoirs963Figure16CO 2 mole fraction in oil at the 400th day after CO 2 injection.The distribution of permeability change caused by asphaltene precipitation.a) At the 150th day after CO 2 injection; b) at the 400th day after CO 2 injection.