- Email: [email protected]

Thermal transients during nonisothermal fluid injection into oil reservoirs Ibrahim Kocabas * Chemical and Petroleum Engineering Department, UAE University, Al Ain, United Arab Emirates

Abstract During hot fluid injection into oil reservoirs, the importance of determining the temperature profile to estimate the thermal efficiency is well known. In addition, the resultant temperature distribution due to cold water injection into a hot reservoir during waterflooding may significantly influence the stress distribution in the reservoir. Owing to these facts, in this work we present a new two-dimensional analytical model for analyzing the thermal transients during nonisothermal fluid injections into oil reservoirs that may provide a better insight into the mechanisms of heat transfer in oil reservoirs. The new model has several distinguishing aspects. Primarily, it is the first analytical solution of an unsteady state twodimensional heat transport process in a laterally/vertically confined layer. In addition, both finite longitudinal and transverse heat dispersions have been accounted for in the model as well as the heat loss to the bounding layers. Thus, the model allows one to observe the roles of both boundary conditions and fluid mechanics controls simultaneously, a feature not possessed by the previous analytical models that assume either boundary conditions or fluid mechanics controls. Hydrodynamic heat dispersion concept has also been incorporated in the new model in order to account for the much larger temperature transition zones observed in the field than that would be caused by pure conduction. Finally, since the solution is free of numerical dispersion and nonphysical oscillations, especially in two-dimensional domains, it serves as a higher stepping stone for validation of numerical models. D 2003 Elsevier B.V. All rights reserved. Keywords: Thermal transients; Convection dispersion equation; Thermal dispersion; Nonisothermal flow; Heat transport

1. Introduction All modern secondary and/or enhanced oil recovery methods utilize the injection of a non-isothermal fluid into an oil reservoir. First of all, in case of a heat carrying fluid injection for a thermal recovery technique, determining the temperature profile and consequently estimating the thermal efficiency are of

* Tel.: +971-3-705-1546; fax: +971-3-762-4262. E-mail address: [email protected] (I. Kocabas). 0920-4105/$ - see front matter D 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.petrol.2003.12.006

great importance. In addition, the temperature dependence of processes such as acid reaction rate and the injected material properties such as polymer rheology makes an accurate estimation of temperature field a requirement in many EOR methods. Finally and most importantly, the geomechanical aspects of nonisothermal flow is of increasing concern in the oil industry such as the case of cold water injection into a hot reservoir, where the thermal and fluid pressure gradients could be used to estimate the net stress distribution, and hence, the possibility of fracture development.

134

I. Kocabas / Journal of Petroleum Science and Engineering 42 (2004) 133–144

A comprehensive modeling of the temperature transients in the porous media can best be achieved by the numerical models. However, in the implementation of numerical models there are some points to be understood and several critical issues to be resolved for the practicing engineer. First of all, gaining a good insight into the physical processes influencing the transport and accurately accounting for them in the numerical model is of vital importance. Analytical models are quite useful to provide the basic insight into the mechanisms of heat transport in porous media. Hence, a number of analytical models have been presented on the temperature distributions due to a purely convective flow in linear (Lauwerier, 1955), and radial (Malofeev, 1963) reservoirs and heat losses to the confining impermeable layers of infinite length. Rubinshtein (1959), on the other hand, presented an expression for the heating efficiency without deriving a temperature distribution for a convective flow in a radial layer with both longitudinal and transverse thermal conductivities in both the flowing and its confining layers. Then, Avdonin provided temperature distribution solutions that incorporate an additional finite longitudinal (Avdonin, 1964a) or transverse (Avdonin, 1964b) thermal conductivity term in the Lauwerier and Malofeev models. Later on, Spillette (1965) summarized the previous studies and presented a numerical solution to the problem of hot water injection into a horizontal layer, where assumptions made by Lauwerier and Advonin are relaxed by incorporating finite longitudinal and transverse thermal conductivities. For both Lauwerier and Malofeev models, Spillette (1965) has also reported the formulas for the thermal efficiency that is defined as the fraction of the total injected energy that remains in the productive layer. Further improvement of solutions to be able to model various cases such as finite confining layers separated by multiple fractures (Gringarten et al., 1975; Satman, 1988), modeling thermal injection backflow tests (Kocabas and Horne, 1990), applications to naturally fractured geothermal reservoirs (Cendejas and Rodriguez, 1994) and development of analogies to tracer transport (Ramirez et al., 1993; Kocabas and Islam, 2000a,b) has continued up to present. The issues of numerical dispersion, grid size and orientation effects and other behavioral errors such as nonphysical oscillations were recognized at very early

