Experimental Study of Carbon Dioxide Diffusion in Oil-Saturated

Aug 27, 2009 - The gas diffusion coefficient is determined by fitting the amount of gas diffused into the ... The density gradient can induce natural ...
0 downloads 0 Views 2MB Size
Download PDF
Ind. Eng. Chem. Res. 2009, 48, 9307–9317

9307

Experimental Study of Carbon Dioxide Diffusion in Oil-Saturated Porous Media under Reservoir Conditions Zhaowen Li†,§ and Mingzhe Dong*,‡ Faculty of Engineering, UniVersity of Regina, Regina, Saskatchewan, S4S 0A2, Canada, and Department of Chemical and Petroleum Engineering, UniVersity of Calgary, 2500 UniVersity DriVe NW, Calgary, Alberta, T2N 1N4 Canada

For both CO2 enhanced oil recovery and CO2 sequestration in oil reservoirs, the diffusion of injected CO2 into oil-saturated porous media is of great significance in project design, risk assessments, economic evaluation, and performance forecast. However, the measurement of CO2 diffusion coefficient in liquid-saturated porous media under reservoir conditions has not been well-established due to the complications caused by densityinduced natural convection in the CO2-oil systems. This paper presented a new method for measuring the effective CO2 diffusion coefficient in oil-saturated porous media under reservoir conditions. A core plug with the two end faces sealed was designed as the physical model for radial diffusion process. The measurement was conducted in a high-pressure diffusion cell with an oil-saturated core sample placed in the middle and the remaining void space of the cell filled by high-pressure CO2 sample. A small-pressure decay technique was employed to record the pressure change of the gas phase during the diffusion measurement. Because the oil phase contained in the porous medium swells as CO2 diffuses into it, the experiment essentially involves both diffusion and swelling-induced convection. To describe these processes involved in the measurement, a mathematical model along with numerical solutions was derived. The effective diffusion coefficient was determined by matching the experimental pressure decay curve with the corresponding mathematical model. The proposed method was demonstrated for oil-saturated Berea core samples at pressures ranging from 2.2 to 6.5 MPa. The measured pressure decay curves showed good agreement with the mathematical predictions using the best-fitted effective diffusion coefficients. The derived method can be readily implemented in laboratories that can handle high-pressure fluids; thus, this study provides a tool for studying CO2 diffusion in oil-saturated porous rocks. Introduction Gas injection into oil reservoirs to enhance hydrocarbon recovery has been widely practiced in the oil industry since the beginning of the last century. More recently, CO2 geological sequestration has been considered as one of the most promising options to mitigate the dramatic increase of CO2 concentration in the atmosphere due to the booming fossil energy consumption since the last century. In both CO2-EOR (enhanced oil recovery) and CO2 geological sequestration processes, reliable data of gas effective diffusion coefficients in liquid-saturated porous rocks under reservoir conditions are required for project design, risk assessment, economic evaluation, and performance forecast. Unfortunately, studies of CO2 diffusion in liquid-saturated porous media under reservoir conditions are very limited due to the complications of high-pressure CO2-liquid systems. For instance, one of the most notorious difficulties is the natural convection caused by the density increase in the liquid phase as a result of CO2 dissolution, which significantly interferes with the diffusion process and makes the conventional experimental schemes unsuitable to measure the CO2 diffusion coefficient in bulk liquids at high pressures.1-5 Much research on the gas diffusion coefficient in bulk liquids has been conducted since the 1930s.2,6-10 Most of them were conducted using a PVT (pressure-volume-temperature) cell, also called an adsorption or diffusion cell. This method is often * To whom correspondence should be addressed. Tel.: 1 (403) 2107642. Fax: 1 (403) 284-4852. E-mail: [email protected]. † University of Regina. ‡ University of Calgary. § Currently with Husky Energy Inc., Calgary, Canada.

referred as the PVT method. In this method, either the change in pressure of the system at a constant volume or the change in gas-phase volume at a constant pressure, or both, is recorded as a function of time. The amount of gas diffused into the liquid phase is calculated from the change in pressure or volume of the gas phase with the equation of state (EOS) of real gases. The gas diffusion coefficient is determined by fitting the amount of gas diffused into the liquid phase (or the change in pressure or volume of the gas phase) versus the time curve with the corresponding diffusion models. The PVT method has been successfully used to measure the diffusion coefficients of different gases in liquids under pressures, such as hydrocarbon gases and nitrogen in oil and water. Because the dissolution of these gases into the liquid reduces the liquid-phase density, there is no concern of density-induced natural convection. However, the PVT method cannot be used for the majority of the CO2-liquid systems. This is because the dissolution of CO2 in liquids, such as oil and water, results in an increase in density of the liquid phase and there will be a density gradient along the distance from the gas/liquid interface to the liquid phase. The density gradient can induce natural convection.1-3 As a result, the mass transfer of CO2 into liquid in the PVT method is interfered with or even dominated by the density-induced convection unless the viscosity of the liquid phase is sufficiently high such that the gravity force caused by the density difference can be balanced by the viscous force.3 For this reason, the results of the diffusion coefficient of CO2 in reservoir liquids obtained with the conventional PVT method may have significant errors. Instead of using a diffusion cell, Grogan et al. proposed a glass capillary tube to measure the diffusion coefficient of CO2

10.1021/ie900145c CCC: $40.75  2009 American Chemical Society Published on Web 08/27/2009

9308

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

