A Model for Direct Estimation of Wetting Phase Relative Permeabilities

Article pubs.acs.org/IECR

A Model for Direct Estimation of Wetting Phase Relative Permeabilities Using a Multistep Drainage Process Shengdong Wang,† Mingzhe Dong,*,† and Jun Yao‡ †

Department of Chemical and Petroleum Engineering, University of Calgary, Calgary, AB, T2N 1N4, Canada College of Petroleum Engineering, China University of Petroleum (Huadong), Qingdao, Shandong, China

‡

ABSTRACT: A new analytical model has been developed for determining the wetting phase relative permeabilities in a multiple gas−liquid drainage process. The multistep drainage process consists of a series of drainage processes with a tiered pressuredifference profile. With this model, the wetting phase relative permeabilities can be readily obtained through regression of the wetting phase recovery history. The assumptions used in deriving the analytical model are examined by numerical simulations. The agreement between the analytical and the numerical results indicates that these assumptions are reasonable and valid. Numerical simulations and laboratory results demonstrate that this model is effective for direct estimation of the wetting phase relative permeabilities in the multistep drainage process.

1. INTRODUCTION CO2 sequestration1 in geological formations provides an appealing option to reduce CO2 concentration in the atmosphere and thereby mitigates the impacts of the greenhouse gas. The sequestration injection of CO2 gas into the geological formation is a process where the nonwetting gas phase displaces the wetting liquid phase(s). Near the end of this process, the capillary pressure of the two phases becomes comparatively high, and the relative permeability of the wetting phase becomes very low. A clear understanding concerning the capillary pressure and the relative permeabilities at low wetting phase saturation can assist to maximize CO2 sequestration capacity in geological formations. The porous plate method,2 which is referred to as the multistep drainage process in this paper, is one of the conventional methods used to measure the drainage-type capillary pressure versus liquid saturation relationship of a porous medium. In this process, the wetting phase recovery history is determined by the capillary pressure and the relative permeabilities, simultaneously. The capillary pressure versus saturation curve is determined by displacing the wetting fluid point by point, at each of which, the pressure equilibrium is established. The challenge is how to determine the relative permeabilities in this process. One-step drainage experiments were conducted using an air−wetting soil column system, combined with a parameter optimization algorithm to simultaneously estimate the capillary pressures and the relative permeabilities.3 A gradient-based automatic history matching method was used by Jennings et al.4,5 with a one-dimensional (1-D) numerical simulator to obtain the wetting phase relative permeabilities in a core sample. Chen et al.6 considered oil−water flow in a multistep drainage process and obtained the relative permeabilities to both oil and water simultaneously by numerical simulation and history matching. However, all the above inverse methods using history matching require an accurate simulator and an effective history matching algorithm. Moreover, the main concern associated with the inverse function estimation method is that © 2012 American Chemical Society

the result of history matching is not unique. Although said concern was more or less mitigated by simultaneously measuring internal capillary pressure7 and extending the onestep drainage process to a multistep drainage process,8 the uniqueness of the matched results still requires systemic examination. Eching et al.9 and Liu et al.10 used the wetting phase production data at the time periods immediately following an increase in nonwetting phase pressure to estimate permeability functions. However, the flow rate changes rapidly at the beginning of a drainage process. If simple analytical equations can be developed to match the entire recovery history, not only will there be no further concern about the uniqueness of the results but the relative permeabilities of the wetting phase can be directly estimated. As early as the 1950s, some models for calculating the relative permeabilities using a capillary pressure curve were developed for imbibition-type displacement based on the assumption that the wetting phase and the nonwetting phase flow independently in the small pores and the large pores, respectively.11,12 This assumption may not always be valid because, at low wetting phase saturation, the wetting phase stays at the corner of the large pore structure and starts contributing to the flow rate. For imbibitions, analytical or empirical equations modeling the volume variation of an imbibed wetting phase into a nonwetting-phase-saturated core sample have been presented since 1950s. Handy’s equation13 is usually used to describe the variations of the volume of the imbibed wetting phase as a function of the imbibition time: Vw 2 = A2

2PcKK rwφSwf t μw

(1)

where A is the cross-section area of the core sample, Vw is the volume of water imbibed into the core sample, and ϕ and μw Received: Revised: Accepted: Published: 15472

March 9, 2012 October 6, 2012 November 1, 2012 November 1, 2012 dx.doi.org/10.1021/ie300638t | Ind. Eng. Chem. Res. 2012, 51, 15472−15483

Industrial & Engineering Chemistry Research

Article

are the porosity of the core sample and the viscosity of water, respectively. t is the imbibition time; Swf is the water saturation behind the imbibition front. Krw and Pc are the relative permeability of water and the capillary pressure at Swf. It should be noted that the imbibition volume of the wetting phase cannot go to infinity when the time goes to infinity.14 Other models describing the process of spontaneous water imbibition into water-wet rocks have been suggested. The first model was presented by Aronofskyn et al.15 The variation of oil recovery as a function of time is governed by R = R 0(1 − e−σt )

(2)

where R is the oil recovery, R0 is the limit toward which the recovery converges, σ is a constant giving the rate of convergence, and t is the imbibition time. This model was based on the following two empirical assumptions: (1) cumulative oil recovery from a small volume of core converges to a finite limit, and (2) the convergence rate and the convergence limit do not change during the imbibition process. On the basis of Aronofsky’s model, several modifications or derivations of the parameter λ were later reported.16−19 Recently, Li and Horne20 developed the following model to describe the relationship between the water recovery history and imbibition time in a free imbibition process where water imbibed into a gas-saturated core sample. (1 − R *)e R * = e−tD