stages of numerical simulation of convective dispersive processes (Peacemean and Rachford, 1962). Later on, significant progress has been made on quantification (Lantz, 1971) and elimination (Leonard, 1981; Leonard et al., 1995; LeVeque, 1996) of numerical dispersion. However, numerical dispersion and nonphysical oscillations still continue to be major problems in numerical modeling of heat transport (Gupta, 1990; Kocabas and Margoub, 2000). Hence, for quantifying the accuracy of numerical models, comparison of the numerical model results with the available analytical models is of vital importance too. Therefore, development of more sophisticated analytical models serves invaluably for both gaining better insight into the transport mechanisms and numerical model validation purposes. In addition, the analytical models are valuable tools for studying the collective role of the parameters influencing the transport (Kocabas and Islam, 2000a,b).

2. Theoretical development The design and evaluation of the hot fluid injection into an oil reservoir depend on the ability to predict the thermal transients that occurs in the field during the project life. Such a prediction requires, in addition to the description of reservoir and thermal parameters, an accurate accounting of the thermal energy utilized in the process. Therefore, an understanding of the heat transfer mechanisms in porous media is essential to the development of the heat balance equation that accurately represents the processes involved. In this work we present a new two-dimensional analytical model for analyzing the thermal transients during nonisothermal fluid injections into oil reservoirs. The new model has several distinguishing aspects regarding both describing the processes and as well as accurately accounting for them in modeling. The detailed features of the new model are presented in the following. 2.1. Heat transfer mechanisms and thermal dispersion concept The proposed theory in this work for heat transfer mechanisms during a hot fluid injection into an oil reservoir has originated from the work of Sauty et al.

I. Kocabas / Journal of Petroleum Science and Engineering 42 (2004) 133–144

(1978). Based on this theory the heat transfer mechanisms during non-isothermal flow in porous in media may be divided into three major types. The first type is the energy transfer due to physical movement of the injected fluid. In other words, the average bulk movement of the injected fluid in the reservoir causes a forced heat convection. The second type of heat transfer is due to hydrodynamic dispersion. The velocity variations in both direction and magnitude in the porous media couple with the thermal conduction and give rise to a phenomenon called hydrodynamic thermal dispersion. The variations in the velocities are due to velocity profile in a single pore, velocity differing from one pore to another and tortuosity of the streamlines. While the variations in the velocity vector cause extra spreading of the transition zone between hot and cold fluids, the thermal conductivity makes this mixing irreversible. In other words, upon the reversal of the flow direction the transition zone does not become smaller, on the contrary, it becomes larger. Such a hydrodynamic dispersion occurs in both longitudinal and transverse directions and the corresponding thermal dispersion coefficients are defined as (Sauty et al., 1978): KL ¼ k þ aL qw cw uL

ð1Þ

Kt ¼ k þ at qw cw ut

ð2Þ

Additionally, heat is transferred by thermal conduction from the reservoir into the virtually impermeable confining layers. This assumption of virtual impermeability allows us to use conductivity in the confining layers rather than dispersion. Finally, the convective transfer between the injected fluid and the original reservoir fluids and solid material should be considered. Fortunately, the velocities in the reservoir are low enough to justify the assumption of instant thermal equilibrium. In fact, it is reported that the temperature equilibrium is reached in 1 s for 1 mm, in 1 min for 1 cm and in 2 h for 10 cm diameter grains (Sauty et al., 1982). Consequently, the concept of thermal dispersion is important to be included in the modeling, and hence, it forms a distinguishing feature of the new model.

135