in light oils under high pressures.7 The capillary tube was enclosed in a windowed high-pressure cell filled with overburden fluids. To reduce the density-induced convection, the capillary tube was positioned horizontally. A CO2 gas slug was created in the middle of the capillary tube filled with oil sample. The positions of the two CO2/oil interfaces in the tube were recorded as a function of time for the diffusion coefficient calculation using a numerical model. As pointed out by the researchers themselves,7 cautious operations were required in the experiment to reduce the convection impact on the measurements. Although this method can be used to measure the diffusion coefficient of CO2 in liquid phase, it is difficult to apply the measured results to liquid-saturated porous media before the characteristics of gas diffusion in liquid-saturated porous media are wellunderstood. Renner used liquid-saturated Berea cores to measure the effective diffusion coefficient of CO2 in porous media.2 In his measurement, the side surface and one end face of the Berea core were sealed with epoxy resin and fiberglass tapes before it was saturated with the test liquid. The gas diffused into the Berea core only from one end face. Due to the limitation of the crosssectional area of the core sample, the diffusive flux of CO2 through the end face was small; therefore, in order to make the pressure decay measurable, a small gas volume has to be used in the measurement. This made the pressure decay process very sensitive to the fluctuations of pressure and temperature change of the environment during the experiment. Renner also found that the diffusion coefficient measured when the core was positioned vertically was much larger than that when the core was positioned horizontally due to the effect of density-induced natural convection in the former case where CO2 diffused into the core sample from the top.2 There are two most commonly seen definitions of diffusive flux in porous media in the literature:11-14 J ) -Deff∇cr

(1.a)

J ) -Deff ′ ∇c

(1.b)

and

where J is the diffusive flux based on the area of a porous medium; c and cr are the concentrations of the diffusing species in liquid phase and in the bulk volume of the porous medium, respectively. Deff and D′eff are the two different effective diffusion coefficients of the gas in the porous medium. The two concentrations, cr and c, are related by cr ) φc

(2)

where φ is the porosity of the porous medium. The two different ′ , are defined, effective diffusion coefficients, Deff and Deff respectively, as follows: Deff ) D/ε

(3.a)

Deff ′ ) φD/ε

(3.b)

and

where D is the diffusion coefficient in bulk liquid phase; ε is the diffusive tortuosity factor of the liquid-saturated porous medium. Given that the diffusion coefficient of a gas in a liquid and the diffusive tortuosity factor of a porous medium are available, the effective diffusion coefficient can be calculated by eq 3.a or eq 3.b. However, the diffusive tortuosity factor of a porous medium is difficult to obtain without the diffusion

measurements using the porous medium. The values of tortuosity factors obtained from other methods might be different from the diffusive tortuosity factors.14 Krooss and Schaefer15 and Krooss and Leythaeuser16 developed an experimental approach to measure effective diffusion coefficients of hydrocarbon gases in water-saturated sedimentary rocks at atmospheric pressure and different temperatures. In their experiments, a thin rock slice (3.10 mm thick) cut from a rock plug was used. The amount of gas diffused across the watersaturated rock sample was monitored through a gas chromatography, and the diffusion coefficient was obtained through the transient diffusion model in a plane sheet. They pointed out that great care had to be taken during the sample preparation and installation, especially for the poorly consolidated rocks.15 They also indicated that their experimental setup was not suitable for measurements under high pressures because a small pressure fluctuation in the gas chamber can cause gas flow through the thin porous plate. Although attempts have been made to measure gas effective diffusion coefficients in oil-saturated porous rocks, a reliable method for high-pressure conditions has not been developed. The objective of this work is to develop an experimental method and mathematical solutions to determine gas effective diffusion coefficients in oil-saturated porous media. Berea core samples which are widely used in petrophysical studies are used as the porous media in this study. The experiment is designed to allow the change in pressure (for a constant gas-phase volume) caused by diffusion to be observable and less sensitive to the environment at the same time. The technique should also be easily implemented in laboratories. To meet those requirements, a core plug with the two end faces sealed was designed as the physical model so that the gas diffuses into the sample only from radial directions. A small-pressure decay technique as used in the PVT method was adopted to record the pressure change of the gas phase during the measurement. This physical design provides a much larger diffusion area and, thus, a more observable diffusive flux than that of Renner’s experiment.2 Moreover, the surface area can be adjusted simply by changing the length of the core sample. To demonstrate this method, the effective diffusion coefficients of CO2 in oil-saturated Berea cores are measured under different pressures. It should be mentioned that more details regarding the experimental design and procedure have been provided in a previous study of CO2 effective diffusion coefficient measurements in water-saturated porous media.5 The major difference between CO2-water and CO2-oil systems is that the swelling effect caused by CO2 dissolution is negligible in the former but significant and has to be considered in the latter.17 Therefore, this work will focus on model development for a swelling system. Experimental Section Physical Model. In the designed physical model shown in Figure 1, a core plug with end faces sealed is used in the measurement. Since the two end faces of the core sample are sealed, the gas diffusion occurs only along the radial direction as illustrated in Figure 1. The side surface has a larger area which can be increased easily by increasing the length of the core sample to allow more measurable diffusive flux of gas into the porous medium. In this experiment, the core sample is positioned vertically in a diffusion cell which is placed in a gas or water bath. The change in gas-phase pressure in the cell is recorded as a function of time to generate a pressure decay curve which is used for calculating the effective diffusion coefficient.

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

Figure 1. Schematic of gas diffusing into an oil-saturated porous column along radial direction.