Figure 1. Cylindrical core disk and membrane configuration.

the wetting phase saturation starts from the previous equilibrium state, the gas phase pressure at the outlet increases to a new level, and the wetting phase recovery history is recorded. Since the saturation variation is small within a drainage step and the porous disk is relatively thin, the following assumptions are made in developing the analytical model: (1) The porous medium is homogeneous. (2) The effect of gravity is negligible compared to the capillary forces, due to the thinness of the core disk. (3) The flow of fluids in the core disk is governed by Darcy’s law. (4) At the end of each drainage step, the saturation of the wetting phase reaches an equilibrium state. (5) The pressure drop of the gas phase along the drainage direction within the porous sample is negligible, due to high gas mobility. (6) Wetting phase relative permeability within a drainage step is considered a constant. (7) Within a drainage step, the capillary pressure is linearly dependent on the wetting phase saturation. Assumptions 1−4 are reasonable for any Darcy’s flow in a porous medium, while rationales of assumptions 5−7 are based on the numerical analysis in the later sections of this paper and will be discussed then. 2.1. Membrane Resistance Negligible. If it is assumed that the hydraulic resistance of the membrane is negligible, an analytical model which describes the correlation of the wetting phase recovery and the drainage time can be derived. As shown in Figure 2, the length of the core disk is L, the thickness of the membrane is Lm, Pin is the nonwetting (gas) phase pressure at the inlet, and Pout is the wetting phase pressure at the outlet. The ith step drainage process starts from the initial equilibrium water saturation Sw(i). The gas phase pressure has a sudden increase from Pin(i) to Pin(i+1) at the beginning of the ith drainage process. The wetting phase saturation correspondingly decreases and finally reaches Sw(i+1) at the end of the ith drainage process. For one-dimension two-phase flow in a homogeneous core disk, combining the material balance equation and the Darcy’s equation gives the following continuity equation for the wetting phase:21

(3)

where R* is the normalized recovery (R* = cRpv). Rpv is the recovery in the units of pore volume, and c is a parameter associated with the ratio of gravity to capillary pressure. However, all of these models are used for a one-step imbibition process. Multistep displacement experiments, either imbibition processes or drainage processes, require a porous plate or a thin membrane to control the pressure for each saturation point. There is no analytical model for the multistep drainage process, especially when considering the resistance of the porous plate or the membrane. This paper presents a new analytical model, which describes the wetting phase recovery history in a multistep drainage process. From the regression of the wetting phase recovery curve at each pressure step, the wetting phase permeabilities at that saturation point can be determined. In order to examine the model and its assumptions, numerical simulations are conducted for one-dimensional two-phase drainage-type displacement. Additionally, numerical simulations are carried out to test and validate the analytical model. Finally, the experimental data reported by Jennings4 are used to further validate this method. Comparisons of the results of the model in this paper and the results presented in Jennings’ work point to the effectiveness and validity of the model proposed in this paper.

2. THEORETICAL DEVELOPMENT As shown in Figure 1, the model for the multistep drainage process is a cylindrical core disk, with its side surface sealed. The outlet of the core disk connects to a thin semipermeable membrane. Within the breakthrough pressure of the membrane for the gas phase (nonwetting phase), only the wetting phase is allowed to flow through the membrane. Gas is injected into the sample through the top surface and the wetting phase is expelled from the lower side. An entire multistep drainage process consists of several single drainage steps. In each step,

∂S ∂ ⎛ KK rw ∂Pw ⎞ ⎜⎜ − ⎟⎟ = −ϕ w ∂t ∂x ⎝ μw ∂x ⎠ 15473

(4)

dx.doi.org/10.1021/ie300638t | Ind. Eng. Chem. Res. 2012, 51, 15472−15483


Article

Figure 2. Conceptual diagram of the analytical model. Figure 3. Capillary pressure considered as a linear function of the wetting phase saturation in one drainage step.

where t is the drainage time; x is the coordinate along the core disk; K and ϕ are the instinct permeability and porosity of the core disk, respectively; μw is the viscosity of the wetting phase; and Krw, Pw, and Sw are the relative permeability, the pressure, and the saturation of the wetting phase in the core disk, respectively, at time t and position x. If λw = (KKrw/μw) is defined as the wetting phase mobility and λw is assumed be a constant within a single drainage step, the following equation can be obtained: 2

λw

∂ Pw ∂x 2

=ϕ

∂Sw ∂t

Pn = Pin

The capillary pressure at x is equal to the pressure difference between the gas phase and the liquid phase. Therefore Pw(x , t ) = Pn − Pc(x , t )

(5)

Pc(x , t ) = Pc(i + 1) +

− Sw

(i + 1)

(Sw(x , t ) − Sw (i + 1))

Substituting eqs 10−13 into eq 14 gives ⎛ ⎞ −ΔPc (Sw(x , t ) − Sw (i + 1))⎟ Pw(x , t ) = Pin − ⎜Pc(i + 1) + ΔSw ⎝ ⎠

(7)

(15)

Substituting eq 15 into eq 5 gives λw

(8)

∂Sw (x = 0, t ) = 0 ∂x

(9)

(10)

ΔSw = Sw (i) − Sw (i + 1)

(11)

(16-1)

(16-2)

Substituting eq 15 into eq 9 and considering

As shown in Figure 3, if the capillary pressure is assumed to vary linearly as saturation changes within a single drainage step, the wetting phase saturation decreases from Sw(i) to Sw(i+1) during one drainage step. The capillary pressure increase and the wetting phase saturation decrease in a drainage step are represented by, respectively ΔPc = Pc(i + 1) − Pc(i)

