Graphical estimation of interface mass-transfer coefficient (k) using pressure-decay data

Chemical Engineering Science 146 (2016) 35–43

Contents lists available at ScienceDirect

Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces

Graphical estimation of interface mass-transfer coefficient (k) using pressure-decay data F.J. Pacheco Roman, S.H. Hejazi n Department of Chemical and Petroleum Engineering, University of Calgary, AB, Canada T2N 1N4

H I G H L I G H T S


An approximate analytical solution for gas-heavy oil diffusion problem is derived using a time dependent non equilibrium boundary condition.
A graphical method is proposed to estimate the interface mass-transfer coefficient.
The validity of the proposed method is verified with synthetic and experimental data.

art ic l e i nf o

a b s t r a c t

Article history: Received 18 October 2015 Received in revised form 16 January 2016 Accepted 25 January 2016 Available online 11 February 2016

A new graphical procedure is proposed to estimate the interface resistivity in terms of interface masstransfer coefficient ðkÞ in gas–oil pressure decay experiments. We provide a time-domain closed-form analytical solution by the use of Integral Method (IM) to the equations governing the gas-phase-pressure decay with the improved time-dependent gas–oil interface boundary condition. A graphical technique, derived from this analytical solution, is proposed to extract the interface mass-transfer coefficient and Henry's constant. Application of this method to synthetic and experimental pressure-decay data shows that the estimated parameters are in agreement with those reported in the literature. The proposed procedure, coupled with a graphical technique to estimate the diffusion coefficient, is an asset for a simple, quick, and reliable estimation of interface mass-transfer coefficient and Henry’s constant. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Gas-heavy oil diffusion process Pressure-decay technique Mathematical modeling Interface mass transfer coefficient

1. Introduction The determination of mass transfer parameters, i.e. diffusion coefficient ðDÞ, interface mass-transfer coefficient ðkÞ and solubility ðH ij Þ of gases in liquids, is relevant for many science and engineering applications. In petroleum engineering, these parameters are used to model and design gas injection processes (Riazi et al., 1994). Thus, the efforts to estimate these properties are extensive. The diffusion coefficient ðDÞ can be inferred experimentally by direct or indirect techniques (Sheikha et al., 2005). The intrusive, expensive, and time-consuming nature of direct methods has been the driving force in the use of indirect methods such as the pressure-decay technique (PDT). Riazi (1996) proposed the PDT and subsequently Zhang et al. (2000) adapted the technique for gas/heavy-oil systems. Basically, the PDT utilizes a pressure/ volume/temperature (PVT) cell that contains a known mass of oil and gas at a constant temperature. There, the gas-phase pressure n

Corresponding author. Tel.: þ 1 403 389 4697. E-mail address: [email protected] (S.H. Hejazi).

http://dx.doi.org/10.1016/j.ces.2016.01.049 0009-2509/& 2016 Elsevier Ltd. All rights reserved.

declines as the gas diffuses into the heavy oil and the mass transfer parameters are inferred from the recorded changes in pressure. Although the PDT was originally conceived to estimate D, posterior works showed that the interface mass-transfer coeffi
 cient ðkÞ and Henry's constant H ij could be determined by considering various boundary conditions in the modeling process. For instance, assuming the concentration of gas in heavy oil proportional to the gas pressure and resistance to mass transfer at the gas/heavy-oil interface leads to a solution from which H ij and k can be estimated, respectively. As classified by Tharanivasan et al. (2004), most models utilized an equilibrium or quasi-equilibrium boundary condition at the gas/ heavy-oil interface, while few considered a non-equilibrium one. The effect of the three boundary conditions on the diffusion-coefficient estimation was researched by Tharanivasan et al. (2004, 2006) who suggested their selection on the basis of the particular gas being used. The equilibrium condition, mathematically expressed in Eq. (1), considers that at all times the interface is saturated with gas at the equilibrium pressure (Zhang et al., 2000; Renner, 1988):
 ð1Þ C ðz ¼ 0; t 4 0Þ ¼ C sat P eq

36

F.J. Pacheco Roman, S.H. Hejazi / Chemical Engineering Science 146 (2016) 35–43

The quasi-equilibrium boundary condition assumes that the interface is saturated with gas at the existing pressure in the gas phase (Eq. (2)) instead of the equilibrium pressure (PachecoRoman and Hejazi, 2015; Upreti and Mehrotra, 2000, 2002; Sheikha et al., 2005, 2006). C ðz ¼ 0; t 40Þ ¼ C sat ½PðtÞ

ð2Þ

The non-equilibrium boundary condition considers resistance to mass transfer at the gas–oil interface. The mathematical expression, Eq. (3), states that the rate of transfer at the interface is proportional to the difference between the saturation concentra
 tion (at equilibrium pressure), C sat P eq , and the existing concentration at the interface, C ðz; t Þjz ¼ 0 :  
  ∂C  D z ¼ 0 ¼ k C sat P eq  C ðz; t Þz ¼ 0 ∂z

ð3Þ

Civan et al. first modeled the pressure-decay experiments with Eq. (3) and proposed a graphical procedure to estimate D and k from limited pressure-decay data (i.e., data not recorded until reaching the equilibrium pressure) (Civan and Rasmussen, 2006). Later on, they included the equilibrium pressure and solubility in the parameter-estimation procedure (Civan and Rasmussen, 2009; Rasmussen and Civan, 2009). Subsequently, Etminan et al. (2013) proposed a new version of the non-equilibrium boundary condition, Eq. (4): D

   ∂C   z ¼ 0 ¼ k C s ðt Þ  C ðz; t Þ z ¼ 0 ∂z

ð4Þ