As discussed in the previous study,5 the advantages of proposed design are (1) the porous column (rock plug) used in the measurement is easier to prepare and install; (2) the side surface provides a greater diffusion flux area and, thus, the measurement is easy to operate and less sensitive to the working environment; (3) during a measurement, the gas, rock sample, and liquid are always under the same pressure; hence, this method can be used under high pressures to simulate actual reservoir conditions. Apparatus. A schematic diagram of the experimental setup is shown in Figure 2. It consists mainly of a diffusion cell, a high-pressure gas cylinder, a high-pressure positive-displacement pump (hand pump), and a data acquisition system. The high-pressure diffusion cell (with a diameter of 6.3 cm and depth of 22 cm) was used to hold the rock plug and the test gas during the experiment. The high-pressure gas cylinder (500 cm3) was used to supply pressurized CO2. Both the diffusion cell and the high-pressure cylinder were placed in the water bath where the temperature was maintained at a desired value with an accuracy of (0.1 °C. A pressure stability test showed that the change in pressure caused by the surrounding temperature fluctuations was within (1 kPa for a period of 5 h when the pressure was set at 10 MPa. The high-pressure hand pump was used to pressurize the CO2 and, with the assistance of a backpressure regulator (BPR), to maintain a constant system pressure when transferring CO2 from the gas cylinder to the diffusion cell. A flask was connected to the bottom of the diffusion cell to collect the oil displaced when loading CO2 into the diffusion cell. The pressure of the diffusion cell was recorded by a Heise digital pressure gauge with a resolution of 1 kPa. Since CO2 has higher solubility in oil than in water, a larger gas-phase chamber is desirable to ensure small pressure decay in a measurement. It should also be mentioned that a relatively large free chamber below the core sample was created by a slotted steel ring to collect the swelled-out oil from the core column as CO2 diffusion proceeds. Procedure. Before a diffusion measurement, the entire system was tested for leakage using nitrogen at about 10 MPa. Cores were cut into the designed lengths as the test samples, cleaned, and then dried in an oven at 110 °C for 3 days. After that, gas permeability of core sample was measured by a PDPK-400 permeameter.18 Epoxy resin was applied to seal the two end faces of the dried core samples, and then the core sample was put into the desiccator for over 1 day to allow the epoxy resin to completely solidify. Following that, the core sample was

9309

weighed, vacuumed, and saturated with oil. The porosity of the core sample was calculated from the weight difference before and after it was saturated with oil. After the above preparations, the following procedure was applied for each measurement: (1) Place the slotted stainless steel ring on the bottom of the cell to support the oil-saturated core sample. Load the core sample and position it vertically in the middle of the cell. Tighten the cap of the cell and inject the oil sample from the bottom of the diffusion cell to displace the air in the cell. (2) Place the diffusion cell and the gas cylinder into the water bath, and wait for over 1 day for the temperature of the system to stabilize at the test temperature. Ensure the valves connecting the diffusion cell and gas cylinder are closed during the warmup stage. Adjust the pressures of the gas cylinder and diffusion cell to a designed test pressure. (3) Vacuum all the lines connecting the gas cylinder, the pressure gauge, and the diffusion cell. Carefully open the valves connecting the diffusion cell and gas cylinder. Introduce CO2 sample into the diffusion cell from the top of the cell to displace the free oil (oil in the cell but not contained in the core) sample in the cell. Control the drainage valve at the bottom of the diffusion cell with the assistance of the BPR to maintain the pressure of the system at the test pressure for the entire CO2 transfer process. Collect the drained oil sample with the flask to determine the volume of the gas phase inside the diffusion cell. (4) After the free oil in the cell is drained, close the valve at the bottom of the diffusion cell and the valve to the CO2 cylinder. Start to record the cell pressure immediately to perform pressure decay measurement. (5) After a pressure decay measurement, release the pressure of the diffusion cell. Detach, open, and clean the cell for the next measurement. A fresh core sample is used for each measurement. Materials. Six core samples were cut from three long Berea core plugs (about 30 cm long each) for the CO2 effective diffusion coefficient measurements. The dimensions and properties of the six core samples are listed in Table 1. Carbon dioxide (CO2) used in the experiment was supplied by Praxair (Regina, Canada) with a purity of 99.9%. NHexadecane (92+ %) (Alfa Aesar, U.S.A.) was used as the oil sample. The molecular weight of the N-hexadecane (92+ %) was 226.45 g/mol, and the density was 0.7733 g/cm3 at 18 °C and 0.7596 g/cm3 at 40 °C.19 The viscosity of the oil was measured to be 3.35 mPa · s at 20 °C and 2.14 mPa · s at 40 °C. The measured viscosity values are rather close to those of pure N-hexadecane under the same conditions.20 Theoretical Development Assumptions. The following assumptions are made for smallpressure decay tests: (1) The CO2 effective diffusion coefficient is constant over the range of the pressure decay measurement. (2) The CO2 concentration in the oil phase at the core surface (r ) ro) is constant in the pressure decay measurement. (3) The core sample is homogeneous. (4) The density-induced natural convection is negligible. (5) The effect of oil evaporation into the gas phase is negligible. (6) The oil and gas are not miscible. Assumption 1 requires the effective diffusion coefficient to be independent of both the gas concentration in oil and the

9310

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

Figure 2. Schematic of the experimental setup for effective diffusion coefficient measurement. Table 1. Dimension, Porosity, and Permeability of Six Core Samples Berea plug

diameter (cm)

permeability (md)

porosity (%)

test sample

length (cm)

Ber-1

5.15

263

18.85

Ber-2

5.05

163

18.20

Ber-3

5.15

160

18.32

Ber-1-1 Ber-1-2 Ber-2-1 Ber-2-2 Ber-3-1 Ber-3-2

9.90 10.25 10.15 9.74 10.21 9.98

pressure during a measurement. The former is usually valid for the dilute solutions, and the latter requires that the pressure drop in a decay test is small. Assumption 2 is acceptable if the pressure change in a pressure decay measurement is small. For the proposed CO2/ oil-saturated porous column, the volume of oil expands as CO2 diffuses into the oil phase. The CO2/oil solution that swells out of the column forms a thin oil film on the column surface and flows downward along the surface to the chamber beneath the core sample. Since the oil film on the core surface is thin and exposed to CO2 phase, it is reasonable to assume that the oil concentration at the column surface is at equilibrium with CO2 phase under the test pressure and temperature. The homogeneity of the core samples were confirmed by measuring the permeabilities along and around the core sample with the PDPK-400 permeameter. When a homogeneous porous column with end faces sealed is positioned vertically, the CO2 concentration profile develops along the radial direction through diffusion, and thus, there is no density gradient in the vertical (gravitational) direction. The only natural convection caused by radial density gradients as CO2 diffuses into the oil phase has been demonstrated to be negligible based on the calculated Rayleigh number Ra e 10.4 The solubility of both CO2 in hexadecane and hexadecane in oil under different pressures was calculated by the tuned EOS with CMG WinProp.21 To tune the EOS, GOR (gas-to-oil ratio) and swelling factor were measured under different pressure levels and were used as input data for Winprop. The default

parameters in Winprop were used for tuning the EOS. More details of the binary CO2-hexadecane systems have been given elsewhere.22 The maximum solubility of hexadecane in CO2 phase was calculated to be 0.12 mol % at the maximum pressure of 8 MPa used in the measurements. The maximum mass of the hexadecane that could be extracted into the CO2 phase for the designed experiment was determined to be approximately 0.27 g, which is negligible. Mathematical Model of Diffusion and Swelling-Induced Convection. The mass transfer of CO2 in an oil-saturated porous column can be considered as a process of diffusion combined with bulk convection (volume expansion) along the radial direction in a one-dimensional radial system with fixed boundary location and constant boundary concentration. The conservation equation of the concentration for diffusing species bears the form of ∂c ) Deff∇2c - ∇ · (cv) ∂t