ΔPc ∂ 2Sw ∂S =ϕ w ΔSw ∂x 2 ∂t

Substituting eq 15 into eq 8 gives

At the outlet of the core disk, if the resistance of the membrane is ignored, the boundary condition is written as Pw(x = L , t ) = Pout

Sw

(i)

(14)

where Pout is the wetting phase pressure at the outlet and Sw(i) the initial equilibrium water saturation at the beginning of this drainage step. As there is no wetting phase flowing into the core disk at the inlet, the boundary conditions for the wetting phase are written as ∂Pw (x = 0, t ) = 0 ∂x

Pc(i) − Pc(i + 1)

(6)

and Sw(x , t = 0) = Sw (i)

(13)

As shown in Figure 3, the capillary pressure at any time t and any location x has the following relationship to the wetting phase saturation:

The initial conditions for the pressure and saturation of the wetting phase are

Pw(x , t = 0) = Pout

(12)

Pc(i + 1) = Pin − Pout

the first boundary condition changes to Sw(x = L , t ) = Sw(i + 1)

(16-3)

Similarly, substituting eq 15 into eq 10, the initial condition changes to Sw(x , t = 0) = Sw(i)

If the resistance to gas phase flow is negligible, the pressure in the gas is constant, that is

(16-4)

Equation 16-1 to eq 16-4 can be normalized by defining the following dimensionless variables: 15474


Industrial & Engineering Chemistry Research SwD =

Sw − Sw (i + 1) ΔSw

(17)

x L

xD =

Article

(18)

and tD =

λ w ΔPc ϕL2ΔSw

λ12 =

⎛ 3π ⎞2 ⎜ ⎟ ≈ 22.21 ⎝ 2 ⎠

λ22 =

⎛ 5π ⎞2 ⎜ ⎟ ≈ 61.69 ⎝ 2 ⎠

The solution of eq 22 can be approximated by taking the first term

t (19)

SwD ̅ ≈

Substituting the following expressions for water saturation, distance, and time into eq 16,

2 −λ 0 2 t D e = λ02

x = x DL ln SwD ̅ ≈−

ϕL2ΔSw tD λ w ΔPc

ΔPc ∂ 2(SwDΔSw + Sw (i + 1)) ∂(SwDΔSw + Sw (i + 1)) = ϕ φL2ΔSw ΔSw ∂(x DL)2 t ∂ λ w ΔPc D

⎛ π2 ⎞ π2 t D − ln⎜ ⎟ 4 ⎝8⎠

or ∂ 2SwD ∂x D2

qw (x = L , t ) = −

∂S = wD ∂t D

(20-1)

Sw(x = L , t ) − Sw (i + 1) Sw (i) − Sw (i + 1)

=0 qwm =

(20-2)

(i + 1) ⎞ ∂SwD ∂ ⎛ S (x = 0, t ) − Sw ⎟⎟ = 0 (x D = 0, t D) = x ⎜⎜ w ∂x D ∂L ⎝ Sw (i) − Sw (i + 1) ⎠

(20-3)

Sw(x , t = 0) − Sw

SwD(x D , t D) =

∑ n=0

KK rw ∂Pw K m Pw(x = L , t ) − Pout =− (x = L , t ) μw Lm μw ∂x

1 ∞

∫0 ∑ λ2 n=0

2

n

∞

∑ n=0

(28)

Considering the definition of SwD,

2 −λ n 2 t D e λn 2 (22)

Pw(L , t ) = Pin − Pc(i + 1) + ΔPcSwD(x D = 1, t D)

(29)

Pw(L , t ) = Pout + ΔPcSwD(x D = 1, t D)

(30)

Substituting eq 30 into eq 28,

where ⎛π ⎞ ⎜ ⎟ ≈ 2.47 ⎝2⎠ 2

λ02 =

(27)

⎛ (S (L , t ) − Sw (i + 1)) ⎞ ⎟⎟ Pw(L , t ) = Pin − ⎜⎜Pc(i + 1) − ΔPc w ΔSw ⎝ ⎠

(21)

cos λnx D e−λn tD dx D =

(26)

Equation 28 is further simplified by applying eq 15 at the boundary. When x = L, eq 15 can be rearranged as

2 2 cos λnx D e−λn tD , where λn λn

The average saturation in the core disk can be expressed as SwD ̅ =

K m Pw(x = L , t ) − Pout μw Lm

as

(20-4)

2n + 1 = π , n = 0, 1, .... 2

(25)

The boundary condition at the interface is, therefore, written

=1

In order to obtain an analytical solution to eqs 20-1 to 20-4, the separation of variables method is applied and the solution is expressed as ∞

KK rw ∂Pw (x = L , t ) μw ∂x

qw (x = L , t ) = qwm

(i + 1)

Sw (i) − Sw (i + 1)

(24)

where Pw(x = L,t) is the wetting phase pressure at the interface of the core disk and the membrane. qwm is the wetting phase flow rate in the membrane. From material balance

The initial condition of eq 16-4 becomes SwD(x D , t D = 0) =

(23)

If the membrane is considered as a homogeneous porous medium, the wetting phase flow rate in the membrane is expressed as

Boundary conditions of eq 16-2 and eq 16-3 are SwD(x D = 1, t D) =

2

( π2 )

2

2.2. Membrane Resistance Considered. If the resistance of the membrane at the outlet boundary is not negligible compared to that of the core disk, the boundary condition at the lower end of the core disk (the interface between the core disk and the membrane) will change. At this interface, the wetting phase flow rate in the core disk is governed by Darcy’s Law22 as

Equation 16-1 then changes to λw