2.2. Two dimensional analytical model The analytical solutions for temperature distributions presented in earlier works (Lauwerier, 1955; Malofeev, 1963; Avdonin, 1964a,b; Gringarten et al., 1975; Satman, 1988; Kocabas and Horne, 1990; Ramirez et al., 1993; Cendejas and Rodriguez, 1994; Kocabas and Islam, 2000a,b) are all onedimensional solutions. In addition they can account for only one type of dominant heat transfer mechanism, and hence, can be divided into two types. The first types of solutions are those that assume the heat losses to impermeable bounding layer, and hence, the boundary conditions dominate. The first and the most well-known of these solutions is due to Lauwerier (1955). The Lauwerier model assumes a linear, one-dimensional incompressible flow in a homogeneous layer, constant physical properties and saturations, infinite thermal conductivity in transverse direction and zero longitudinal thermal conductivity in the reservoir, finite and equal transverse thermal conductivities and zero longitudinal thermal conductivities in the bounding layers, and a constant injection temperature. The solution to the Lauwerier problem for the radial flow geometry was presented by Malofeev (1963). The two other solutions regarding the temperature distributions due to hot liquid injection were presented by Avdonin (1964a,b). His first work included a constant longitudinal conductivity for both linear and radial geometries. All other assumptions were identical to those made by Lauwerier. Later, he presented solutions for a finite transverse conductivity in the reservoir but the longitudinal one was assumed to be negligible. Such solutions are applicable to relatively thin productive layers. Then, slight modifications to Avdonin’s first model were presented by Gringarten et al. (1975), Satman (1988), and Kocabas and Islam (2000a,b). The second type of models assumes that the reservoir is thick enough to neglect the heat losses to the impermeable bounding layer. However, the reservoir is assumed to consist of two superposed continua namely a highly permeable fracture network and a relatively small permeability matrix blocks. Such models are considered to be controlled by fluid mechanics. The major heat transfer mechanism is the one that occurs between the fractures and matrix

136

I. Kocabas / Journal of Petroleum Science and Engineering 42 (2004) 133–144

blocks. Such types of models are based on tracer transport models (Ramirez et al., 1993) or geothermal models (Cendejas and Rodriguez, 1994) and are applicable to relatively thick and, also, densely fractured reservoirs. In the new analytical model presented here, the assumptions in Lauwerier’s model have been relaxed to include both a finite longitudinal dispersion (usually neglected in earlier models) and a transverse thermal dispersion in the productive layer, making the solution two-dimensional. In order to build a conceptual model, first we assume a linear flowing unit confined by two impermeable layers. We also assume an incompressible fluid flow with constant linear flow velocities in x and z directions leading to constant longitudinal and transverse thermal dispersion coefficients. In addition, the assumption that the velocity in z direction is negligible compared to that in x direction allows us to neglect the convective transport in z direction and yet preserve the thermal transverse dispersion term in the governing equations. Finally, we assume constant and equal heat conduction coefficients in the bounding layers. 2.3. Governing equations and development of solution The two-dimensional governing differential equations of the heat transport during nonisothermal fluid injection into an oil layer where no phase change and instantaneous equilibrium of injected and resident liquids are assumed are as follows: Bðqr cr T Þ þ j ðqf cf ðq=AÞT Þ j KjT ¼ 0 Bt

ð3Þ

and Bðqm cm Tm Þ j km jTm ¼ 0; Bt

ð4Þ

where q=A ¼ /u

Based on the assumptions of the conceptual model described in the previous section, the governing equations become: qr cr

BT BT B2 T B2 T þ /uqf cf KL 2 Kt 2 ¼ 0; Bt Bx Bx Bz

qm cm

BTm B2 Tm km ¼ 0: Bt Bz2

ð5Þ

qr cr ¼ /qf cf þ ð1 /Þqs cs

ð6Þ

qm cm ¼ /m qw cw þ ð1 /m Þqs cs

ð7Þ

ð9Þ

The initial and boundary conditions are specified as follows: T ¼ Tm ¼ T0 T ¼ Ti

at

t¼0

at

ð10Þ

x¼0

ð11Þ

Both mediums are assumed to be semi-infinite: T !0 Tm ! 0

as

xD ! l

ð12Þ

as zD ! l

ð13Þ

Then the two equations are coupled through the equality of temperatures and fluxes at the boundary. T ¼ Tm

Kt

at

z¼b

BT BTm ¼ km Bz BZ

at

ð14Þ

z¼b

ð15Þ

The mathematical formulation is simplified using following dimensionless variables: TD ¼

T0 T T0 Ti

/qf cf ux ; KL /qf cf uz zD ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ KL Kt xD ¼

and qf cf ¼ qw cw Sw þ qo co So

ð8Þ

h1 ¼

qm cm Kt ; qr cr km

and TmD ¼

tD ¼

T0 Tm T0 Ti

ð/qf cf uÞ2 t ; qr cr KL

ð16Þ

and ð17Þ

h2 ¼

pﬃﬃﬃﬃﬃ km ; and h ¼ h2 h1 Kt

ð18Þ