(4)

where v is the volume average velocity caused by the swelling of oil phase due to gas diffusion. By applying ∇ · (cv) ) c(∇ · v) + v · ∇c

(5)

eq 4 can be written into ∂c ) Deff∇2c - c(∇ · v) - v · ∇c ∂t

(6)

Noting that both the velocity and the concentration gradients in the θ and z directions are zero, eq 6 can be further written in scalar form as c∂ ∂c ∂c + (ru) + u ) Deff∇2c ∂t r ∂r ∂r

(7)

Equation 7 is a nonlinear partial differential equation (PDE) because the coefficient of (∂c/∂r) (velocity u) depends on the

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

9311

concentration development of the diffusing species in the porous column. For a nonswelling system, the velocity equals zero (no convection), and then eq 7 is simplified into Deff ∂ ∂c ∂c ) r ∂t r ∂r ∂r

( )

(8)

Equation 8 is the radial diffusion equation. The analytical solution of eq 8 with its application in determining the gas diffusion coefficient for a nonswelling system has been discussed in the previous study.5 The effort of this section is to derive a numerical solution to eq 7. The boundary conditions for such a system are c ) c o,

at r ) ro ∂c ) 0, at r ) 0 ∂r

u ) 0,

(9)

The initial conditions are u ) 0, u ) 0,

c ) c i, c ) c o,

for r < ro for r ) ro at t ) 0

(10)

For the diffusion coefficient measurement in this study, the liquid phase that saturates the porous rock column is initially free of gas, i.e., ci ) 0. Numerical Solution. By introducing the reference units ro2/Deff, ro, Deff/ro, and co - ci for time, length, velocity, and concentration, respectively, a new set of dimensionless variables are defined as τ)

t , ro /Deff 2

jr )

r , ro

uj )

u , Deff /ro

jc )

c - ci co - ci (11)

Then, eq 7 can be reduced to the following dimensionless form jc ∂ ∂cj ∂cj + (rjuj) + uj ) ∇2jc ∂τ jr ∂rj ∂rj

(13)

where λ ) uj -

1 jr

derivative of time. A Gauss-Seidel iteration scheme was used to calculate the concentration and velocity distributions for each time step. For each new time step, the concentrations at the previous time step were used as the initial values to calculate the velocity distributions (the calculation of velocity from concentration field is derived in the following section) for the new time step. With the new calculated velocities, the concentration field was updated. The velocities were then computed again with the updated concentration field. The above iteration was repeated until the maximum relative error of the concentration at each radial node was less than the error tolerance e (e ) 10-4 in this study). Calculation of Convective Velocity Caused by Volume Swelling. In the calculation of convective velocity caused by volume swelling, it is assumed that the dissolution of gas causes the swelling of the oil phase and consequently a radial expansion of the oil phase toward the surface. This expansion will generate a velocity pointing to the surface. For a given annular volume element from r to r + ∆r of an oil-saturated porous column, as shown in Figure 3, the volume increase of the oil phase due to swelling within a time interval ∆t is

[(

(12)

For the convenience of discretization, the term ∇2jc is expanded, and eq 12 can be rearranged into ∂uj jcuj ∂cj ∂2jc ∂cj + jc + +λ ) 2 ∂τ ∂rj jr ∂rj ∂rj

Figure 3. Annular volume element of an oil-saturated porous column.

∆Ve ) φh[π(r + ∆r)2 - πr2] f

c ) 0 for jr < 1 c ) 1 for jr ) 1 at τ ) 0

)]

(17)

(14) ∆ur+∆r,t+∆t )

∆Ve 2π(r + ∆r)h∆t

(18)

Substitution of eq 17 into eq 18 yields (15) ∆ur+∆r,t+∆t )

The initial conditions, given by eq 10, are reduced to uj ) 0, uj ) 0,

(

where ∆Ve is the volume increase of the oil phase within the circular element and f is the swelling factor as a function of concentration, which is determined experimentally. An averaged velocity generated at r + ∆r by the volume swelling of this circular element within ∆t is

The boundary conditions, given by eq 9, are also reduced to jc ) 1 at jr ) 1 ∂cj ) 0 at jr ) 0 uj ) 0, ∂rj

)

cr,t+∆t + cr+∆r,t+∆t 2 cr,t + cr+∆r,t f 2

(16)

The dimensionless diffusion-convection equation, eq 13, is solved by fully implicit finite difference method (FDM). In discretization, central difference approximations were used for both the first- and second-order derivatives of concentration and velocity, and forward difference was used for the first-order

[(

)

cr,t+∆t + cr+∆r,t+∆t φr∆r f (r + ∆r)∆t 2 cr,t + cr+∆r,t f 2

(

)]

(19)

The actual velocity at r + ∆r and time t + ∆t is the summation of the velocities generated by all the elements from r ) 0 to r + ∆r, i.e., r+∆r

ur+∆r,t+∆t )

∑ ∆u

r+∆r,t+∆t

r)0

(20)

9312

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

By applying the reference units ro2/Deff, ro, Deff/ro, and co ci for time, distance, velocity, and concentration, respectively, eqs 19 and 20 can be converted into the following dimensionless forms: ∆ujjr+∆rj,τ+∆τ )

[(

)

jcjr,τ+∆τ + jcjr+∆rj,τ+∆τ φrj∆rj f (co - ci) (rj + ∆rj)∆τ 2 jcjr,τ + jcjr+∆rj,τ f (co - ci) (21) 2

(

)]

r+∆r

ujjr+∆rj,τ+∆τ )

∑ ∆uj