π

e −( 2 ) t D

Taking the logarithm of the both sides of the above equation gives

Sw = SwDΔSw + Sw (i + 1)

t=

2

KK rw ∂S Km ΔPcSwD(x D = 1, t D) = − ΔPc wD (x D = 1, t D) Lm L ∂x D 15475


Industrial & Engineering Chemistry Research ⎛ K ⎞ ∂S ⎜ ⎟ SwD(x D = 1, t D) = −⎜ KL ⎟K rw wD (x D = 1, t D) m ∂x D ⎝ Lm ⎠

Article

(31)

If a dimensionless parameter α, which essentially means the conductivity ratio of the core disk to the membrane, is defined as: ⎛ K ⎞ ⎜ ⎟ α = ⎜ KL ⎟ m ⎝ Lm ⎠

(32)

The new boundary condition is written as SwD(x D = 1, t D) = −αK rw

∂SwD (x D = 1, t D) ∂x D

(33)

The partial differential equation and its boundaries' conditions are summarized as ∂ 2SwD ∂x D

2

=

∂SwD ∂t D

Figure 4. Graphic demonstration for the solutions for the eigenvalues λ in eq 38: (a) the first three solutions, (b) the first solution in domain [π/4,π/2].

(34-1)

SwD(x D = 1, t D) = −αK rw

∂SwD (x D = 1, t D) ∂x D

(34-2)

∂SwD (x D = 0, t D) = 0 ∂x D

(34-3)

SwD(x D , t D = 0) = 1

(34-4)

Similar to eq 22, the solution of eq 40 can also be approximated by taking the first term. However, there is no direct analytical solution for the first eigenvalue λ0. Here, an approximation method to calculate λ0 is provided. Noticing that the conductivity of the core disk is usually smaller than the membrane and Krw is always less than 1, we can say that αKrw < 1. If αKrw < 1 at the domain of [π/4,π/2], as shown Figure 4b, the function ctgλ0 is approximated as π ctgλ 0 ≈ − λ 0 (41) 2 Substituting eq 41 into eq 38, we get π 1 λ0 ≈ 2 1 + αK rw (42)

The solution for eq 34-1 is in the form of: 2

SwD(x D , t D) = (A cos λx D + B sin λx D) e−λ tD

(35)

Applying boundary condition 34-3 to eq 35, 2

( −Aλ sin(λ · 0) + Bλ cos(λ · 0)) e−λ tD = 0

(36)

B has to be zero to make eq 36 satisfied for every λ. Equation 35 is thus changed to

Equation 40 changes to

2

SwD(x D , t D) = A cos λx D e−λ tD

(37)

SwD ̅ ≈

Applying boundary condition 34-3 to eq 37 gives 2

2

A cos λx D e−λ tD = αK rwλA sin λx D e−λ tD

ctgλ = αK rwλ

(38)

ln SwD ̅ ≈−

Equation 38 is a nonlinear equation, whose solution can be obtained numerically. Figure 4a depicts the solutions of eq 38 schematically, which are the cross-points of the straight line y = αKrwλ and the curves y = ctgλ. The general solution to eq 37 can be expressed as: ∞

∑ n=0

2 2 cos λnx De−λn tD λn

(39)

where λn is the n solution of eq 38. The average wetting phase saturation over the entire core disk is: SwD ̅ =

∫0 ∑ n=0

2 2 cos λnx D e−λn tD dx D = λn

∞

∑ n=0

π2 π2 t − ln D 4(1 + αK rw )2 8(1 + αK rw )2

(44)

2.3. Applications of Analytical Models. Equations 22 and 40 presented in the previous sections are models with and without considering the membrane resistance effect, respectively. By taking the first term as an approximation, both models change to linear equations, eqs 24 and 44. It is also noticed that when αKrw approaches zero, eq 44 converges to eq 24, which means that, if the resistance of the membrane is negligible, these two equations are consistent with each other. For the application of this model, at each drainage step we define another dimensionless time, tDk, as

th

1 ∞

(43)

Taking the logarithm of the both sides of the above equation gives

That is

SwD(x D , t D) =

8(1 + αK rw )2 −π 2 /4(1 + αK rw)2 tD 2 −λ 0 2 t D = e e λ02 π2

2 −λn 2tD e λn 2

t Dk =

⎛ π ⎞2 K ΔPc ⎜ ⎟ t ⎝ 2 ⎠ φL2μ ΔSw w

(45)

The dimensionless wetting phase saturation Sw̅ D is defined as

(40) 15476



Article

(46)

Sw|t = 0 = Swi

(55)

where Vw is the cumulative produced volume of the wetting phase and Vt is the total volume expelled by the nonwetting phase at the end of this drainage step. Equation 24 and eq 44 then become, respectively

Pw|t = 0 = Pout

(56)

SwD ̅ = 1 − Vw /Vt

⎛ π2 ⎞ −ln SwD ̅ ≈ K rwt Dk + ln⎜ ⎟ ⎝8⎠

where Swi is the initial wetting phase saturation and Pout is the wetting phase pressure at the outlet. After the drainage process starts, at the inlet, the boundary conditions are written as

(47)

∂Pw ∂x

and −ln SwD ̅ ≈

K rw (1 + αK rw )2

t Dk

2.4674 + ln 2(1 + αK rw )2

∂Pn ∂x

(52)

(53)

Sw + Sn = 1

(54)

Pout = Pout′ + qw |e2 μw ηm

(61)

ΔPm μw qw

(62)