I. Kocabas / Journal of Petroleum Science and Engineering 42 (2004) 133–144

In terms of the dimensionless variables the governing equations and initial and boundary conditions reduce to: BTD BTD B2 TD B2 TD þ ¼ 0; BtD BxD Bx2D Bz2D

h1

BTm B2 Tm ¼ 0; BtD Bz2D

TD ¼ TmD ¼ 0 TD ¼ 1

xD ! l

as

TmD ! 0

tD ¼ 0

zD ! l

as

at

zD ¼ zDb

BTD BTmD ¼ h2 BzD BzD

at

zD ¼ zDb

ð25Þ

ð26Þ

ð19Þ

ð20Þ

xD ¼ 0

at

TD ! 0

at

TD ¼ TmD

137

In order to derive the solution, we have first performed a dependent variable transformation for dimensionless temperatures in both regions: TD ¼ vexpðxD =2Þ;

ð27Þ

TmD ¼ vm expðxD =2Þ:

ð28Þ

ð21Þ ð22Þ ð23Þ ð24Þ

Then we have applied the Laplace and Fourier sine transforms to the equations of v and vm. Finally, applying the initial and boundary conditions we have obtained the following solution in the Laplace-Fourier space (Kocabas, 2001).

" # pﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ hs exp s þ k zD þ exp s þ k zD 2x 1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ v¨ ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ psðs þ kÞ s þ k þ hs exp s þ k zDb s þ k hs exp s þ k zDb

ð29Þ

where k = x2 + a. To simplify the derivations we have assumed that the parameter h has a unit value. This assumption is made in order to obtain a real space solution. After manipulating Eq. (29) considerably and using tables of Laplace transforms together with the convolution theorem, we have obtained the following Fourier sine

space solution for the dimensionless temperature in the oil layer:

1 xD tD pﬃﬃﬃﬃ erfc TD ¼ 2 2 tD

xD þ t D pﬃﬃﬃﬃ þ expðxD Þerfc ð30Þ TDsinv2 2 tD

where l Z X

TDsinv2 ¼

n¼0

0

tD

expððð2n þ 1ÞzDb zD Þ2 =4sÞ þ expððð2n þ 1ÞzDb þ zD Þ2 =4sÞ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ Fsinv ds p sðtD sÞ

ð31Þ

and Fsinv ¼

Z 0

l

2

2

x þa x þa ðt ðt exp sÞ I sÞ sin xxdx: D D nþ1=2 2 2 ðx2 þ aÞ1=2 1

ð32Þ

138

I. Kocabas / Journal of Petroleum Science and Engineering 42 (2004) 133–144

To carry out the above inverse Fourier sine transform (Eq. (32)), we have utilized the following relation which is selected as the most convenient one among the available Modified Bessel Function In + 1/2 representations: Inþ1=2 ðyÞ ¼

rﬃﬃﬃﬃﬃ 2y 1 p 2y 9 8 X n ðn þ mÞ! > > m > > ð2yÞ expðyÞ > > = < m!Cðn m þ 1Þ m¼0 : n X > > ðn þ mÞ! n > > > ð2yÞm expðyÞ > ; : ð1Þ m!Cðn m þ 1Þ m¼0

Then we have performed a term-by-term inversion using the tables of Fourier sine and cosine integral transforms, and the theory of convolution. For example when n = 0, the first approximation to the term TDinv2 in the solution becomes: 0 TDinv2 ¼

Z

expððzDb zD Þ2 =4sÞ þ expððzDb þ zD Þ2 =4sÞ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ p sðtD sÞ 0

9 8 > > erfc s pxﬃﬃﬃD expðx Þerfc s þpxﬃﬃﬃD > > D = 1< 2 s 2 s

ds 2> t x t þ x > > erfc D pﬃﬃﬃﬃﬃD þ expðx Þerfc D pﬃﬃﬃﬃﬃD > ; : D 2 tD 2 tD tD

ð34Þ

ð33Þ For n = 1, the second order approximation to the term TDinv2 in the solution is obtained as: Z

expðð3zDb zD Þ2 =4sÞ þ expðð3zDb þ zD Þ2 =4sÞ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ p sðtD sÞ 0

9 8 1 s xD s xD s þ xD 1 s þ xD > > > > p ﬃﬃ ﬃ p ﬃﬃ ﬃ þ þ Þerfc expðx erfc D = < 2 t s tD s 2 2 s 2 s D