jr+∆rj,τ+∆τ

(22)

r)0

For each time step, eq 22 is used to calculate the dimensionless convective velocity in eq 13. Calculation Procedure of the CO2 Effective Diffusion Coefficient in Oil-Saturated Porous Media. The mass conservation of the diffusing species is used to relate the amount of gas lost in the gas phase and the amount of gas entered into the oil phase. Here, the oil phase includes both the oil contained in the porous column and the oil swelled out of the column caused by CO2 dissolution. As the oil swells into the gas phase, the volume of the gas phase decreases. With the EOS, the following equation is obtained for calculating the amount of gas lost in the gas phase, given the pressure drop, ∆P, in the gas phase and the volume of the oil swelled out are known: ∆ng )

∆PV + ∆VPo - ∆V∆P ZRT

(23)

where Z is the gas deviation factor, R is the gas constant, T is the absolute temperature, ∆ng is the amount of gas lost in the gas phase, Po is the initial pressure of the gas phase when the gas diffusion into the porous column starts, and ∆V is the volume decrease of the gas phase in the diffusion cell, or the volume of the oil swelled out of the column. ∆V can be calculated from the velocity of the oil phase at the surface of the porous column using the following equation: τ

∆V ) 2πr2ohφ

∑ uj

jr)1,τ∆τ

(24)

τ)0

From the mass balance, the amount of gas lost in the gas phase, ∆ng, given by eq 23, equals the amount of gas diffused into the oil phase, q, at time t. Here, q includes the gas dissolved in both the oil contained in the porous column and the oil swelled out. Replacing ∆ng in eq 23 with q, eq 23 can be rearranged into ∆P )

qZRT - ∆VPo V - ∆V

Figure 4. Comparison of concentration profiles from numerical and analytical solutions for a nonswelling system at different times.

(25)

Equation 25 is used to match the plot of experimentally measured pressure drop ∆PEx versus t1/2 thereby determining the CO2 effective diffusion coefficient in an oil-saturated porous column. q and ∆V are determined thorough numerical solutions. Results and Discussion Validation and Analysis of the Numerical Solution. Two attempts were made to examine the numerical solution and computer code. First, the numerical solution was compared with the analytical method for a nonswelling case which has an analytical solution, i.e., the analytical solution of eq 8.5 Second, the impact of swelling was examined through analyzing the

Figure 5. Dimensionless velocity profiles from numerical solution for a system with swelling factor of 1.26.

velocity distribution for a swelling case and comparing concentration profiles and adsorption rates of systems with different swelling factors to the nonswelling base case. For the nonswelling case, the selected concentration profiles from both the numerical and analytical solutions are plotted in Figure 4. It can be seen that the numerical solution agrees with the analytical method well throughout the entire domains of distance and time (0 e jr e 1; 0 e τ e 1). The dimensionless velocity distribution for a swelling case with swelling factor of 1.26 is plotted with respect to dimensionless radius at different times as shown in Figure 5. The peak velocity occurs near the surface of the porous column at the very early stage when the concentration gradient is the sharpest and thus the diffusive flux is the fastest. As the diffusion proceeds, the peak velocity point moves toward the center of the column as the sharpest concentration gradient develops toward the center. Another reason contributing to the decrease in velocity toward the surface is that the fluxing area becomes larger toward the surface, and thus velocity is reduced for the same volume expanded fluids passing through the larger area. Selected concentration profiles calculated numerically at dimensionless time τ ) 0.1 and τ ) 0.2 for different swelling factors are shown in Figure 6. Here, f is the swelling factor defined as the ratio of the volume of oil containing the dissolved gas to the volume of oil without gas under the same pressure and temperature conditions. It is shown that the volume swelling-induced convection tends to retard the gas mass transfer into the porous column and a greater swelling factor yields a greater retarding effect on mass transfer rate. Physically, this

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

Figure 6. Comparison of concentration profiles for systems with different swelling factors.

Figure 7. Comparison of the gas absorption rate by the oil phase relative to a nonswelling system, q/qp vs τ1/2, with different swelling factors.

can be explained by that the direction of the velocity induced by volume swelling is opposite to the diffusion and, thus, tends to slow down the concentration development. To compare the models with and without swelling effect, the gas absorption rate by the total oil phase relative to a nonswelling reference case is necessary. The total gas that can be absorbed by the oil in the porous column at infinite time for the nonswelling case is denoted as qp. The dimensionless amount of gas absorbed (relative to its corresponding nonswelling case) at time t is denoted as q/qp. Figure 7 shows the numerically calculated q/qp versus τ1/2 for two swelling cases (f ) 1.26 and f ) 1.69) compared with a nonswelling case (f ) 1.0). It can be found that the gas transfer rate into the total oil phase is higher for the swelling case (f ) 1.26 and f ) 1.69) than that of the nonswelling case (f ) 1). A greater swelling factor yields a faster gas transfer rate and a greater cumulative amount of gas transferred into the oil phase. This is because the undersaturated oil layer near the column surface is pushed to the surface of the porous column by swelling, where the saturation with CO2 is instantaneous. Note that q/qp reaches unity at τ1/2 ) 1 for the nonswelling case. However, for the swelling cases of f ) 1.26 and f ) 1.69, q/qp approaches, but is slightly less than, 1.26 and 1.69, respectively, at τ1/2 ) 1. At time τ1/2 ) 1, q/qp ) 1.253 for the swelling case with f ) 1.26 and q/qp ) 1.682 for the case with f ) 1.69. This means that although a greater swelling factor yields a faster mass transfer rate of gas into oil phase, it takes a longer time for the total oil phase to become saturated because

9313