4. VALIDATION OF ASSUMPTIONS A program was developed using Visual C++ 6.0 to solve the above numerical model. The equations were solved using the IMPES (Implicit Pressure Explicit Saturation) method and a one-dimensional model of 27 grid blocks. The simulation was run at a desktop computer with a 2.0 GHz CPU and 4 GB of memory for a couple of minutes. To benchmark the results of the C++ program, saturation and pressure profiles in the core disk before the nonwetting phase front reaches the membrane were examined by using a commercial simulation package (CMG-IMEX). In order to confirm the validity of the analytical model and the rationality of its assumptions, numerical simulations for a single drainage step were carried out. At the beginning of this drainage step, the system was assumed to have reached an equilibrium state, and the average saturation of the wetting phase was 0.40. To start the drainage process, the nonwetting phase pressure at the outlet was increased from 0.20 to 0.32 atm. The wetting phase saturation was, consequently, reduced from 0.40 to 0.36, which was in the middle section of the wetting phase saturation of the capillary pressure curve. The recovery history of the drainage process was recorded. As depicted in Figure 5, the capillary pressure curve and the relative permeability curves were calculated using Corey’s equations.23 The coefficients used in Corey’s equations are listed in Table 1. Other properties of the core disk and the fluids adopted in the numerical simulations are listed in Table 2.

To solve the above equations, the two following equations are required: Pc = Pn − Pw

(60)

where ΔPm is the viscous pressure drop to yield a flow rate of qw in the membrane. Before the multistep drainage experiment, this conductivity can be measured.

3. NUMERICAL MODELING In order to validate the assumptions of the above analytical model, and carry out numerical simulations, a numerical model was developed. In this model, the drainage process was considered a one-dimension two-phase flow model, considering the wetting phase, the nonwetting phase, and the hydraulic resistance of the membrane:

⎛ KK ⎞ ϕ∂(Sw ) ∇·⎜⎜ rw ∇Pw ⎟⎟ + qw = μ ∂t ⎝ w ⎠

(59)

Pw|e2 = Pout′

ηm =

(50)

(51)

=0 e2

where e1 and e2 represent the inlet and outlet boundaries, respectively; Pin is the nonwetting phase pressure at the inlet; Pout′ is the wetting phase pressure at the interface between the membrane and the core disk; and ηm is the hydraulic conductivity, which is defined as

(49)

⎛ KK ⎞ ϕ∂(Sn) ∇·⎜⎜ rn ∇Pn⎟⎟ + qn = μ ∂t ⎝ n ⎠

(58)

At the outlet, the boundary conditions are given as

If αk is small enough, k can be directly used as an estimation of the wetting phase relative permeability. Otherwise, the following equation is required to revise the permeability to the wetting phase: K rw = k(1 + αk)2

(57)

Pn|e1 = Pin

(48)

The slopes of the linear regression are dominated by the permeability of the wetting phase at this step. In order to estimate the relative permeability, tDk and Sw̅ D are plotted as a semilog graph. The slop of the curve k is determined by linear regression of the experiment points. To estimate the effect of the membrane, a dimensionless parameter is calculated: ⎛ K ⎞ ⎜ ⎟ αk = ⎜ KL ⎟k m ⎝ Lm ⎠

=0 e1

where K is the permeability of the core disk; Krn and Krw are relative permeabilities for the nonwetting and the wetting phase, respectively; μn and μw are viscosities of the nonwetting and the wetting phase, respectively; ϕ is the porosity of the core disk; Sn and Sw are saturations of the nonwetting and wetting phase, respectively; and Pn and Pw are pressures of the nonwetting and wetting phase, respectively. For one drainage step, when t = 0, the initial wetting phase saturation in the core disk is Swi and wetting phase pressure is equal to the pressure at the outlet, Pout. Therefore, the initial conditions are written as 15477



Article

Figure 5. The capillary pressure curve and the relative permeability curves calculated using Corey’s equation.

Table 1. Coefficients Used in Corey’s Equation to Calculate the Relative Permeability Curves coefficients used in Corey’s equation

value

Snr = nonwetting phase residual saturation Swr = wetting phase residual saturation Nn = exponent of nonwetting phase Nw = exponent of wetting phase Krnmax = nonwetting phase maximum relative permeability Krwmax = wetting phase maximum relative permeability

0.2 0.15 2.0 2.0 0.9 0.4

Table 2. Properties of the Core Disk and the Fluids Used in Numerical Simulation parameter

value

core disk length core disk cross-sectional area core disk porosity core disk permeability wetting phase viscosity nonwetting phase viscosity membrane conductivity

1.00 cm 16.00 cm2 20% 100 md 1.00 cP 0.02 cP 60 d/cm

Figure 6. Comparisons of the wetting phase recovery history of the simulation model, the full analytical model and the one term approximation model.

in developing the analytical model that the nonwetting phase resistance is negligible. Furthermore, the results of the one-term approximation model show a little difference at the very beginning of the process, but at later stages, the approximation model shows the same results as the numerical solution and the full analytical model. The one term approximation model is a reasonable approximation of the full analytical model. Moreover, when both the cumulative production volume of the wetting phase and time in Figure 6a are normalized and replotted as −ln(SwD) versus tD graph in Figure 6b, a perfect linear relationship between these two variables shows up. As we discussed in the above section, the slope of the straight line is determined by the wetting phase relative permeability. 4.2. Membrane Resistance. In the conventional porous plate method, the porous plate is used to prevent the nonwetting phase from breaking through and being produced from the outlet. This porous plate usually has high gas entry pressure and low permeability. Although Jennings4 replaced the porous plate with a plastic membrane and reduced the experimental time to some extent, the resistance of the membrane may still not be negligible. Therefore, so as to investigate the impacts of the membrane, the membrane hydraulic conductivity was taken as a reasonable number (60 md/cm) in the numerical simulations. The cumulative production volume of the wetting phase obtained from numerical simulation is plotted as circles in Figure 7a. For