In Eq. (4), C s ðt Þ is the time-dependent interface concentration at equilibrium with the vapor pressure in the gas phase, whereas C ðz; t Þjz ¼ 0 is the concentration at the interface evaluated from the oil phase. In contrast to Eq. (3), Eq. (4) considers a time-dependent concentration at the interface rather than a constant value at equilibrium pressure. On the basis of this modified boundary condition, Etminan et al. (2014) proposed nonlinear regression to estimate H ij , D, and k. They reported that assuming a constant saturation concentration underestimates the rate of gas dissolution which consequently underpredicts diffusion coefficients. Although useful, nonlinear regression requires computational expertise as various numerical algorithms are required. This demonstrates the difficulty in evaluating these parameters and the need of producing simple procedures to obtain quick, yet acceptable estimates. There are various graphical or history-matching techniques to evaluate D while a straightforward graphical technique to estimate k is lacking. Therefore, we focus on the development of an analytical-graphical approach for the easy estimation of k from PD data. In this study, we obtain a closed-form solution in the time domain utilizing the Integral Method (IM) (Goodman, 1964) for the time-dependent interface boundary condition. This solution allows the development of a graphical technique to extract the interface mass-transfer coefficient and Henry's constant provided the diffusion coefficient is known. In fact, the proposed solution, which is derived based on the time dependent non-equilibrium boundary condition, is a generalization of our previous model developed for quasi-equilibrium boundary condition in the absence of interface mass resistance (Pacheco-Roman and Hejazi, 2015). The structure of this paper is as follows: first, we present the derivation of the analytical solution (i.e., forward problem) and its validation. Then, the graphical method is presented and applied to synthetic and experimental data to estimate mass-transfer coefficients (i.e., the inverse problem). Finally, we present the results and discussion followed by the conclusions.

Fig. 1. Pressure-decay cell schematic.

2. Forward problem: deduction of the mathematical model 2.1. Mathematical formulation of the problem The diffusion of gas into oil in the pressure/volume/temperature (PVT) cell as schematized in Fig. 1 is modeled using Fick’s second law (Eq. (5)) (Fick, 1855): ∂C ∂2 C ¼D 2 ∂t ∂z

ð5Þ

The use of Fick's second law to model the pressure-decay experiments has been discussed in different studies (Zhang et al., 2000; Sheikha et al., 2005; Tharanivasan et al., 2006; PachecoRoman and Hejazi, 2015; Etminan et al., 2013). The concentration of gas in heavy oil is negligible at the beginning of the process. Therefore, the initial condition is expressed as: C ðt ¼ 0; z Z 0Þ ¼ 0

ð6Þ

The interface boundary condition (Eq. (4)) assumes resistance to mass transfer at the gas-heavy oil interface ðz ¼ 0Þ which can be expressed in the form of Eq. (7). The detailed derivation of Eq. (7) is presented in Appendix A (Etminan et al., 2013) and the parameters are defined in the nomenclature section.    ψ Hij dC  ψ Hij ∂2 C  ∂C  ¼  ð7Þ z ¼ 0 z ¼ 0 ∂z  D dt  k ∂t∂zz ¼ 0 where:

ψ¼

LM ZRT

ð8Þ

The bottom of the cell is finite and closed ðz ¼ hÞ, hence a no flow boundary condition is considered:   ∂C ¼0 ð9Þ ∂z z ¼ h 2.2. Solution of the mathematical formulation The mathematical model given by Eqs. (5)–(7) and (9) has an exact solution in the Laplace domain (Etminan et al., 2013). The solution is evaluated at the interface ðz ¼ 0Þ and expressed in terms of pressure with Henry's Law to obtain Eq. (10): h

pffiffiffi pffiffis  pffiffi pffiffiffii Hij ψ s  2h Ds 1 þe  2h D  Dk  Ds De D Pi   pffiffis

ð10Þ P ðsÞ ¼ H ψ ffiffi ffi p H ij ψ H ij ψ pffiffisffi ij  2h D H ij ψ s D s þð1 þ k sÞ D þ e D s  ð1 þ k sÞ D The above solution requires a numerical-inversion technique to obtain the values of pressure in the real-time domain. Therefore, it cannot be used for developing a graphical technique to estimate interface mass transfer coefficient. For this purpose, an alternative closed-form analytical solution for the problem is presented.

F.J. Pacheco Roman, S.H. Hejazi / Chemical Engineering Science 146 (2016) 35–43

37

We solve the mathematical formulation given by Eqs. (5)– (7) and (9) by use of the Integral Method (Goodman, 1964). The IM allows to obtain an approximate analytical solution that is also evaluated at the interface ðz ¼ 0Þ and expressed in terms of pressure vs. time (Eq. (11)). The detailed derivation of this solution and its validation are presented in Appendix B and C, respectively.  13 2 0

1 þ ψ hHij
 Pi 4 h t A5 ¼ m2;2 exp m1;2 t þ b2;2 P ðt Þ ¼ 1þ exp@ 

2 h 5 h ψ Hij 1 þ ψhHij þ k 12 D

Therefore, k is determined from the slope m1;2 : As mentioned before, there are some methods dependent on PDT data to evaluate D, while techniques to estimate k are scarce. Assuming that we know the value of D, the interface mass transfer coefficient is solved from Eq. (13) as:       m 1;2  ð15Þ k ¼  

1 1 5 h  þ þ 12 Dm1;2  h ψ H ij

ð11Þ

Since the diffusion coefficient is required to evaluate the interface mass-transfer coefficient ðkÞ we propose the estimation of D by the use of some other graphical methods available in the literature. For this purpose, we use the method developed by Sheikha et al. (2006). This technique requires plotting of early time h i pffiffi  1 P ðt Þ vs. t . Diffusion coefficient is data in the form of erfc Pi
2 evaluated through the slope of straight line, m, as D ¼ mψ H ij .