the swelling tends to retard the gas entering into the porous column, which is shown in Figure 6. CO2 Effective Diffusion Coefficient in Oil-Saturated Rock Samples. Six measurements of CO2 effective diffusion coefficient (EDCO2/O-1 to EDCO2/O-6) were conducted at 40 °C and at pressures ranging from 2.3 to 6.3 MPa. The recorded pressure decay curves for measurements EDCO2/O-1 through EDCO2/O-6 are shown in Figure 8a-f. The effective diffusion coefficient, Deff, is obtained when the experimental plot of ∆P versus t1/2 was matched by eq 25 of which q and ∆V are determined by numerical solutions. The very early part of the ∆P versus t1/2 plot from experiment was slightly deviated from the theoretical prediction due to the initial effect in the measurement, which has also been found by other researchers.2,3,23 This small deviation from the mathematical prediction at the very early few minutes can be explained by (1) the mathematical model assumes an instantaneous equilibrium condition at the CO2/oil interface, while time is needed in the experiment for the CO2/oil interface to reach equilibrium after CO2 is loaded; (2) the experimental operations including stopping the pump, closing the valve, and starting to record pressure decay also cause certain deviations. Therefore, this early deviation should be considered as a system error of the experiment. The early part of the experimental plot is disregarded, and Deff is determined with the remaining part of the experimental plot. To make it more explicit, the experimental plot can be easily corrected to exclude the initial effect by subtracting a small pressure difference from the measured values data, and then the experimental plot was matched by the predicted one. The details of this procedure have also been given in the previous study.5 The experimental, predicted, and corrected curves of pressure drop ∆P versus square root of time t1/2 for all six measurements are shown in Figure 9a-f, respectively. Before determining the effective diffusion coefficient Deff using the numerical solution, an approximate value of Deff can be estimated first through the analytical solution without considering the swelling-induced convection, i.e., the analytical solution of eq 8. The results determined are listed in the column “analytical” in Table 2. The detailed procedure of determining the value of Deff for a nonswelling system has been provided in the previous study.5 This approximated value of Deff is then used as the initial value in the numerical matching process. An ultimate effective diffusion coefficient Deff is obtained when the corrected experimental plot is best matched by eq 25. The results of the measurements for CO2 in six oil-saturated rock samples are summarized in Table 2. The concentration of CO2 at the core surface, co, and the gas compressibility factor, Z, are evaluated at the time-averaged pressure of each pressure decay measurement. In Table 2, the column of f(co) gives the swelling factors of the oil at concentration co (the concentration at the surface of the rock column). The PVT properties of the CO2-oil system, including co, f(co), and Z are calculated with PengRobinson equation of state (CMG WinProp, 2004) which was tuned with the measured data.22 The total amount of CO2 diffused into the core (or lost in the gas phase) was calculated using eq 23 and is listed in the column “∆ng” of Table 2. In the column “Deff”, the analytical results were calculated using the analytical solution of eq 8, whereas the numerical results were obtained through eq 25 where q and ∆V are determined thorough the numerical solutions derived in the previous section. It is shown in Table 2 that the measured effective diffusion coefficient falls in the range of 6 × 10-6 to 8 × 10-6 cm2/s. For each pair of core samples cut from the same core plug, the diffusion coefficient increases slightly with increasing pressure.

9314

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

Figure 8. Recorded pressure decay curves for measurement EDCO2/O-1 through EDCO2/O-6.

However, Deff for measurement EDCO2/O-6 is higher than that of both measurements EDCO2/O-3 and EDCO2/O-4; even the test pressure for EDCO2/O-6 was lower. This is not surprising because the core samples used in these experiments are different. The core sample used in measurement EDCO2/O-6 was cut from Berea plug Ber-3, whereas the core samples used in measurements EDCO2/O-3 and EDCO2/O-4 were cut from Ber-2. They may have different tortuosity factors. The results in Table 2 may suggest that the CO2 effective diffusion coefficient in the tested oil is only weakly dependent on pressure within the pressure range used in this study. Renner measured the diffusion coefficients for CO2 in liquid decane under conditions similar to those in this study with a pressure range of 1.5-5.86 MPa and temperature at 38 °C.2 The reported results varied from 1.97 to 5.05 × 10-5 cm2/s. If it is assumed that the core samples used in this work have similar diffusive characteristics and the diffusion coefficient for CO2 in decane is close to that in

hexadecane, the diffusive tortuosity factors of the Berea cores used in this work are estimated to be in the range of 3-6 with eq 3.a. Examination of Effects of Pressure Variation, Natural Convection, and Swelling. Pressure Variation. The maximum pressure drop of all the pressure decay measurements was less than 6.8% of the time-averaged pressure. The maximum relative error caused by the solubility and compressibility factor variation due to the pressure decrease in a measurement was estimated to be less than 7.4%. The measured results also show that the effective diffusion coefficient may be only slightly dependent on the pressure. During a pressure decay measurement, the relative change of the effective diffusion coefficient is less than 2.5%. Therefore, assuming a constant effective diffusion coefficient in a small-pressure decay process is reasonable. Natural Convection. To examine the effect of natural convection on the determined effective diffusion coefficients,

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

9315

Figure 9. Experimental and theoretical ∆P vs t1/2 plots for measurements EDCO2/O-1 through EDCO2/O-6. Table 2. Summary of Measurements of the CO2 Effective Diffusion Coefficient in Oil-Saturated Rocks Deff (10-6 cm2/s)

pressure (kPa) 3

measurement no.

sample

initial

final

average

V (cm )

EDCO2/O-1 EDCO2/O-2 EDCO2/O-3 EDCO2/O-4 EDCO2/O-5 EDCO2/O-6

Ber-1-1 Ber-1-2 Ber-2-1 Ber-2-2 Ber-3-1 Ber-3-2

2369 6518 4416 5015 5850 3175

2242 6138 4165 4838 5484 2968

2283 6259 4244 4925 5600 3031

188.64 180.81 164.23 189.72 191.56 142.63

-3

c0 (10

the Rayleigh numbers for the experimental conditions were calculated for all the measurements. The range of the calculated Rayleigh numbers varies from 0.4 to 0.7. The effect of natural convection on the diffusion process is negligible for such low Rayleigh number.4 Swelling. The oil-phase swelling caused by gas dissolution has two opposite effects on the pressure decay process compared to a nonswelling base case. On one hand, oil swelling increases

3

mol/cm )

1.05 3.77 2.20 2.70 3.20 1.42

Z

f(co)

0.886 0.634 0.821 0.734 0.689 0.848

1.05 1.24 1.12 1.16 1.19 1.07