4.1. Gas Phase Flow Resistance. In order to investigate the impacts of the nonwetting phase resistance, the viscosity of air under standard conditions was used as the nonwetting phase viscosity. The viscosity of the nonwetting phase was then taken as 0.02 cP, and its relative permeabilities are shown in Figure 5. The wetting phase permeability was taken as a constant, while the membrane hydraulic conductivity was assigned a huge number to make the resistance of the membrane negligible. The cumulative production volume of the wetting phase obtained from simulation is plotted as circles in Figure 6a. For comparisons, the calculated cumulative production volume of the wetting phase using the full analytical solution, eq 22, and the one term approximation solution, eq 24, were also plotted. These results were shown as a dash line and solid line, respectively. It can be seen that the cumulative production volume of the wetting phase obtained from numerical simulation considering the nonwetting phase resistance and the curve using the analytical model without considering the nonwetting phase resistance are almost identical. This implies that the wetting phase recovery history is not sensitive to the nonwetting phase if its mobility is large enough and thus confirms the assumption 15478



Article

the wetting phase and drainage time are normalized and plotted as −ln(SwD) versus tD graph in Figure 7b. It can be seen from Figure 7a that, if the resistance of the membrane is neglected, the wetting phase is recovered faster and the slope of the straight line in Figure 7b is bigger. It means that without considering the resistance of the membrane, using eq 22, the relative permeability of the wetting phase is underestimated from the wetting phase recovery history. The estimated permeability combines the resistance of the core disk and the membrane together. However, when the membrane resistance is considered, using eq 40, the results of the analytical model are consistent with the numerical model. In the case when the resistance of the membrane is not negligible, the model considering the membrane should be used to revise the results. 4.3. Flow Functions. The real capillary pressure is not linearly dependent on the wetting phase saturation, and the real relative permeabilities of the wetting phase also vary at different saturations. In the numerical simulations, the relative permeabilities of the wetting phase and the capillary pressure were calculated using monotone cubic spline interpolation.24 The monotone cubic spline interpolation can keep the monotonicity and smoothness of the flow functions, simultaneously. Four simulations were run for four single drainage steps, and the corresponding saturation ranges were 0.45−0.40, 0.40−0.35, 0.35−0.30, and 0.30−0.25, respectively. For the analytical model, the arithmetic mean of the wetting phase permeabilities was used as the constant relative permeabilities, Krw. The capillary pressure gradient was calculated by the ratio of the capillary pressure change to the saturation change in the corresponding drainage step. The results from simulations and from the analytical model are shown in Figure 8, which depicts the relationship of the logarithm of the dimensionless average saturation versus the dimensionless time. It can be seen that both the simulation results and the analytical results show a straight line relationship, but the curves for the analytical model do not match with the curves of the simulation model. This means the slopes of these straight lines are not good approximations of the

Figure 7. Comparison of the wetting phase recovery history of the simulation model with the membrane resistance, the one term approximation model and the complementary model.

comparisons, the cumulative production volume of the wetting phase calculated using the analytical solution neglecting the membrane, eq 22, and the analytical solution considering the membrane, eq 40, are plotted together as a dash line and solid line, respectively. Both the cumulative production volume of

Figure 8. The linear relationship between the dimensionless time and logarithm of the dimensionless average saturation for the analytical model and numerical models. From the left to the right, the curves are wetting phase recovery for the drainage steps with the wetting phase saturation ranges for 0.45−0.4, 0.4−0.35, 0.35−0.3, and 0.3−0.25. 15479



Article

Figure 9. Linear relationship between the dimensionless time and logarithm of the dimensionless average saturation for the analytical model and numerical models. From the left to the right, the curves are for the drainage steps with a wetting phase saturation of 0.45−0.4, 0.4−0.35, 0.35−0.3, and 0.3−0.25.

Figure 10. Comparison of the results estimated directly using the analytical model and the original simulation values without the membrane resistance. (1) Solid line: the relative permeability curve used in the numerical simulation. (2) Circles: estimated relative permeabilities at each drainage step.

used to generate these results in the simulation are 0.059, 0.038, 0.021, and 0.009. All relative errors are smaller than 0.5%.

arithmetic mean of the wetting phase relative permeabilities in the corresponding drainage step. From Figure 8, it can be noted that the linear relationship shows up at the later stage of a drainage process. Using the average saturation and the capillary gradient mounted on the entire saturation range may not be reasonable. Therefore, the wetting phase permeabilities and the capillary pressure gradient at the ending saturation were used in the analytical model. The results were replotted in Figure 9. The simulation results do not match the analytical results, but the slopes of the curves, directly determined by the relative permeability, are consistent. The fitted results are shown as the linear equations in Figure 9. The slopes of these equations are the estimated relative permeabilities, and the corresponding relative permeabilities

5. COMPUTATIONAL TESTS To further validate the effectiveness of the analytical model, a virtual numerical multistep drainage experiment was carried out using the numerical model developed above. The simulated multistep drainage process consists of 10 drainage steps, and the initial wetting phase saturation is 0.7. The hydraulic conductivity of the membrane is 1.0 Darcy/cm. The maximum αKrw is 0.025. Since αKrw is less than 0.1, the relative permeabilities were estimated using the no-membrane model. The relative permeabilities used as the input in the simulations and the estimated relative permeabilities obtained by analytical model are both plotted in Figure 10. The values in 15480



Article