ds > t xD 1 t D xD t D þ xD 1 t D þ xD > > ; : D pﬃﬃﬃﬃ þ pﬃﬃﬃﬃ > þ erfc þ expðxD Þerfc 2 2 tD s 2 tD tD s 2 tD

1 0 ¼ TDinv2 þ TDsinv2

tD

The equations for n = 2 are quite complex and do not lend themselves to a practical solution easily. However, as discussed in the average temperatures section, the accuracy of this second order approximation is quite adequate for the computational purposes, and for most of the time we have utilized this form to compute the temperature profiles.

3. Two-dimensional temperature profiles in the oil layer The two-dimensional temperature distributions give us a general idea about the temperature surface in the oil layer. Such information may be used to design and interpret the fluid injection operations into oil reservoirs. It may also be used to justify/discard a one-dimensional model. Figs. 1, 2 and 3 show temperature gradients developed along the transport direction as well as across half of the injection plane of the reservoir, where the other half would be the exact mirror image.

ð35Þ

In the two dimensional temperature distributions, the following features are of great importance, namely, the extent of the transition region in both lateral/vertical and longitudinal directions and magnitude of gradients in each direction. These features may be inferred from the number of grids with similar shape, size and orientation. The square shaped small grids indicate regions of constant temperatures at or close to the injection temperature. The rectangular shaped grids indicate regions of temperature variations, the larger the aspect ratio of a rectangle the greater the temperature gradient. Finally, the orientation of rectangles coincides with the direction of the gradient. Fig. 1 shows that the temperature surface for a large Peclet number and a small L/b ratio consists of three distinct regions. The small size square grids occupy the largest area indicting that most of the reservoir is heated to the injection fluid temperature. Then, there exists a narrow region of high lateral/vertical temperature gradients next to the confining layer boundary. Finally, there is a region of high temperature gradients

I. Kocabas / Journal of Petroleum Science and Engineering 42 (2004) 133–144

139

Fig. 1. Dimensionless temperature distributions across half of the injection plane, Pe = 100, KL/Kt = 10, L/b = 10.

along the flow direction. From these observations, we expect a high heating efficiency. Fig. 2 shows the effect of increasing dispersivity (smaller Peclet numbers) on the temperature distribution. There are also three distinguishable regions in Fig. 2. Here the increasing dispersivity reduces the portion of the reservoir heated close to the injection fluid temperature. The two transition regions near the confining layer boundary and along the flow direction occupy relatively larger portions than that in Fig. 1. These regions are characterized by large rectangles with small aspect ratios. Thus, while the transition

zone gets larger the temperature gradients become less pronounced. Fig. 3 shows the effect of a large L/b ratio representing a long thin reservoir. In this case, the transition zone develops only along the flow direction extending over the whole flow distance. The smoothed transverse gradients indicate the total dominance of heat losses through the confining layer boundary leading to a low heating efficiency. More importantly, however, is that two-dimensional temperature distribution surfaces serve a higher stepping stone for validating numerical models. Nu-

Fig. 2. Dimensionless temperature distribution across half of the injection plane, Pe = 100, KL/Kt = 10, L/b = 10.

140

I. Kocabas / Journal of Petroleum Science and Engineering 42 (2004) 133–144

Fig. 3. Dimensionless temperature distributionsacross half of the injection plane, Pe = 100, KL/Kt = 10, L/b = 100.

merical models are known to suffer at varying degrees from errors such as numerical dispersion, and grid size orientation effects and other behavioral errors such as nonphysical oscillations. In fact these errors may become quite difficult to evaluate for multidimensional models (Kocabas and Margoub, 2000). Therefore, this is the first time a two-dimensional analytical model is presented for confined reservoirs, and hence, it gives us the opportunity to match two dimensional error free temperature surfaces serving to better evaluate some of the errors involved in numerical models.

4. Average temperature profile in the oil layer Obtaining an expression of average temperature profile is important for determining the heating efficiency, which is the most important parameter of thermal recovery methods. Particularly for the present model, the average temperature serves as an accurate one-dimensional solution. In addition, we can infer the collective roles of parameters and geometric dimensions on the heating efficiency. The average temperature across the half width of the reservoir (which is also equal to the average temperature distribution across the whole width) is found using the following relation: 1 T¯ D ¼ zDb

Z

ZDb

TD ðzD Þ dzD 0

ð36Þ

Applying this relation to the solution, one obtains:

1 xD t D pﬃﬃﬃﬃ erfc T¯ D ¼ 2 2 tD

xD þ tD pﬃﬃﬃﬃ þ expðxD Þerfc ð37Þ T¯ Dsinv2 2 tD where l Z tD 1 X T¯ Dsin v2 ¼ zDb n¼0 0

nzDb ðn þ 1ÞzDb pﬃﬃﬃ erfc pﬃﬃﬃ erfc s s pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ p sðtD sÞ

Fsinv ds:

ð38Þ

Similar to the derivation of the two dimensional temperatures, the following approximate solutions are derived for the average temperatures along the transport direction. For n = 0, the first approximation to the term T¯Dinv2 in the solution is: 0 T¯ Dsinv2

zDb erf pﬃﬃﬃ 1 s pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ zDb 0 p sðtD sÞ

9 8 s xD s þ xD > > > pﬃﬃﬃ expðxD Þerfc pﬃﬃﬃ > erfc > > = < 2 s 2 s 1 ds

> 2> > > t x t þx > ; : erfc D pﬃﬃﬃﬃﬃD þ expðxD Þerfc D pﬃﬃﬃﬃﬃD > 2 tD 2 tD Z

tD

ð39Þ

I. Kocabas / Journal of Petroleum Science and Engineering 42 (2004) 133–144

Fig. 4. The new solution for L/b = 20 and various KL/Kt values.

141

Fig. 6. Comparison of Lauwerier, Avdonin and 2D models for KL/ Kt = 1 and various L/b values.

For n = 1, the second order approximation to the solution is:

1 T¯ Dsin v2

zDb 2zDb Z tD erfc pﬃﬃﬃ erfc pﬃﬃﬃ 1 s s 0 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ¼ T¯ Dsinv2 þ zDb 0 pðtD sÞ

9 8 1 s xD s xD s þ xD 1 s þ xD > > > > p ﬃﬃ ﬃ p ﬃﬃ ﬃ = < 2 þ t s erfc 2 s t s þ 2 expðxD Þerfc 2 s D D

ds

> 1 t xD t D xD tD þ xD 1 t D þ xD > > > ; : þ D pﬃﬃﬃﬃ þ pﬃ þ expðxD Þerfc erfc 2 2 tD s 2 tD tD s 2 tD

The average temperature equations are first used to investigate the accuracy of the first and second order approximations of the full solution (n = 0 and n = 1, respectively). Figs. 4 and 5 show the temperature profiles for a range of parameters for the cases of n = 0 and n = 1 where curves of the first and second approximations coincide except the lowest pairs of curves in both figures. The concurrence of the first and second order approximations indicates that the approximate solution for those cases is highly accu-

Fig. 5. The new solution for KL/Kt = 1 and various L/b values.

ð40Þ

rate. The lowest pairs of curves in both figures correspond to the combinations of a large L/b with a low KL/Kt. In the case of a large L/b, however, one would expect that the Lauwerier and Avdonin assumptions to hold, and hence, those solutions to approach the full solution. In fact, Figs. 6 and 7 show the comparison of the second order approximation with the Lauwerier and Avdonin solutions. The lower pair of curves in Figs. 6 and 7 corresponding to L/

Fig. 7. Comparison of Lauwerier, Avdonin and 2D models for L/ b = 1 and various KL/Kt values.

142

I. Kocabas / Journal of Petroleum Science and Engineering 42 (2004) 133–144

b = 30 and 20 respectively and KL/Kt = 1, show that the differences between the Avdonin and the two dimensional models are negligible. Thus, as Avdonin solution approaches the full solution, its proximity to the second order approximation as shown in Figs. 6 and 7, further builds the confidence on the accuracy of the second order approximation. Therefore, considering that curves in Figs. 4, 5 and 6 correspond to a wide range of practical values, we can be confident that the second order approximation is sufficiently accurate for computational purposes. In Figs. 6 and 7, a comparison of the Luwerier and Avdonin models with the two dimensional model indicates that for a small KL/Kt value, as the aspect ratio (i.e. large L/b) increases the Lauwerier and Avdonin solutions approximate the two dimensional solution fairly closely. On the other hand, as aspect ratio gets smaller or KL/Kt increases even in the case of fairly large L/b (i.e. 20) values, use of the two dimensional solution becomes a necessity for accurate temperature profile estimations. The average temperature distributions for a wide range of values of the aspect ratio, L/b and the Peclet number, Pe, are given in Figs. 8 and 9, respectively. Fig. 8 shows the effect of L/b ratio on the average temperature distribution for a large Peclet number. Qualitatively, heating efficiency is proportional to the area under the curve of average temperature. Thus, one can see that as L/b ratio increases, heating efficiency decreases rapidly. If one were to define the effective heating of the reservoir as increasing the reservoir temperature to 0.5 TD with 2PV injection, Fig. 8 shows that even for a small dispersivity case, L/ b should be at least 20. Naturally, increasing disper-