3. Inverse problem: extraction of mass transfer parameters The inverse problem consists of the estimation of Henry's constant and interface mass-transfer coefficient by the solution proposed in Eq. (11) along with the data from pressure-decay tests.

3.4. Accuracy and Sensitivity of k to diffusion coefficient

3.1. Development of the graphical method The inverse problem is solved by developing a graphical approach derived from the manipulation of the IM solution (Eq. (11)). The first step of our method requires plotting the experimental data as suggested by the linear Eq. (12): 3 2



  1 þ ψ hHij ψ h dP ðt Þ 12Dk 5t þ ln  ln  ¼  4

¼ m1;2 t  b1;2 2 2 h P i dt 12Dh þ5h k þ5h k

12 D

ð12Þ Eq. (12) is obtained by differentiating Eq. (11) with respect to time, applying natural logarithm, and simplifying accordingly. The purpose of this plot is to determine the slope m1;2 (Eq. (13)) to be used for creating a second plot.

 1 þ ψ hHij  ð13Þ m1;2 ¼ 

2 h þ5h k 12 D
 The next step is to plot P ðt Þ vs. exp  m1;2 t , as given by Eq. (11) and evaluate the slope and intercept, m2;2 and b2;2 , respectively. 3.2. Estimation of Henry’s constant The slope ðm2;2 Þ and the intercept ðb2;2 Þ of the resulting straight line of this second plot only contain one unknown parameter (i.e., Henry's constant). Therefore, Henry's constant can be estimated either from the slope or the intercept. However, a previous study shows that the intercept rather than the slope yields more accurate estimations (Pacheco-Roman and Hejazi, 2015). The intercept in Eq. (11) represents the pressure as time tends to infinity (i.e., at equilibrium conditions), which aligns with Henry’s law stating the ultimate amount of gas dissolved in heavy oil at equilibrium. Thus, the intercept is used for the estimation of Henry's constant.   h b2;2 H ij ¼ ð14Þ ψ P i  b2;2 3.3. Estimation of mass transfer coefficient Once the value of H ij is determined, k could be determined from the slope or intercept of Eq. (12). However, the slope m1;2 is proportional to the rate of gas-phase pressure change or equivalently to the rate of gas diffusion in the heavy oil, dictated by the diffusion and the interface mass-transfer coefficients (Pacheco-Roman and Hejazi, 2015). On the other hand, the intercept of Eq. (12) does not contain a physical meaning.

We apply the proposed inverse solution for graphical estimation of k to synthetic pressure-decay data produced with the exact solution (Eq. (10)). To generate the data set, available experimental parameters for the system carbon dioxide-heavy oil shown in Table 1 are used (Etminan et al., 2014). Circles in Fig. 2 shows the generated pressure data. Henry's constant estimated with the proposed method in Section 3.2 is 78357 Pa∙m3 =kg which is in agreement with the value originally used for synthetic data generation. Eq. (15) is used to estimate k with the slope m1;2 and original value of D. The evaluated k as 2:473  10  7 m=s agrees well with the original mass transfer coefficient. In fact for the exact value of D, our proposed technique can estimate k with less than 3% error which demonstrates the strength of the proposed methodology. In the second step, we recalculate k with the estimated diffusion coefficient using Sheikha et al. (2006) method. Diffusion coefficient estimated by this technique is 1.07  10  10 m2/s which is 20% less than the original value. Using this D; interface mass-transfer coefficient is estimated as 1.170  10  7 m/s; 54% less than the original value. To demonstrate the influence of estimated k on the pressure profiles, we plot the pressure decay curve in Fig. 2 along with the synthetic data. Since the estimated k from the value of D obtained through Sheikha et al. (2006) method is less than the original value, there is more resistance in mass transport. Thus, the pressure profile does not decay as fast as the original one declines. We also perform sensitivity analysis to distinguish the effect of k vs. D on the pressure-decay curves. Assuming that the exact values of H ij and D are known, we assign different values for k, Table 1 Pressure-decay parameters of the system CO2/heavy oil used to generate the synthetic data. Parameter

Value

Unit

Diffusion coefficient, ðDÞ

1.339  10-10 2.546  10-7 30:52 78,400 44.01 0.7779 3.5299 297.05 0.01605 0.03565

m2 =s m=s Dimensionless Pa m3/kg kg/kmol dimensionless MPa K m m

Mass transfer coefficient ðkÞ Biot mass transfer number ðBiÞ Henry's constant, ðHij Þ Molecular mass, ðMÞ Compressibility factora, ðZÞ Initial pressure ðP i Þ Temperature ðTÞ Thickness of bitumen layer, ðhÞ Gas zone in PVT cell ðLÞ

a Calculated with the Peng–Robinson EOS at the initial pressure (Peng and Robinson, 1976).

38

F.J. Pacheco Roman, S.H. Hejazi / Chemical Engineering Science 146 (2016) 35–43

3.6 3.5

k from D of Sheikha et al. 2006

3.4

Pressure, MPa

validating the accuracy of the results is more meaningful than comparing specific values reported. In the next section, we apply the developed graphical technique to estimate the mass transfer parameters for some experimental data available in the literature.

Synthetic Data (original k and D)

k with the exact value of D

3.3 3.2

4. Results and discussions: application of the graphical method to experimental data

3.1 3.0 2.9 2.8

0

200

400

600

800

1000

1200

Time, hr Fig. 2. Comparison of pressure-decay profiles generated with the original and estimated parameters.

3.6 3.5

k = 2.546E-9 k = 2.546E-8 k = 2.546E-7 k = 2.546E-6 k = 2.546E-5 k = 2.546E-4

k= 0

Pressure, MPa

3.4 3.3 3.2

Increasing k

3.1 3.0 2.9

k →∞

2.8 0

200

400

600

800

1000

1200

Time, hr Fig. 3. Sensitivity analysis of k in the original pressure decay profile.

smaller and larger than the original one from k ¼ 0 to k-1. The corresponding pressure-decay curves are depicted in Fig. 3. A value of k ¼ 0 means that there is no mass transfer in the system. Hence, the gas-phase pressure does not decay. Fig. 3 demonstrates that as we increase the value of k, a typical pressuredecay curve results. The pressure-decay curves are less sensitive to the value of k when we use values larger than original k. When k becomes very large, i.e. k-1; there is no resistance to mass transfer which is acceptable at late times (Tharanivasan et al., 2004). For values of k, which are one order of magnitude or more larger than the reference value and also when k-1, the predicted pressure data perfectly match the synthetic pressure data at late times, while there is a slight underestimation at early times. The deviations from the synthetic data become significant as the order of magnitude of k becomes smaller than the original value, which means that there is an increasing resistance to mass transfer at the interface. It should be noted that because k and D have very small values, their precise estimation by the use of history matching or objective-function minimization is difficult. Tharanivasan et al. (2004) recognized the difficulty of finding an optimal value for D and k as both parameters are also coupled in their solution in terms of the Biot mass-transfer number. Etminan et al. (2014) also mentioned the problem of finding these two parameters for which values are very small. They acknowledged that parameters with large magnitudes such as Henry's constant are the most desirable so that the inverse problem is not very sensitive to measurement errors. In summary, comparing the order of magnitude of the interface mass transfer coefficients for

We apply our method to estimate Henry's constant ðH ij Þ and the interface mass-transfer coefficient ðkÞ to the experimental data of Zhang et al. (2000), Tharanivasan et al. (2006) and Behzadfar and Hatzikiriakos (2014) for carbon dioxide-bitumen. Properties of this gas-bitumen system and other experimental conditions of these tests are presented in Table 2. The inverse solution, previously presented for the synthetic data, is applied to these experimental data. Thus, the method is shown concisely with just some comments to consider when we deal with this type of data. First, the data are plotted according to Eq. (12) as shown in Fig. 4. The experimental data are smoothed before performing the numerical differentiation. The slope obtained from the plot of Fig. 4 from each data set is used to construct the second plot according to Eq. (11) and presented in Fig. 5: The intercept of this second plot along with Eq. (14) is used to calculate H ij for each experimental data set. Subsequently, the diffusion coefficient is calculated from the graphical method of Sheikha et al. (2006) as shown in Fig. 6. Finally, the interface mass-transfer coefficient is determined with Eq. (15) considering the estimated Henry's constant and diffusion coefficient. Table 3 presents the results of the estimated H ij ; D, and k in addition to the reported values from the literature for these experimental data by Civan and Rasmussen (2006), Tharanivasan et al. (2006) and Etminan et al. (2013). It is important to remind that to evaluate these coefficients, Etminan et al. (2013) use history matching with the timedependent non-equilibrium boundary condition (Eq. (4)), Civan and Rasmussen (2006) use their developed graphical methods and Tharanivasan et al. (2006) use history matching, both considering the non-equilibrium boundary condition with a constant con
 centration value C sat P eq , (Eq. (3)). Please note that the estimate values for k is some cases are one order of magnitude different from other references. This is mainly due to the differences in the diffusion coefficients used in these studies. On the other hands, all reported values are estimations and depends upon the assumptions and used boundary conditions. To evaluate the validity of estimated parameters, we regenerate the pressure decay curves Table 2 Experimental pressure-decay parameters for CO2/bitumen. Parameter

Zhang et al. (2000)

Tharanivasan et al. (2006)

Behzadfar and Hatzikiriakos (2014)

Unit

Molecular mass, ðMÞ Compressibility factora, ðZÞ Initial pressure ðP i Þ Temperature ðTÞ Thickness of bitumen layer, ðhÞ Gas zone in PVT cell ðLÞ

44.01

44.01

44.01

kg/kmol

0.7765

0.7282

0.7636

–

3.435 294.15 0.70

4.169 297.05 0.0435

5.008 323.15 0.0314

MPa K m

0.18

0.1165

0.0559

m

a Calculated with the Peng–Robinson EOS at the initial pressure (Peng and Robinson, 1976).

F.J. Pacheco Roman, S.H. Hejazi / Chemical Engineering Science 146 (2016) 35–43

-21.0

ln[(-(ψh/Pi )dP/dt]

using the estimated ones. Fig. 7a and b represent the experimental pressure-decay curves and the ones generated by the estimated and reported H ij , k and D in the exact solution (Eq. (10)). The pressure-decay profiles from our estimated results are in agreement to the experimental data, which supports the accuracy of the estimated parameters. Thus, it validates the potential of our proposed method in the estimation of H ij and k. On the basis of the results obtained with the synthetic and experimental data, our graphical procedure is deemed capable for an easy and reliable estimation of H ij and k in comparison to other available methods. For instance, the implementation of the method of Etminan et al. (2014) requires the use of various numerical algorithms to determine these mass-transfer parameters. The graphical method of Civan and Rasmussen (2009) has the advantage of using time-limited data in their proposed solution (i.e., data not reaching the equilibrium pressure) and asserts to use a nonconstant-compressibility factor. However, the assumption of constant saturation at the gas–oil interface is a drawback according to Etminan et al., (2013). The use of Henry's Law in our graphical procedure limits its application to gases sparingly soluble in solvents under given conditions of pressure and temperature as discussed in other works (Sheikha et al., 2005; Etminan et al., 2013). It is widely accepted that resistance to mass transfer at the carbon dioxide/ heavy-oil interface should not be neglected, while it can be neglected in methane/heavy-oil systems (Tharanivasan et al., 2004; Etminan et al., 2013). Therefore, we only analyzed data for carbon dioxide/-heavy oil. Also, the selection of only a short portion of data might not be representative in real-experimental-data analyses because of noise present in the system. Therefore, our graphical approach to estimatek andD requires early-time and sufficient late-time data. In spite of the simplifying assumptions considered in the mathematical formulation of the problem, it has been shown that the present graphical method can use experimental data and provide acceptable results for the estimation of interface masstransfer coefficient.

Data (Zhang et al., 2000) Data (Tharanivasan et al., 2006) Data (Behzadfar et al., 2014) Representative straight line

-22.0 -23.0

y = -0.0091x - 21.506 R² = 1

-24.0

y = -0.0022x - 24.153 R² = 1

-25.0

y = -0.0117x - 23.643 R² = 0.9512

-26.0 -27.0 0

100

200

300

400

500

600

Time, hr Fig. 4. Plot of experimental pressure-decay data according to Eq. (12).

5.3 Data (Zhang et al., 2000) Data (Tharanivasan et al., 2006)

Pressure, MPa

4.8

Data (Behzadfar et al. 2014) y = 0.6694x + 4.3053 R² = 0.9954

Representative straight line

4.3

3.8

y = 0.5056x + 3.4834 R² = 1

3.3

y = 0.5925x + 2.8425 R² = 1

2.8 0.0

0.2

0.4

0.6

0.8

1.0

exp(m 1t) Fig. 5. Plot of experimental pressure-decay data according to Eq. (11).

0.12

Data (Zhang et al., 2000) Data (Tharanivasan et al., 2006)

erfc-1[P(t)/Pi ]

0.10

5. Conclusions

Data (Behzadfar et al., 2014) Representative straight line

0.08

y = 0.0122x -0.0304 R² = 0.9999

0.06 0.04

y = 0.0096x -0.0113 R² = 0.9991

0.02 0.00

39

y = 0.0071x -0.0094 R² = 0.9985

0

2

4

6

8

10

12

Square root of time, hr 0.5 Fig. 6. Sheikha et al.'s graphical method (2006) to estimate D from experimental data.

In this study, the diffusion of gases into heavy oil is modeled by use of Fick's second law and a time-dependent boundary condition considering the gas–oil-interface resistivity. The mathematical formulation of the problem solved by the Integral Method (IM) provides a closed-form analytical solution in the time domain suitable for the prediction of the declining pressure as a function of time in the pressure-decay technique (PDT). The resultant solution derived from a time-dependent interface concentration is adequate for the development of a graphical procedure in which the mass-transfer coefficient, in addition to Henry's constant, are estimated with the diffusion coefficient obtained from an auxiliary graphical method. The results obtained from synthetic and experimental pressure-decay data

Table 3 Comparison of parameters estimated though the current graphical procedure Parameter

Experimental data of Zhang et al. (2000)

Current procedure

H ij , (Pa m3/kg) 80,506 1.746  10  6 k, (m/s) 4.663  10  9 D, (m2/s)

Experimental data of Tharanivasan et al. (2006)

Experimental data of Behzadfar and Hatzikiriakos (2014)

Etminan et al. Civan and (2013) Rasmussen (2006)

Zhang et al. Current (2000) procedure

Etminan et al. (2013)

Tharanivasan et al. 2006

Current procedure

Behzadfar and Hatzikiriakos (2014)

85,500 1.548  10  7 6.202  10  8

NA 1 4.8  10  9

78735 4.369  10  6 3.558  10  10

NA 4 1.31  10  6 5.7  10  10

159,265 3.79  10  7 9.33  10  10

NA NA 1.08  10  9

NA 3.26  10  7 1.13  10  8

77575 5.398  10  8 6.841  10  10

40

F.J. Pacheco Roman, S.H. Hejazi / Chemical Engineering Science 146 (2016) 35–43

3.5

4.2 Experimental data (Zhang et al.,2000) Estimated D, k & Hij

3.3 3.2 3.1

3

Estimated D, k & Hij

4 3.9 3.8 3.7 3.6

2.9 2.8

Experimental data (Tharanivasan et al. 2006)

4.1

Pressure, Mpa

Pressure, MPa

3.4

0

5

10

15

3.5

20

0

10

20

Time, days

30

40

Time, days

5.1 Experimental data (Behzadfar et al. 2014)

5

Pressure, Mpa

Estimated D, k & Hij

4.9 4.8 4.7 4.6 4.5

0

10

20

30

40

50

60

Time, hours Fig. 7. Comparison of experimental and predicted pressure-decay curves of CO2/bitumen.

analyses demonstrate the strength of our proposed graphical method for the estimation of k in gas/heavy-oil diffusion process without making it an adjustable parameter.

mg mi mo m1,2 m2,2

Nomenclature A b b1,2 b2,2 Bi B1 C Cs C sat C si D h H ij k L m M md mf

cross-sectional area of the pressure cell, L2 ; m2 intercept of a graphical method intercept of straight line from first plot of the graphical method, dimensionless intercept of straight line from second plot of the graphical method, m=Lt2 ; Pa Biot mass transfer number calculated as hk=D, dimensionless constant of the forward solution given by
  P i hψ m1;2 =D h þ ψ H ij , m=L3 ; kg=m3 concentration of gas in heavy oil, m=L3 ; kg=m3 concentration of gas at the mass transfer interface, m=L3 ; kg=m3 saturation concentration of gas in heavy oil, m=L3 ; kg=m3 initial concentration of gas at the interface, m=L3 ; kg=m3 diffusion coefficient, L2 =t; m2 =s thickness of heavy oil layer in the pressure cell, L; m Henry's constant of solute i in solvent j, L2 =t2 ; Pa  m3 =kg interface mass transfer coefficient, m=s height of the gas zone in the pressure cell, L; m slope of a graphical method molar mass of gas, m=n; kg=kmol dissolved mass of gas in the heavy oil, m; kg final mass of gas in the gas cap, m; kg

nm P P P avg P eq Pi R s t T Vg z Z

mass of gas in the gas cap, m; kg initial mass of gas in the gas cap, m; kg mass of heavy oil, m; kg slope of the straight line from the first plot of the graphical method, t  1 ; s  1 slope of straight line from the second plot of the graphical method, m=Lt2 ; Pa number of moles, n; kmol pressure, m=Lt2 ; Pa pressure in the Laplace domain average pressure, m=Lt2 ; Pa equilibrium pressure, m=Lt2 ; Pa initial pressure, m=Lt2 ; Pa universal gas constant, R ¼ 8314Pa  m3 =kmol  K Laplace transform variable time, t; s temperature, T; K volume of gas, L3 ; m3 depth, L; m gas compressibility factor, dimensionless

Greek symbols

ψ

constant value calculated as LM=ZRT, t2 =L; kg=Pa  m2

Acknowledgment The authors appreciate the constructive comments received from anonymous journal reviewers. The first author expresses his

F.J. Pacheco Roman, S.H. Hejazi / Chemical Engineering Science 146 (2016) 35–43

gratitude to the Mexican National Council for Science and Technology (CONACYT), the Secretariat of Public Education (SEP) and the Mexican Government for the awarded scholarships to pursue his graduate studies. The authors also acknowledge financial support from the Natural Sciences and Engineering Research Council of Canada (NSERC).

Appendix A. Transformation of the non-equilibrium boundary condition (Eq. (4)) The general form of the interface boundary condition (Eq. (4)) is rewritten as Eq. (7) following the procedure of Etminan et al. (2013). First, Eq. (4) in the main text is derived with respect to time and the terms rearranged to obtain Eq. (A.1):   ∂C s ∂C  D ∂2 C  ¼ z ¼ 0  ∂t k ∂t∂zz ¼ 0 ∂t

ðA:1Þ

Eq. (A.1) is written uniformly in terms of C (concentration of gas at the interface in contact with heavy oil) through gas-phase mass-balance equations by use of the real-gas equation of state
 PV g ¼ Znm RT . Then, the mass of gas leaving the gas cap is expressed as: ∂mg V g M dP ¼ ZRT dt ∂t

ðA:2Þ

The mass of gas transferred or diffused into the heavy oil is be expressed in terms of Fick's first law as:  ∂md dC  ¼ DA  dz z ¼ 0 ∂t

ðA:3Þ

Equating Eqs. (A.2) and (A.3) yields to:  V g M dP dC  ¼  DA  ZRT dt dz z ¼ 0

ðA:4Þ

The pressure in Eq. (A.4) is expressed in terms of C s (concentration of gas at the interface in contact with the gas cap) by use of Henry's law: P ¼ H ij C s

ðA:5Þ

Thus, substituting Eq. (A.5) into Eq. (A.4) and rearranging the terms yields to Eq. (A.6):  dC s DAZRT dC  ¼ H ij V g M dz z ¼ 0 dt

ðA:6Þ

Therefore, the left hand side of the boundary condition (Eq. (4)) can be expressed as:    DAZRT dC  ∂C  D ∂2 C  z¼0 ¼ z¼0   H ij V g M dz ∂t k ∂t∂zz ¼ 0

   H ij V g M dC  H ij V g M ∂2 C  ∂C  z¼0 ¼ z¼0   ∂z DAZRT dt kAZRT ∂t∂zz ¼ 0

Appendix B. Derivation of the Finite-Acting Solution by the Integral Method (IM) for the Non-equilibrium Boundary Condition at the Interface The integration of Eq. (5) over the domain of interest ð0; hÞ results in Eq. (B.1), referred to as mass-balance integral: Z ∂C ∂C 1 d h  ¼ C ðz; t Þdz ðB:1Þ ∂z z ¼ h ∂z z ¼ 0 D dt 0 The IM is characterized by the assumption of a trial function for the dependent variable of the differential equation, which in this case is the concentration C. Precisely, one drawback of the IM is the arbitrary selection of the approximating function, which influences the accuracy of the method (Mitchel and Mayer, 2012). In this case, an incomplete cubic polynomial is chosen for the concentration profile Cðz; tÞ in the form of: C ðz; t Þ ¼ a0 ðt Þ þ a1 ðt Þz þ a3 ðtÞz3

ðA:8Þ

ðB:2Þ

Eq. (B.2) is used along with the inner and outer boundary conditions (Eqs. (7) and (9)) to obtain the following equations: a1 ðt Þ ¼

ψ H ij D

a3 ðt Þ ¼ 

a00 ðt Þ 

ψ H ij k

a01 ðt Þ

ðB:3Þ

a1 ðt Þ 3h

ðB:4Þ

2

Then, the coefficients of the concentration profile can be rewritten only in terms of a1 ðt Þ by substituting Eqs. (B.3) and (B.4) into Eq. (B.2):   Z t D D z3 D þ z  2 þ a0 ð0Þ  a1 ð0Þ a1 ðt Þdt þ a1 ðt Þ ðB:5Þ C ðz; t Þ ¼ k k ψ Hij 0 3h The first term of Eq. (B.1) is zero according to the outer boundary condition (Eq. (9)). Then, the concentration profile given by Eq. (B.5) is introduced in both sides of the mass-balance integral (Eq. (B.1)) to obtain the following:

 1 þ ψ hHij ¼0 ðB:6Þ a0 1 ðt Þ þ a1 ðt Þ

2 h þ5h k 12 D Therefore, the IM has reduced the PDE problem to a first order ODE for the coefficient a1 ðt Þ, (Eq. (B.6)), so:
 ðB:7Þ a1 ðt Þ ¼ B1 exp  m1;2 t Where:

 1 þ ψhHij  m1;2 ¼

2 h þ5h k 12 D

ðB:8Þ

Integration of Eq. (B.3) gives the expression for the concentration at the interface C ð0; t Þ ¼ ao ðt Þ: Z t D D D a ðt Þdt þ a1 ðt Þ þ a0 ð0Þ  a1 ð0Þ ðB:9Þ a0 ðt Þ ¼ k k ψ H ij 0 1 Substituting Eq. (B.7) into Eq. (B.5) and simplifying accordingly:




 1 1 DB1 DB1 þ C s ðt Þ ¼ a0 ðt Þ ¼ DB1 exp  m1;2 t  þ a0 ð0Þ  ψ H ij m1;2 k k ψ H ij m1;2

ðB:10Þ

ðA:7Þ

Finally, Eq. (7) is obtained when the terms in Eq. (A.8) are rearranged and the rests of parameters grouped as constant ψ :

41

An expression in terms of pressure is obtained from the inner
 boundary condition (Eq. (4)) by use of Henry's Law P ¼ H ij C : P ðt Þ ¼ 

DHij a1 ðt Þ þ H ij a0 ðt Þ k

ðB:11Þ

Then, Eqs. (B.10) and (B.7) are substituted into (B.11) and the constants a0 ð0Þ and B1 are determined by evaluating the resultant expression at time zero and time infinity.

42

F.J. Pacheco Roman, S.H. Hejazi / Chemical Engineering Science 146 (2016) 35–43

3.5

Since, the pressure at time zero is equal to the initial pressure, the constant ao ð0Þ is expressed as:

P i DB1 ðB:12Þ a0 ð0Þ ¼ þ H ij k

P hψ m1;2  B1 ¼  i D h þ ψ H ij

ðB:13Þ

The expression for P eq is found from mass conservation law (Eq. (B.14)) which states that the initial mass of gas in the pressure/ volume/temperature (PVT) cell is equal to the final mass of gas remaining in the gas cap plus the mass of gas diffused into heavy oil as determined by Henry's law. mi ¼ mf þ md

ðB:14Þ

PiV g M ¼ P i Aψ ZRT

The saturation concentration according to Henry's law is used to determine the final mass of gas dissolved in the volume of heavy oil ðV o ¼ AhÞ: md ¼

P eq Ah H ij

ðB:17Þ

P i ψ H ij h þ ψ H ij

ðB:18Þ

Once the constants B1 and a0 ðt Þ are known, these are substituted in Eq. (B.11) to obtain the analytical solution presented in the main text as Eq. (11):  13 2 0

1 þ ψhHij Pi 4 h  t A5 ðB:19Þ P ðt Þ ¼ 1þ exp@ 

2 h ψ Hij 1 þ ψhH þ5h ij

k

12 D

When k tends to infinity the non-equilibrium model (Eq. (11)) reduces to the one previously developed under equilibrium conditions (Pacheco-Roman and Hejazi, 2015).

Appendix C. Validation of the forward solution The validity of the IM solution (Eq. 11) is verified by generating the pressure-decay profile and comparing it to the one obtained by the exact solution in the Laplace domain (Eq. (10)). We use the Table C1 Pressure-decay parameters of carbon dioxide/heavy oil used to generate Fig. C1. Parameter Diffusion coefficient, ðDÞ Mass transfer coefficient ðkÞ Biot mass transfer number ðBiÞ Henry’s constant, ðHij Þ Molecular mass, ðMÞ Compressibility factora ðZÞ Initial pressure ðP i Þ Temperature ðTÞ Thickness of bitumen layer, ðhÞ Gas zone in PVT cell ðLÞ

0

100

Value

200

300

400

500

Fig. C1. Pressure-decay curves of carbon dioxide/heavy oil from the exact solution (Eq. (10)) and the IM solution (Eq. 11).

t = 500 hrs

35.0 30.0

t = 250 hrs

25.0 20.0

t = 100 hrs

15.0

Exact solution

10.0

IM solution

t = 10 hrs

5.0

therefore: P eq ¼

3.0

Time, hr

Concentration. kg/m3

ðB:16Þ

3.1

2.8

ðB:15Þ

mf ¼ P eq Aψ

3.2

2.9

The initial and final mass of gas in the gas cap is given by these equations according to the real gas law: mi ¼

IM Solution

3.3

Pressure, MPa

Evaluation as time tends to infinity, which is equal to the equilibrium pressure P eq , gives the expression for B1 :

Exact solution by Stehfest algorithm

3.4

0.0

0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

Thickness of bitumen layer, m Fig. C2. Concentration profiles from the exact and IM solution of carbon dioxide along heavy oil layer for 10, 100, 250 and 500 h.

Stehfest algorithm (Stehfest, 1970) to obtain pressures in the realtime domain. The parameters to produce the curves (Table C1) are taken from a pressure-decay test reported in the literatura (Tharanivasan et al., 2006). Fig. C1 shows that the pressure-decay curves generated from the IM solution (solid line) accurately matches the ones generated by the Laplace solution (circles), which demonstrates the capability of approximate IM solution to be used instead of the exact one. Thus, the development of this analytical solution (Eq. 11) to predict measurements (i.e., pressure as a function of time) is considered suitable for the purpose of inverse-problem development. We also examine the capability of this approximate analytical solution by generating concentration profiles. Fig. C2 presents the concentration profiles of gas in heavy oil derived with the exact (circles) and the IM solution (solid lines) at different times, which demonstrates the capability of approximate IM solution.

Unit 8

6.202  10 1.548  10  7 0.174 85,500 44.01 0.7761 3.435 294.15 0.07 0.18

m2/s m/s dimensionless Pa m3/kg kg/kmol Dimensionless MPa K m m

a Calculated with the Peng–Robinson EOS at the initial pressure (Peng and Robinson, 1976).

References Behzadfar, E., Hatzikiriakos, S.G., 2014. Diffusivity of CO2 in bitumen: pressuredecay measurements coupled with rheometry. Energy Fuels 28 (02), 1304–1311. http://dx.doi.org/10.1021/ef402392r. Civan, F., Rasmussen, M.L., 2006. Determination of gas diffusion and interface-mass transfer coefficients for quiescent reservoir liquids. SPE J. 11 (01), 71–79. http://dx.doi.org/10.2118/84072-PA. Civan, F., Rasmussen, M.L., 2009. Rapid simultaneous evaluation of four parameters of single-component gases in nonvolatile liquids from a single data set. Chem. Eng. Sci. 64 (23), 5084–5092. http://dx.doi.org/10.1016/j.ces.2009.08.021. Etminan, S.R., Maini, B.B., Chen, Z., 2014. Determination of mass transfer parameters in solvent-based oil recovery techniques using a non-equilibrium boundary

F.J. Pacheco Roman, S.H. Hejazi / Chemical Engineering Science 146 (2016) 35–43

condition at the interface. Fuel 120, 218–232. http://dx.doi.org/10.1016/j. fuel.2013.11.027. Etminan, S.R., Pooladi-Darvish, M., Maini, B.B., Chen, Z., 2013. Modeling the interface resistance in low soluble gaseous solvents-heavy oil systems. Fuel 105, 672–687. http://dx.doi.org/10.1016/j.fuel.2012.08.048. Fick, A., 1855. On liquid diffusion. Philos. Mag. Ser. 4 10 (63), 30–39. http://dx.doi. org/10.1080/14786445508641925. Goodman, T.R., 1964. Application of integral. Methods Transient Nonlinear Heat 1, 51–122. http://dx.doi.org/10.1016/S0065-2717(08)70097-2. Mitchel, S.L., Myers, T.G., 2012. Application of heat balance integral methods to onedimensional phase change problems. Int. J. Differ. Equ. 2012, 1–22. http://dx. doi.org/10.1155/2012/187902. Pacheco-Roman, F.J., Hejazi, S.H., 2015. Estimation of solubility and diffusivity of gases in heavy oils using late-time pressure-decay data: an analytical-graphical approach. SPE J., 1–12. http://dx.doi.org/10.2118/170957-PA (Preprint). Peng, D.Y., Robinson, D.B., 1976. A new two-constant equation of state. Ind. Eng. Chem. Fundam. 15 (1), 59–64. http://dx.doi.org/10.1021/i160057a011. Rasmussen, M.L., Civan, F., 2009. Parameters of gas dissolution in liquids obtained by isothermal pressure decay. Fluid Mech. Transp. Phenom. 55 (1). http://dx. doi.org/10.1002/aic.11669. Renner, T., 1988. Measurement and correlation of diffusion coefficients for CO2 and rich-gas applications. SPE Reserv. Eng. 3 (2), 517–523. http://dx.doi.org/10.2118/ 15391-PA. Riazi, M.R., 1996. A new method for experimental measurement of diffusion coefficients in reservoir fluids. J. Pet. Sci. Eng. 14 (3–4), 235–250. http://dx.doi.org/ 10.1016/0920–4105(95)00035-6. Riazi, M.R., Whitson, C.H., Da Silva, F., 1994. Modelling of diffusional mass transfer in naturally fractured reservoirs. J. Pet. Sci. Eng. 10 (3), 239–253. http://dx.doi. org/10.1016/0920–4105(94)90084-1.

43

Sheikha, H., Mehrotra, A.K., Pooladi-Darvish, M., 2005. Development of graphical methods for estimating the diffusivity coefficient of gases in bitumen from pressure-decay data. Energy Fuels 19 (5), 2041–2049. http://dx.doi.org/10.2118/ 75135-MS. Sheikha, H., Mehrotra, A.K., Pooladi-Darvish, M., 2006. An inverse solution methodology for estimating the diffusion coefficient of gases in athabasca bitumen from pressure-decay data. J. Pet. Sci. Eng. 53 (3–4), 189–202. http://dx.doi.org/ 10.1016/j.petrol.2006.06.003. Stehfest, H., 1970. Algorithm 368: numerical inversion of laplace transforms. Commun. ACN 13 (1), 47–49. http://dx.doi.org/10.1145/361953.361969. Tharanivasan, A.K., Gu, Y., Chaodong, Y., 2006. Measurements of molecular diffusion coefficients of carbon dioxide, methane, and propane in heavy oil under reservoir conditions. Energy Fuels 20 (6), 2509–2517. http://dx.doi.org/10.1021/ ef060080d. Tharanivasan, A.K., Yang, C., Gu, Y., 2004. Comparison of three different interface mass transfer models used in the experimental measurement of solvent diffusivity in heavy oil. J. Pet. Sci. Eng. 44 (3–4), 269–282. http://dx.doi.org/ 10.1016/j.petrol.2004.03.003. Upreti, S.R., Mehrotra, A.K., 2000. Experimental measurement of gas diffusivity in bitumen: results for carbon dioxide. Ind. Eng. Chem. Res. 39 (4), 1080–1087. http://dx.doi.org/10.1021/ie990635a. Upreti, S.R., Mehrotra, A.K., 2002. Diffusivity of CO2, CH4, C2H6 and N2 in Athabasca Bitumen. Can. J. Chem. Eng. 80 (1), 116–125. http://dx.doi.org/10.1002/ cjce.5450800112. Zhang, Y.P., Hyndman, C.L., Maini, B.B., 2000. Measurement of gas diffusivity in heavy oils. J. Pet. Sci. Eng. 25 (1–2), 37–47. http://dx.doi.org/10.1016/ S0920-4105(99)00031-5.