Figure 11. Comparison of the original assumed values and the results estimated directly using the analytical model with and without membrane resistance revising. (1) Solid line: the relative permeability curve used in the numerical simulation. (2) Circles: estimated relative permeabilities at each drainage step without membrane resistance. (3) Rectangles: estimated relative permeabilities at each drainage step considering membrane resistance.

permeabilities to water were obtained by automatic history matching.

the simulation are plotted as a solid line, while the estimated values are plotted as circles. It can be seen that, using the direct estimation method, the relative permeabilities obtained are consistent with the assumed values. However, as αKrw increases, the no-membrane model without the membrane resistance may underestimate the relative permeabilities. In the second numerical simulation, the conductivity of the membrane was reduced to 0.1 Darcy/cm. The corresponding maximum αKrw is 0.25, which is larger than 0.1. Another 10step drainage process was simulated, and the wetting phase relative permeabilities were calculated using the membrane model. Figure 11 shows that, at a relative permeability range less than 0.35, the values estimated by both the model neglecting the membrane and the model considering the membrane are consistent with the values assumed in the simulation. But, as the relative permeabilities increase, the error caused by the resistance of the membrane is more and more remarkable. However, with the model considering the membrane resistance, the relative permeabilities are almost the same as the input values. Several more simulations were then run, and the results indicated that an estimated αKrw value can be used to determine which model should be used to estimate the relative permeability. If the estimated αKrw is less than 0.05, the model neglecting membrane is reliable enough to estimate the relative permeability. If the estimated αKrw is larger than 0.05, the model considering the membrane is required to estimate the wetting phase relative permeabilities from the wetting phase recovery data.

Figure 12. Capillary pressure measured in Jennings’ experiment.

In order to calculate the relative permeabilities to water, a graph of −ln(SwD) versus t was plotted for each step, as shown in Figure 13. It can be seen that in each step, the curve has a linear part, except for step 2. It is shown that the drainage time for step 2 is shorter than the other three steps; thus the valid period of the model developed in this study has not been reached yet. If the drainage time was longer, the linear relationship should be more remarkable. Using the slopes at each plot and both the model neglecting the membrane and the model considering the membrane, the relative permeabilities to water were calculated. The results are plotted in Figure 14. The consistent results showed with the numerical simulations, at a low saturation range, the three relative curves have the same value, while at high water saturation range, the results considering membrane resistance have a better result. If we assume that Jennings’ simulation results can represent the true values of the relative permeabilities to water, it is indicated that the direct estimation method presented in this paper is valid

6. EXPERIMENTAL VALIDATION In order to further validate the effectiveness of the model, Jennings’ experimental data4 were used to calculate the wetting phase relative permeabilities, using the model developed in this study. The capillary pressure curves measured in Jennings’ experiment are shown in Figure 12. In his study, a four-step drainage experiment was carried out, and the relative 15481



Article

Figure 13. Capillary pressure measured in Jennings’ experiment: (a) step one, (b) step two, (c) step three, (d) step four.

permeability of the wetting phase is the permeability at the ending saturation. (3) A dimensionless parameter was created to determine if the resistance of the membrane is negligible. If it is not, a complementary model is developed to modify the wetting phase permeabilities and get better results. (4) Numerical simulations and laboratory experimental data verification demonstrated that this method is effective for the direct estimation of the wetting phase relative permeability.

■

AUTHOR INFORMATION

Corresponding Author

*Tel.: 1 (403) 210-7642. Fax: 1 (403) 284-4852. E-mail: [email protected].

Figure 14. Comparisons of the results obtained using Jennings’ automatic history matching method and the direct estimation methods presented in this study.

Notes

The authors declare no competing financial interest.

■

and effective, because neither simulation nor automatic history matching are required.

ACKNOWLEDGMENTS The authors acknowledge, with thanks, the support from a research Grant of Natural Sciences and Engineering Research Counsel of Canada.

7. CONCLUSIONS An analytical model to describe the wetting phase recovery history in the multistep drainage process was developed in this study. This model can be used to directly estimate the wetting permeabilities by linear regression of the wetting phase recovery history in each single step. The results of the model are consistent with the results of the one-dimensional numerical simulations. Comparisons of the analytical model and numerical simulation indicate the following: (1) The resistance of the nonwetting phase is negligible when estimating the relative permeability of the wetting phase at low wetting phase saturation. (2) The capillary pressure gradient around the ending saturation should be used, and the estimated relative

■

15482

NOMENCLATURE A = Cross-section area of the core sample c = Gravity to capillary pressure factor i = Drainage step index k = Slope of the fitted straight line K = Instinct permeability of the core disk Km = Instinct permeability of the membrane Krn = Relative permeability to the nonwetting phase Krw = Relative permeability to the wetting phase L = Length of the core disk Lm = Length of the membrane dx.doi.org/10.1021/ie300638t | Ind. Eng. Chem. Res. 2012, 51, 15472−15483


Article