Fig. 9. Effect of Peclet number on average temperature distribution, for PV = 0.9, K/Kt = 10 and L/b = 20.

sivity will increase the requirement of L/b ratio or PV injection. Fig. 9 shows the effect of dispersivity on the average temperature distributions for an L/b ratio of 10 and indicates that as the dispersivity increases the transition zone will get larger, reducing the effectiveness of the reservoir heating. Finally, one can infer that as L/b ratio becomes greater or Peclet number gets smaller, the transition zone along the flow direction increases. However, as for their collective role on increasing the transition zone along the flow direction, L/b ratio dominates, smoothing the temperature gradient across the reservoir width and spreading the transition zone over the whole flow distance. A major advantage of the present solution is that it provides a more accurate thermal efficiency estimate than the one-dimensional models that use infinite transverse conductivity. However, neither the development of an analytical expression nor a numerical integration of the average temperature of the new model profile for the thermal efficiency calculations was attempted here and should be dealt with elsewhere.

5. Conclusions

Fig. 8. Effect of L/b ratio on average temperature distribution, for PV = 0.9, Pe = 100, and K/Kt = 10.

A new two-dimensional analytical solution of the convection – dispersion equation based on a constant linear velocity and constant longitudinal and transverse dispersion coefficients is presented. The thermal transients during nonisothermal fluid injection into a linear confined oil reservoir have been studied

I. Kocabas / Journal of Petroleum Science and Engineering 42 (2004) 133–144

through these solutions. The model presented for studying the thermal transients has several significant features such as ability to generate error free twodimensional temperature surfaces and employing the thermal dispersion concept. Furthermore, it serves to observe the roles of both fluid mechanics through dispersion parameters and boundary conditions via heat loss terms at the confining boundaries. Thus, it has the potential to be used in all modern nonisothermal liquid injection processes for design and interpretation purposes as well as identifying the numerical dispersion errors involved in two dimensions for comprehensive numerical simulators. Nomenclature a a constant equal to 0.25 A area open to flow c heat capacity k thermal conductivity K thermal dispersion In + 1/2 Modified Bessel Function of order (n + 1/2) L length of the reservoir m dummy variable (integer) n dummy variable (integer) Pe longitudinal Peclet number, Pe=xDxf PV pore volume of water injected, PV=ut/L=tD/Pe q injection rate into the reservoir s Laplace transform variable S saturation t time variable of the transport equations T temperature T¯D dimensionless average temperature u interstitial flow velocity x space variable along the flow direction y dummy variable z space variable perpendicular to the flow direction a dispersivity / porosity v transformed temperature variable v¨ Laplace and Fourier Sine transformed v v˙ Laplace transformed v g space convolution variable s time convolution variable pﬃﬃﬃﬃﬃ h dimensionless parameter equal to h2 h1 h1 dimensionless parameter equal to (qmcmKt)/ (qr cr km) h2 dimensionless parameter equal to km/Kt

C q k x j

143

Gamma function density parameter equal to (x2 + a) Fourier sine transform variable gradient operator

Subscripts b boundary D dimensionless i inlet L longitudinal m confining layer 0 initial o oil r reservoir t transverse s solid grains sinv inverse Fourier sine w water

Acknowledgements Funding of this research was provided by the Abu Dhabi National Oil Company (ADNOC). Many fruitful discussions with Dr. Hassan Hejase of UAE University have significantly facilitated our work. Both contributions are gratefully acknowledged.