-2

∆ng (10

1.08 5.00 2.14 2.10 4.51 1.42

mol)

analytical

numerical

6.50 7.58 5.60 6.30 7.58 6.80

6.54 8.01 5.98 6.40 7.95 6.90

the total gas mass transfer rate into the oil phase, as shown in Figure 7, which tends to cause pressure decay faster in the measurement. On the other hand, the volume of oil swelled out of the core decreases the gas-phase volume, which tends to make the pressure decay slower. For instance, for the measurement EDCO2/O-2 with a swelling factor of 1.24, the total amount of CO2 entered into the oil phase at the end of the measurement is increased by 20.3% compared with the nonswelling base case.

9316

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009

Such increase in the amount of CO2 entered into the oil phase would cause the same percentage of increase of total pressure drop (∆P) in the measurement if the gas-phase volume is assumed to be constant. However, a volume of about 2.25 cm3 oil was estimated to be swelled out of the core sample at the end of the measurement. The contribution of swelled out oil to gas-phase pressure change can be approximated by the term (∆VPo)/(V - ∆V) in eq 25. The calculation shows that the swelled-out oil tends to reduce the total pressure drop by 22.6% compared with the nonswelling base case. Therefore, these two effects of oil swelling on pressure decay cancel each other to a large extent for a measurement. The net effect of swelling is to make the pressure drop slightly smaller for a swelling case than that of the nonswelling base case. As a result, the effective diffusion coefficient calculated with the analytical model (nonswelling model) should be slightly smaller than that of the calculated with the numerical solution (swelling model). The above analysis was confirmed by the measured effective diffusion coefficients shown in Table 2. The values of Deff determined by analytical solution are slightly lower than those determined by numerical solution (with the maximum relative difference of about 5%). It should be mentioned, however, that the values of Deff determined by the numerical solution (swelling model) can be different from those determined by the analytical solution (nonswelling model) depending largely on the experimental conditions and different gas-liquid systems. The analytical model can always be used as the first step to quickly estimate a value of Deff for a given measurement, whereas the numerical solution should be used to determine the actual value of Deff and to examine the acceptability of the analytical solution. Applications. The proposed design in this study provides a method measuring the CO2 effective diffusion coefficient in liquid-saturated porous media, which is an important input parameter for large-scale modeling of CO2-EOR or the CO2 sequestration process. For most gas-liquid systems, the conventional PVT method, which has been widely used to measure gas diffusion coefficient, has difficulties being applied to most CO2-liquid systems due to the interference of density-induced natural convection. The proposed model eliminates the effect of natural convection, and the effective diffusion coefficient of CO2 in liquid-saturated porous media can be measured. The measured effective diffusion coefficient Deff (eq 3.a) already includes the diffusive tortuosity and thus can be used directly in modeling CO2 transport in liquid-saturated porous media without the necessity to evaluate the diffusive tortuosity factors of the porous media to convert the diffusion coefficient in the liquid phase to that in porous media. In addition, the proposed method can be used to measure the diffusion coefficient of CO2 in a bulk liquid phase if the diffusive tortuosity factor of the porous medium can be measured. In practice, the diffusive tortuosity factors of porous media can be determined by comparing the measured diffusion coefficients of some gases, such as CH4, in a bulk liquid phase and in liquidsaturated porous medium. Since there is no density-induced natural convection issues for most CH4-liquid systems, the measurement of CH4 in liquid phase can be done with the conventional PVT method. The effective diffusion coefficient of CH4 in liquid-saturated porous medium can be done with the proposed method in this study. With the measured effective diffusion coefficient Deff and the diffusive tortuosity factor ε,

the diffusion coefficient of CO2 in bulk liquids can be determined according to eq 3.a. Conclusions and Recommendations An experimental method along with mathematical model and numerical solutions was developed for measuring CO2 effective diffusion coefficient in oil-saturated porous media under highpressure high-temperature conditions. The proposed design ensures the experiment to be less sensitive to the environments, and the resulted pressure decay due to the diffusion process measurable. The effect of swelling-induced convection as CO2 diffuses into the oil phase was considered in the mathematical model. The overall swelling effect on the pressure decay process was tested to be insignificant for the designed system under the experimental conditions in this study. However, the swelling effects may have different impacts for various gas-liquid systems under different experimental conditions. Thus, the derived numerical solution should always be used to examine the applicability of the analytical solution which does not consider the swelling. The proposed experimental method can be readily implemented in most PVT laboratories. It provided a tool for further studies of gas diffusion in liquid-saturated porous media. The effective diffusion coefficients measured through the proposed method already incorporated the characteristics of porous media and thus can be used directly in modeling CO2 transport in reservoirs. The proposed method also provides an alternative to back out the diffusion coefficient of CO2 in the bulk liquid phase from its effective diffusion coefficient in liquid-saturated porous media and the diffusive tortuosity factor of the porous media. Acknowledgment The financial support from the Natural Sciences and Engineering Research Council (NSERC) of Canada and the Petroleum Technology Research Center (PTRC), Regina, Saskatchewan is gratefully acknowledged. Nomenclature c ) concentration of gas in the liquid phase (mol cm-3) ci ) initial concentration in the liquid phase (mol cm-3) co ) concentration of gas in the liquid phase at the gas/liquid interface (mol cm-3) cr ) concentration of gas based on bulk volume of the porous medium (mol cm-3) D ) diffusion coefficient of gas in liquid (cm2 s-1) Deff ) gas effective diffusion coefficient in the porous medium in eq 1.a (cm2 s-1) Deff ′ ) gas effective diffusion coefficient in the porous medium in eq 1.b (cm2 s-1) e ) maximum relative error tolerance in iteration calculations f ) swelling factor g ) gravitational acceleration (980 cm s-2) h ) height of the porous column (cm) n ) total mass flux (mol cm-2 s-1) P ) pressure (MPa or kPa) Po ) initial pressure of the gas phase in the diffusion measurement (kPa) q ) quantity of gas entered into a porous column at time t (mol) q∞ ) quantity of gas entered into a porous column at infinite time (mol) r ) radial distance of a porous column (cm)