(8) Van Dam, J. C.; Stricker, J. N. M.; Droogers, P. Inverse method to determine soil hydraulic functions from multi-step outflow experiments. Soil Sci. Soc. Am. J. 1994, 58, 647−652. (9) Eching, S. O.; Hopmans, J. W.; Wendroth, O. Unsaturated hydraulic permeability from transient multi-step outflow and soil wetting phase pressure data. Soil Sci. Soc. Am. J. 1994, 58, 687−695. (10) Liu, Y. P.; Hopmans, J. W.; Grismer, M. E.; Chen, J. Y. Direct estimation of air−oil and oil−wetting phase capillary pressure and permeability relations from multi-step outflow experiments. J. Contam. Hydrol. 1998, 32, 223−245. (11) Wyllie, M. R. J.; Spangler, M. B. The application of electrical resistivity measurements to the problem of fluid flow in porous media. Bull. A.A.P.G. 1952, 36 (2), 359−403. (12) Wyllie, M. R. J.; Gardner, G. H. F. The generalized KozenyCarmen equation - its Application to problems of multi-phase flow in porous Media. World Oil 1958, 146, 121−146. (13) Handy, L. L. Determination of effective capillary pressures for porous media from imbibition data. Trans. AIME 1960, 219, 75−82. (14) Li, K.; Horne, R. N. Characterization of spontaneous wetting phase imbibition into gas-saturated rocks. SPEJ, Soc. Pet. Eng. J. 2001, 6 (4), 375. (15) Aronofsky, J. S.; Masse, L.; Natanson, S. G. A model for the mechanism of oil recovery from the porous matrix due to wetting phase invasion in fractured reservoirs. Trans. AIME 1958, 213, 17−19. (16) Mattax, C. C.; Kyte, J. R. Imbibitions oil recovery from fractured water-drive reservoirs. Trans. AIME 1962, 255 (June), 177−184. (17) Lefebvre Du Prey, E. Gravity and capillary effects on imbibition in porous media. Soc. Pet. Eng. J. 1978, 18 (June), 195−206. (18) Reis, J. C. An analysis of oil expulsion mechanisms from matrix blocks during steam injection in the naturally fractured reservoir. InSitu 1992, 16 (1), 43−73. (19) Ma, S.; Morrow, N. R.; Zhang, X. Generalized scaling of spontaneous imbibition data for strongly water-wet systems. J. Pet. Sci. Eng. 1997, 18 (3/4), 165−178. (20) Li, K.; Horne, R. N. Generalized scaling approach for spontaneous imbibition: an analytical model. SPEREE 2006, 251−258. (21) Gottfried, B. S.; Guilinger, W. H.; Snyder, R. W. Numerical solutions of the equations for one-dimensional multi-phase flow in porous media. SPE J. 1966, 6 (1), 62−72. (22) Darcy, H. Les Fontaines Publiques de la Ville de Dijon: Dalmont, Paris, 1856; p 647. (23) Corey, A. T. The interrelation between gas and oil relative permeabilities. Prod. Monthly 1954, 19 (1), 38−41. (24) Fritsch, F. N.; Carlson, R. E. Monotone piecewise cubic interpolation. SIAM J. Num. Anal. 1980, 17, 238−246.

Pc = Capillary pressure between the wetting phase and the nonwetting phase Pin = Constant pressure at the inlet Pn = Pressure of the nonwetting phase Pout = Constant pressure at the outlet Pout′ = Pressure of the wetting phase at the interface between the membrane and the core disk Pw = Pressure of the wetting phase qn = Flow rate of the nonwetting phase in the core disk qw = Flow rate of the wetting phase in the core disk qwm = Flow rate of the wetting phase in the membrane R = Oil recovery in an imbibition process R* = Normalized oil recovery R0 = Maximum oil recovery Rpv = The recovery in the units of pore volume Sn = Saturation of the nonwetting phase Sw = Saturation of the wetting phase SwD = Normalized saturation of the wetting phase Sw̅ D = Normalized average saturation of the wetting phase Swf = Water saturation behind an imbibition front Swi = Initial wetting phase saturation t = Drainage/imbibition time tD = Normalized drainage/imbibition time tDk = Normalized drainage time without the relative permeability to the wetting phase Vt = Total cumulative produced volume the wetting phase in one single drainage step Vw = Cumulative volume the wetting phase x = Coordinate along the core disk xD = Normalized coordinate along the core disk Greek Symbols

α = The conductivity ratio of the core disk to the membrane ηm = Hydraulic conductivity of the membrane λw = Mobility of the wetting phase λn = Mobility of the nonwetting phase μn = Viscosity of the nonwetting phase μw = Viscosity of the wetting phase ϕ = Porosity of the core disk

■

REFERENCES

(1) Bachu, S. Sequestration of CO2 in geological media in response of climate change: road map for site selection using the transform of the geological space into the CO2 phase space. Energy Convers. Manage. 2002, 43, 87−102. (2) Bruce, W. A.; Welge, H. J. Restored-state method for determination of oil in place and connate wetting phase. Oil Gas J. 1947, 46, 223−238. (3) Parker, J. C.; Kool, J. B.; Van Genuchten, M. T. Determining soil hydraulic properties from one-step outflow experiments by parameter estimation: II. Experimental studies. Soil Sci. Soc. Am. J. 1985, 49, 1354−1359. (4) Jennings, J. W. A method for measuring capillary pressure and relative permeabilities based upon a transient analysis of a rapid restored-state experiment. PhD dissertation 1983, Texas A&M University, College Station, TX. (5) Jennings, J. W.; McGregor, D. S.; Morse, R. A. Simultaneous Determination of Capillary Pressure and Relative Permeability by Automatic History Matching. SPE Form. Eval. 1988, 3 (2), 322−328. (6) Chen, J. Y.; Hopmans, J. W.; Grismer, M. E. Parameter estimation of two-fluid capillary pressure−saturation and permeability functions. Adv. Water Res. 1999, 22 (5), 479−493. (7) Eching, S. O.; Hopmans, J. W. Optimization of hydraulic functions from transient outflow and soil wetting phase pressure data. Soil Sci. Soc. Am. J. 1993, 57, 1167−1175. 15483


A Model for Direct Estimation of Wetting Phase Relative Permeabilities

Recommend Documents