References Avdonin, N.A., 1964a. Some formulas for calculating the temperature field of a stratum subject to thermal injection. Neft Gaz 7 (3), 37 – 43 (in Russian). Avdonin, N.A., 1964b. On the different methods of calculating the temperature fields of a stratum subject during thermal injection. Neft Gaz 7 (8), 39 – 44 (in Russian). Cendejas, F.A., Rodriguez, H.R., 1994. Heat transfer processes during low or high enthalpy fluid injection into naturally fractured reservoirs. Proc. Nineteenth Workshop on Geothermal Engineering. Stanford University, Stanford, CA. Gringarten, A.C., Witherspoon, P.A., Ohnishi, Y., 1975. Theory of heat extraction from hot dry rock. J. Geophys. Res. 80 (8), 1120 – 1124. Gupta, A.D., 1990. Accurate resolution of physical dispersion in the multidimensional numerical modeling of miscible and chemical displacement. SPE Reserv. Eng., Nov., 581 – 588. Kocabas, I., 2001. Modeling of heat/tracer transport in oil reservoirs. Externally Funded Project Report IV submitted to Exploration and Production Division ADNOC, UAE Univ., Al Ain, UAE.

144

I. Kocabas / Journal of Petroleum Science and Engineering 42 (2004) 133–144

Kocabas, I., Horne, R.N., 1990. A new method of forecasting the thermal breakthrough time during reinjection in geothermal reservoirs. Proc. Sixteenth Workshop on Geothermal Engineering. Stanford University, Stanford, CA. Kocabas, I., Islam, M.R., 2000a. Concentration and temperature transients in heterogeneous porous media, Part I: Linear transport. J. Pet. Sci. Eng. 26, 211 – 220. Kocabas, I., Islam, M.R., 2000b. Concentration and temperature transients in heterogeneous porous media: Part I. Radial transport. J. Pet. Sci. Eng. 26, 221 – 233. Kocabas, I., Margoub, A., 2000. Improved numerical schemes for transport equations in oil reservoirs. Externally Funded Project Report III submitted to Exploration and Production Division ADNOC, UAE Univ., Al Ain, UAE. Lantz, R.B., 1971. Quantitative evaluation of numerical diffusion (truncation error). Soc. Pet. Eng. J., 315 – 320. Lauwerier, H.A., 1955. The transport of heat in an oil layer caused by the injection of hot fluid. Appl. Sci. Res. 5 (2 – 3), 145 – 150. Leonard, B.P., 1981. A survey of finite differences with upwinding for numerical modelling of the incompressible convective diffusion equation. In: Taylor, C., Morgan, K. (Eds.), Computational Techniques in Transient and Turbulent Flow, vol. 2. Pineridge, Swansea, pp. 1 – 35. Leonard, B.P., MacVean, M.K., Lock, A.P., 1995. The flux integral method for multidimensional convection and diffusion. Appl. Math. Model. 19, 333 – 342. LeVeque, R.J., 1996. High-resolution conservative algorithms for

advection in incompressible flow. SIAM J. Numer. Anal. 33 (2), 627 – 665. Malofeev, G.E., 1963. Calculation of the temperature distribution in a formation when pumping hot fluid into a well. Neft Gaz 9 (2), 31 – 45 (in Russian). Peacemean, D.W., Rachford, H.H., 1962. Numerical calculation of multidimensional miscible displacement. Soc. Pet. Eng. J., 327 – 338. Ramirez, J., Samaniego, F.V., Rivera, J.R, Rodriguez, F., 1993. Tracer flow in naturally fractured reservoirs. SPE Paper #25900, Proceedings, Rocky Mountain Regional/Low Permeability Reservoirs Symposium held in Denver, CO, USA, April 12 – 14. Rubinshtein, L.I., 1959. The total heat losses in injection of a hot liquid into a stratum. Neft Gaz 2 (9), 110 – 116 (in Russian). Satman, A., 1988. Solutions of heat- and fluid flow problems in naturally fractured reservoirs: Part 1-Heat flow problems. SPE Prod. Eng., 463 – 466 (Nov.). Sauty, J.P, Gringarten, A., Landel, P.A., 1978. The effect of thermal dispersion on injection of hot water in aquifers. Proc. The Second Invitational Well Testing Symposium. Div. of Geothermal Energy, U.S. Dept. of Energy, Berkeley, CA, pp. 122 – 131. Sauty, J.P, Gringarten, A.C, Menjoz, A., Landel, P.A., 1982. Sensible energy storage in aquifers: 1. Theoretical study. Water Resour. Res. 18 (2), 245 – 252. Spillette, A.G., 1965. Heat transfer during hot fluid injection into an oil reservoir. In: Dietz, D.N. (Ed.), Thermal Recovery Techniques. SPE Reprint Series, vol. 10, pp. 21 – 26.