Ind. Eng. Chem. Res., Vol. 48, No. 20, 2009 ro ) radius of a porous column (cm) jr, jz, uj, Vj ) dimensionless variables corresponding to t, r, z, u, and V, respectively R ) universal gas constant (8314 kPa cm3 mol-1 K-1) t ) time (s or min) T ) temperature (K) u ) velocity components in r directions in a cylindrical system (cm s-1) V ) volume of the gas phase in the diffusion cell (cm3) Z ) gas compressibility factor Greek Symbols ∆ng ) quantity of gas lost in the gas phase caused by diffusion (mol) ∆P ) pressure drop in the gas phase (kPa) ∆V ) volume of oil swelled out of the porous column in the diffusion measurement (cm3) ∆Ve ) volume increase of the oil phase within a circular element of the column (cm3) ε ) diffusive tortuosity of porous rock φ ) porosity of porous medium τ ) dimensionless time

Literature Cited (1) Chukwuma, F. O. Mass Transfer of Gaseous Carbon Dioxide into a Quiescent Liquid Hydrocarbon. Ph.D. Thesis, the University of Tulsa, Tulsa, OK, 1983. (2) Renner, T. A. Measurement and Correlation of Diffusion Coefficients for CO2 and Rich Gas Applications. SPE ReserVoir Eng. 1988, (May), 517– 523. (3) Tan, K. K.; Thorpe, R. B. Gas Diffusion into Viscous and NonNewtonian Liquids. Chem. Eng. Sci. 1992, 47, 3565–3572. (4) Li, Z.; Dong, M.; Shirif, E. Transient Natural Convection Induced by Gas Diffusion in Liquid-Saturated Vertical Porous Columns. Ind. Eng. Chem. Res. 2006, 45 (9), 3311–3319. (5) Li, Z.; Dong, M.; Li, S.; Dai, L. A New Method for Determination of Gas Effective Diffusion Coefficient in Brine-Saturated Porous Rocks under High Pressures. J. Porous Media 2006, 9 (5), 445–461. (6) Pomeroy, R. D.; Lacey, W. N.; Scudder, N. F.; Stapp, F. P. Rate of Solution of Methane in Quiescent Liquid Hydrocarbons. Ind. Eng. Chem. 1933, 25, 1014–1019. (7) Grogan, A. T.; Pinczewski, V. W.; Ruskauff, G. J.; Orr, F. M., Jr. Diffusion of CO2 at Reservoir Conditions: Models and Measurement. SPE ReserVoir Eng. 1988, (Feb), 93–102.

9317

(8) Riazi, M. R. A New Method for Experimental Measurement of Diffusion Coefficients in Reservoir Fluids. J. Pet. Sci. Eng. 1996, 14, 235– 250. (9) Zhang, Y. P.; Hyndman, C. L.; Maini, B. B. Measurement of Gas Diffusivity in Heavy Oils. J. Pet. Sci. Eng. 2000, 25, 37–47. (10) Upreti, S. R.; Mehrotra, A. K. Experimental Measurement of Gas Diffusivity in Bitumen: Results for Carbon Dioxide. Ind. Eng. Chem. Res. 2000, 39, 1080–1087. (11) Bear, J. Dynamics of Fluids in Porous Media; American Elsevier Publishing Company: New York, 1972. (12) Leythaeuser, D.; Schaefer, R. G.; Yukler, A. Diffusion of Light Hydrocarbons through Near Surface Rocks. Nature 1980, 284, 522–525. (13) Krooss, B. M.; Leythaeuser, D.; Schaefer, R. G. The Quantification of Diffusive Hydrocarbon Losses through Cap Rocks of Natural Gas ReservoirssA Reevaluation: Reply. AAPG Bull. 1992, 76, 1842–1846. (14) Cussler, E. L. Diffusion: Mass Transfer in Fluid Systems, 2nd ed.; Cambridge University Press: New York, 1997. (15) Krooss, B. M.; Schaefer, R. G. Experimental Measurement of the Diffusion Parameters of Light Hydrocarbon in Water-Saturated Sedimentary RockssI. A New Experimental Procedure. Org. Geochem. 1987, 11, 193– 199. (16) Krooss, B. M.; Leythaeuser, D. Experimental Measurement of the Diffusion Parameters of Light Hydrocarbon in Water-Saturated Sedimentary RockssII. Results and Geochemical Significance. Org. Geochem. 1988, 12, 91–108. (17) Simon, R.; Graue, D. J. Generalized Correlations for Predicting Solubility, Swelling and Viscosity Behavior of CO2-Crude Oil Systems. J. Pet. Technol. 1965, (Jan), 102–106. (18) CoreLab, Operations Manual, Profile Permeameter PDPK-400; Core Laboratories Instruments: Houston, TX, 1996. (19) Bolotnikov, M. F.; Neruchev, Y. A.; Melikhov, Y. F.; Verveyko, V. N.; Verveyko, M. V. Temperature Dependence of the Speed of Sound, Densities, and Isentropic Compressibilities of Hexane + Hexadecane in the Range of (293.15 to 373.15) K. J. Chem. Eng. Data 2005, 50 (3), 1095– 1098. (20) Viswanath, D. S.; Natarajan, G. Data Book on the Viscosity of Liquids; Hemisphere Publishing Corporation: New York, 1989; p 380. (21) WinProp-Phase Property Program 2004 User’s Guide; Computer Modeling Group (CMG) Ltd.: Calgary, Alberta, Canada, 2004. (22) Li, Z. Study of Gas Diffusion in Liquid-Saturated Porous Media for Oil Recovery and CO2 Sequestration. Ph.D. Thesis, University of Regina, Regina, Saskatchewan, Canada, 2006. (23) Reamer, H. H.; Opfell, J. B.; Sage, B. H. Diffusion Coefficients in Hydrocarbon Systems Methane-Decane-Methane in Liquid Phase. Ind. Eng. Chem. 1956, 48, 275–282.

ReceiVed for reView January 27, 2009 ReVised manuscript receiVed July 14, 2009 Accepted August 4, 2009 IE900145C