Application of iterated Laplace transformation to tracer transients in heterogeneous porous media

Journal of the Franklin Institute 348 (2011) 1339–1362 www.elsevier.com/locate/jfranklin

Application of iterated Laplace transformation to tracer transients in heterogeneous porous media Ibrahim Kocabas Chemical Engineering Department, American University of Sharjah, P.O. Box 26666, Sharjah, UAE Received 31 May 2009; received in revised form 4 April 2010; accepted 6 April 2010 Available online 22 April 2010

Abstract The Laplace transformation technique has been widely applied to modeling of tracer transport in oil and geothermal reservoirs, and in groundwater aquifers. However, mathematical models of many flow and transport problems could only be obtained as Laplace space solutions, and hence, their computations had to involve a numerical inversion technique. In this work, we employ the iterated Laplace transformation technique to develop novel closed form solutions to the tracer transport models in heterogeneous media. Two types of configurations have been considered: tracer transport in single fracture located in low-permeability matrix and tracer transport in a double porosity medium consisting of flowing and dead-end pore systems. In addition, both linear and radial flow geometries have been considered for both configurations. Applications of iterated Laplace transform technique to these four types of models are presented as fundamental examples and their numerical results were used as benchmarking for the numerical inversion results from Stehfest and Dubner and Abate algorithms. As the technique is quite versatile, we expect that the method should gain widespread acceptance to develop solutions to a wide range of problems in flow and transport in porous media and improve the application of nonlinear regression technique to these solutions. This work has achieved four important objectives: first, two novel Laplace transform relations that are useful in tracer studies are presented. Second, the present work serves to verify/invalidate the results of numerical inversion algorithms. In addition, it provides better insight into tracer transport

Tel.: þ971 6 515 2973; fax: þ971 6 515 2979.

E-mail address: [email protected] 0016-0032/$32.00 & 2010 The Franklin Institute. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.jfranklin.2010.04.002

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

1340

mechanisms. Finally, it serves as a powerful tool of design and interpretation of tracer tests. All four objectives are illustrated in this work. & 2010 The Franklin Institute. Published by Elsevier Ltd. All rights reserved. Keywords: Iterated Laplace transform; Tracer transport; Tracer transients; Convection dispersion equation solutions; Solute transport

1. Introduction This work involves the application of iterated Laplace transform technique to solve a number of important flow and transport problems in porous media. The following is a brief summary of the iterated Laplace transform technique.

1.1. Iterated Laplace transform The double Laplace transform is defined as [1] Z 1 Z 1 f ðs; pÞ ¼ est1 ept2 F ðt1 ; t2 Þdt2 dt1 0

ð1Þ

0

The iterated Laplace transform is a special case of double Laplace transform in which the Laplace space variables s and p are equal. Z 1 Z 1 f ðs; sÞ ¼ est1 est2 F ðt1 ; t2 Þdt2 dt1 ð2Þ 0

0

The most useful property of this iterated Laplace transform is that it results in generalized convolution principle as follows Z t Z 1 st f ðs; sÞ ¼ e F ðt; ttÞdt ð3Þ 0

0

The proof of the generalized convolution principle as a result of iterated Laplace transform is presented by Sneddon [1]. One of the earliest applications of iterated Laplace transform appeared in a work by De Smedt and Weirenga [2]. The same year, application of iterative Laplace transform theory has also been detailed in a notable work by De Smedt et al. [3]. De Smedt et al. [4] have revisited the iterated Laplace transform in 2006 where they have applied it to develop an instantaneous injection solution for solute transport in rivers. According to Eq. (3), inverting f(s,s) requires application of the two Laplace transform inversion successively Z t 1 1 1 L f ðs; sÞ ¼ L2 fL1 f ðs; sÞg ¼ F ðt; ttÞdt ð4Þ 0

Carrying out the first iteration of inversion with respect to some arbitrarily selected s terms results in an intermediate function L1 1 f ðs; sÞ ¼ gðt; sÞ

ð5Þ

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

1341

Repeating the inversion with respect to the remaining s terms and the corresponding real space variable is taken as ðttÞ yields L1 2 gðt; sÞ ¼ F ðt; ttÞ

ð6Þ

The fact that the inversion has been iterated requires a convolution, and hence, result becomes Z t L1 f ðs; sÞ ¼ F ðt; ttÞdt ð7Þ 0

The simplicity and versatility of the technique makes iterated Laplace transform a powerful tool of solving quite complex problems. Also note that using iterated Laplace transform technique, we can derive some novel Laplace inversion relations, which are useful for solving some important tracer transport problems in porous media. Eqs. (8) and (11) are novel Laplace inversion relations that are presented for the first time in this work.   Z t pffiffi bt F ðt=aÞ 1 1 f ðaðs þ b sÞÞ ¼ dt ð8Þ L erfc pffiffiffiffiffiffiffiffi s a 2 tt 0 This relation can be derived using iterated Laplace transforms as follows. The first iteration of inversion is carried out to obtain t pffi pffiffi 1 1 L1 f ðaðs þ b sÞÞ ¼ gðt; sÞ ¼ eb st F ð9Þ 1 s as a The second iteration of inversion yields   pffi t bt F ðt=aÞ 1 1 bt s e F ð Þ ¼ erfc pffiffiffiffiffiffiffiffi L2 as a a 2 tt

ð10Þ

The convolution after the second iteration yields Eq. (8). Using iterated Laplace transform, we may also derive Eq. (11) pffiffi     Z t  pffiffiffiffiffiffiffiffi sÞ 2 tt t2 t 1 L f ðs þ 3=2 ¼ f ðtÞdt ð11Þ pffiffiffi exp  t erfc pffiffiffiffiffiffiffiffi s 4ðttÞ p 2 tt 0 1.2. Tracer transport in heterogenous porous media Heterogeneous porous media as manifested in petroleum and geothermal reservoirs, or groundwater aquifers can exist in various configurations such as a high conductivity flow path located in a low-permeability matrix such as a fault zone in a reservoir or a fractures network and a rock matrix system, as in naturally fractured reservoirs. Because of a high permeability and storativity contrast between the fractures and the matrix, such systems are referred to as heterogeneous porous media. Regardless of existing configuration, models of heterogeneous porous media have mostly been employed in design and interpretation work only by using a numerical Laplace inversion algorithm such as Stehfest [7] or Dubner and Abate [8] algorithms. These two methods are considered to be current as they have been improved several times over the years regarding their convergence and efficiency and have been widely employed in flow and transport problems in porous media. These arguments have been supported with ample evidence in the proceeding Section 2.9.

1342

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

In the following sections, we employ iterated the Laplace transforms to develop novel solutions to tracer transport models of these two major configurations of heterogeneous porous media and use them to validate/invalidate the numerical inversion results.

2. Mathematical models of tracer transport The transport of tracers in the above-mentioned configurations of heterogeneous porous media have usually been modeled using two different approaches namely, considering a boundary conditions dominated transport and considering a fluid mechanics dominated transport. The boundary conditions dominated transport approach assumes that the system is composed of a fluid flowing layer located in a matrix where the fluid is assumed to be virtually stagnant. A single fracture located in a low-permeability matrix and a long high-permeability layer confined by two low-permeability layers are examples of such systems. During the flow of tracer containing fluid in the fracture the tracer is exchanged between these two elements of the systems through the fracture matrix interface. The fluid mechanics dominated approach, on the other hand, assumes that the system is composed of two superposed continua in space, in other words a double porosity reservoir. The first continuum consists of a flowing zone and the second continuum consists of a stagnant fluid containing zone. The porous medium consisting of a dense fracture network located in a low-permeability matrix, such as a naturally fractured reservoir or a porous medium with dead-end pores, is a physical example of fluid mechanics dominated transport approach. In modeling the tracer transport, two flow geometries may be considered, namely linear flow and radial flow. Both boundary and fluid flow dominated systems in each flow geometry will be considered in the following sections.

2.1. Boundary conditions dominated systems in linear flow geometry The mathematical model of tracer transport in a boundary conditions dominated system is constructed with the flowing assumptions: (i) the porous system consists of a single fracture located in low-permeability matrix of infinite extent, horizontal and constant thickness; (ii) a single-phase steady-state incompressible traced fluid flow takes place in the fracture, and the fluid in the matrix is virtually immobile; (iii) a convective transport and a longitudinal dispersion occurs in the fracture, and only a diffusive transport takes place in the matrix in a perpendicular direction to flow; (iv) the transverse dispersion is of infinite magnitude instantaneously equalizing the concentrations across the fracture and (v) the exchange of tracer between the fracture and the matrix occurs at the matrix/fracture interface and is proportional to the transverse concentration gradients in the matrix. Then, the dimensionless governing equations consist of two coupled linear differential equations representing the transport in fracture and matrix [9,10] @CD @CD @2 CD @C mD 1=2 þ  þ l ¼0 ð12Þ @tD @xD @yD yD ¼0 @x2D

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

@Cm @2 Cm  2 ¼0 @tD @yD

1343

ð13Þ

Tracer can be injected in three modes, namely: continuous injection, slug injection and spike injection. For all three types of modes the initial and the boundary conditions are common except the inlet boundary condition. The common conditions are specified as follows: CD ¼ CmD ¼ 0 at tD ¼ 0

ð14Þ

Both media are assumed to be semi-infinite CD -1 as xD -1

ð15Þ

CmD -1 as yD -1

ð16Þ

The continuity of concentrations at the fracture/matrix interface requires CD ¼ CmD at yD ¼ 0

ð17Þ

Dimensionless variables are defined as follows: xD ¼

l¼

ux ; D

tD ¼

fm Dm D ; b2 u2

u2 t ; D

Pe ¼

uy yD ¼ pffiffiffiffiffiffiffiffiffiffiffi DDm uL ; D

Pt ¼

b2 u fDm L

ð18Þ

ð19Þ

Note that l=1/(PePt) 2.2. Continuous injection solution The continuous injection solution involves a flux of constant concentration fluid entering the high permeability streak, and hence, the inlet condition is set as CD ¼ 1 at xD ¼ 0

ð20Þ

where the dimensionless concentrations are defined as CD ¼

CC0 Cm C0 ; and CmD ¼ Ci C0 Ci C0

ð21Þ

The resulting solution from this inlet boundary condition is referred to as flux–flux concentration solution. The continuous injection solution is most useful in interpretation of laboratory experiments. Solving these two coupled dimensionless linear differential equations with the boundary conditions given by Eqs. (14)–(20), result in the Laplace space solution given by Eq. (22) (see Appendix A)  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   pffiffiffiffiffi xD 1 C D ¼ exp 1 1 þ 4ðs þ lsÞ ð22Þ s 2

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

1344

Application of an iterated Laplace transform inversion to Laplace space solution will yield the real space closed form solution (also see Appendix A) 8  9 xD ðxD tÞ2 > > > > p ffiffiffiffiffiffiffi exp > > > Z tD > = < 2 pt3 4t ! p ffiffi ffi dt ð23Þ CD ðt; tD tÞ ¼ lt > 0 > > > > > erfc pffiffiffiffiffiffiffiffiffiffiffi > > ; : 2 tD t Eq. (22) has been derived in many transport studies [9,10] and mostly been computed using numerical Laplace transform inversion techniques such as either Stehfest [7], Dubner Abate or Crump algorithms [8]. The importance of Eq. (23) is that it shows the influence of each transport mechanism distinctly. 2.3. Instantaneous injection solution In field experiments, a continuous injection is hardly ever practiced, instead an instantaneous injection is most commonly applied. For an instantaneous injection of a mass of tracer m into the fluid flux entering the system, the inlet boundary condition is specified as CD ¼ dðtD Þ at xD ¼ 0

ð24Þ

where the dimensionless flux concentration is defined as CD ¼

C ðm=QÞðu2 =DÞ

ð25Þ

Solving the related equations with this inlet boundary condition yields the following Laplace space solution:   qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffi xD C D ¼ exp ð1 1 þ 4ðs þ lsÞÞ ð26Þ 2 9 8   xD ðxD tÞ2 > > > > p ffiffiffiffiffiffiffi exp > > > Z tD > = < 2 pt3 4t p ffiffi ffi   dt CD ¼ 2 lt lt > 0 > > > > > qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp  > > : 4ðtD tÞ ; 2 pðtD tÞ3

ð27Þ

The instantaneous injection is not only convenient because it requires only small amount of tracer to be injected, but also useful as its return profile provides better resolution for the influence of heterogeneity for the nonlinear regression analysis of tracer return profiles. 2.4. Fluid mechanics dominated systems in linear geometry The mathematical model of tracer transport in a fluid mechanics dominated system is constructed with the following assumptions: (i) the reservoir consists of a dense fracture network located in low-permeability matrix of infinite extent, horizontal and constant thickness, and hence the pore system can be divided into a flowing zone and a dead-end pores; (ii) a single-phase steady-state incompressible fluid flow takes place in the fractures

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

1345

and the fluid in the matrix is virtually immobile; (iii) a convective transport and a longitudinal dispersion occurs in the fracture, and only a diffusive transport takes place in the matrix in a perpendicular direction to flow and (iv) while the flow in the system occurs only via fractures (i.e. flowing pores), the matrix (i.e dead-end pores) acts as a uniformly distributed source feeding into the fractures. Then, the dimensionless governing transport equations consist of two coupled linear differential equations [10] o

@CD @CmD @CD @2 CD þ ð1oÞ þ  2 ¼0 @tD @tD @xD @xD

ð28Þ

@CmD ¼ lðCD CmD Þ ð29Þ @t The initial and the boundary conditions of this system are also given by those of the single fracture heterogeneity system, namely Eqs. (14)–(20), except Eq. (16). As Eq. (29) does not require a boundary condition, Eq. (16) does not apply for this system. The dimensionless variables are the same as defined by Eqs. (18)–(21), with the exception of the parameter l given by Eq. (30), which is the ratio of the Damkohler number defined by Eq (31) to the Peclet number KD NDa l¼ 2 ¼ ð30Þ u Pe ð1oÞ

NDa ¼

KL u

ð31Þ

2.5. Continuous injection solution Solving these two coupled dimensionless linear differential equations with the boundary conditions given by Eqs. (14)–(21), will result in the Laplace space solution given by Eq. (32) (see Appendix B)  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi xD  1 C D ¼ exp ð1 1 þ 4ðos þ f ðsÞÞÞ ð32Þ s 2 where f ðsÞ ¼ l

l2 ð1oÞs þ l

ð33Þ

Application of an iterated Laplace transform inversion to Laplace space solution will yield the real space closed form solution (also see Appendix B) CD ¼ CD1 þ CD2 where

CD1

ð34Þ

! 3 xD ytD =o pffiffiffiffiffiffiffiffiffiffiffi þ 7 ð1yÞ erfc 6 exp 2 7 2 tD =o 16 6 !7 ¼ 6 7   26 þ ytD =o 7 5 4 exp xD ð1 þ yÞ erfc xD p ffiffiffiffiffiffiffiffiffiffiffi 2 2 tD =o 2

x

D



ð35Þ

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

1346

y¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 4l

  xDpffiffioffi ðxD t=oÞ2 pffiffiffiffiffiffiffi exp G2 dt 4t=o 0 2 pt3 8

pffiffiffiffiffiffiffiffiffiffiffi pffiffiffi  9 exp ð lt=o ZÞ2 > > > > > > > Z lðtD tÞ=ð1oÞ > = < expð2pffiffiffiffiffiffiffiffiffiffiffiffiffi Zlt=oÞ dZ G2 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi > > 0 > > I1 ð2 Zlt=oÞ > > > > pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; : Zo=ðltÞ Z

ð36Þ

tD

CD2¼

ð37Þ

ð38Þ

Note that equivalent forms of this solution has been developed earlier by many authors [2–4,11–14,16] as early as the work of Lindstrom and Narashiman, in 1973 [11]. Also note that equivalency of these solutions have been shown by Goltz et al. [15]. In this case using this or previously developed solutions become a matter of preference. We content that the present solution should be preferred as it is developed using a fundamental inversion technique, and computationally it is as efficient as any other method. 2.6. Instantaneous injection solution The instantaneous injection solution with appropriate initial and boundary conditions for the fluid mechanics dominated system is given by the following components:     xDpffiffioffi ðxD tD =oÞ2 ltD CD1 ¼ qffiffiffiffiffiffiffiffi exp exp  ð39Þ 4tD =o w 2 pt3D Z

tD

CD2 ¼ 0

  xDpffiffioffi ðxD t=oÞ2 pffiffiffiffiffiffiffi exp G2 dt 4t=o 2 pt3

ð40Þ

where

 pffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2  9 8 > > = < exp  lt=o lðtD tÞ=ð1oÞ pffiffiffi G2 ¼

pffiffiffi pffiffiffi Z > > ; : exp 2 Z I1 2 Z tD t

ð41Þ

and Z¼

l2 tðtD tÞ oð1oÞ

ð42Þ

2.7. Boundary conditions dominated systems in radial geometry For radial transport several considerations are of great importance. First, the transport in diverging flow in radial geometry may be developed using two assumptions: (1) a constant dispersion coefficient, D; (2) a velocity-dependent dispersion coefficient and (3) the dispersivity may be considered as scale dependent, which may also lead to a constant dispersion coefficient in radially diverging flow. All cases have been thoroughly treated by

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

1347

Kocabas and Islam [17] where solutions are presented in Laplace space. Only the constant dispersion case will be treated in this work. Second, despite some authors have used the boundary conditions dominated system for radial geometry for solute/tracer transport, Chen [18], we content that the boundary conditions dominated system is more suitable for heat transport than tracer transport. In other words, sensible energy injection into horizontal oil reservoirs confined by impermeable layers [19], or into aquifers [20,21]is possible to model as a boundary dominated system. As the existence of a horizontal fracture in porous media specially in oil reservoirs is unlikely, the tracer transport in a horizontal fracture will not be considered in this work. However, as temperature may be considered as a physical tracer, the boundary conditions dominated system will also be treated in terms of concentrations as dependent variables in the present work. This option is adapted for the purpose of reducing number of variables in the present work. A second issue related to solutions of convection dispersion equations is the necessity to distinguish between the resident and flux concentrations concepts [10,17,22,23]. Kocabas and Islam [17] have classified the radial convection dispersion equation solutions according to the two concentration variables. In that work all solutions were presented in Laplace space. Kocabas and Islam [17] have also stated that the suitable solutions for practical applications are those called the flux concentration solutions. In the following the real space solutions for flux concentrations obtained by iterated Laplace transform technique are presented for the first time in this work. The boundary condition dominated system is visualized as follows: a horizontal and thin permeable layer is confined by two impermeable layers of infinite extent. A convective and dispersive transport takes place radially in the permeable layer. The concentrations across the vertical thickness of the permeable layer are assumed to be quickly equalized by transverse dispersion. Thus, an infinite transverse dispersion is assumed in the permeable layer. The tracer/solute is assumed to be lost into the confining layers from the permeable layer via diffusion. The two coupled linear differential equations govern the transport in the flowing layer and the confining layers. The complete set of differential equations in terms of resident and flux concentrations is presented by Kocabas and Islam [17]. In this work we will consider only the dimensionless flux concentration equations, which are given as follows: !   @2 C D 2 þ Pe @CD @CD pffiffiffi@CmD   þ l ð43Þ 2rD @rD @tD @yD @r2D yD ¼0

@2 CmD @CmD ¼ @tD @y2D

ð44Þ

Eq. (43) assumes that there is a transient diffusion between the permeable layer and the confining layers. For a flux of constant concentration fluid entering a line source well, the initial and the boundary conditions in the system are stated as follows: CD ¼ CmD ¼ 0 at tD ¼ 0

ð45Þ

CD -0 as rD -1

ð46Þ

CD ¼ 1 as rD -0

ð47Þ

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

1348

CmD -0 as yD -1

ð48Þ

CD m ¼ CD at yD ¼ 0

ð49Þ

Dimensionless variables are defined as follows: CD ¼

CC0 ; Ci C0

CmD ¼

Cm C0 Ci C0

r Dt ; tD ¼ 2 ; rw rw rffiffiffiffiffiffiffi Dm y Dm r2w and l ¼ yD ¼ D rw Dh2 rD ¼

ð50Þ ð51Þ

ð52Þ

Solving these two coupled linear differential equations in Laplace space with the appropriate initial and boundary conditions will result ðr2 ðs þ fs ÞÞv=2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ðs þ fs ÞrD Kv C D ¼ D v1 ð53Þ s2 GðvÞ The first iteration of inversion yields (see Appendix C)  v  2 t1 r2D r expðfs tÞ CD ¼ exp  D s GðvÞ 4t 4t

ð54Þ

The second inversion is carried out and the result is convoluted to yield (Appendix C pffiffiffi !  2 Z tD 1  2 v t rD rD lt CD ¼ ð55Þ exp  erfc pffiffiffiffiffiffiffiffiffiffiffi dt 2 t GðvÞ 4t 4t D t 0 2.8. Fluid mechanics dominated systems in radial geometry The dimensionless governing equations for constant dispersion case in terms of flux concentration variables consist of two coupled linear differential equations representing the transport in flowing and immobile portions of the porous medium   @2 CD 2 þ Pe @CD @CD @CmD  ¼o þ ð1oÞ ð56Þ 2 2r @r @t @tD @rD D D D ð1oÞ

@CmD ¼ lðCD CmD Þ @tD

ð57Þ

Eq. (57) assumes that there is a pseudo steady-state interporosity flow between the flowing and the immobile portions of the porous medium. Note that for a flux of constant concentration fluid entering a line source well, the initial and the boundary conditions in the system are also given by Eqs. (45)–(47). Dimensionless variables are also defined by Eqs. (50) and (51), only the parameter l is defined differently as follows: l¼

Km r2w D

ð58Þ

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

1349

Solving these two coupled linear differential equations in Laplace space with the appropriate boundary condition will result in (see Appendix D) qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðr2 os þ r2 fs Þv=2 Kv C D ¼ D v1 D r2D os þ r2D fs ð59Þ s2 GðvÞ The first iteration inversion is given by (see Appendix D)     v  t1 or2D fs t or2D CD ðt; sÞ ¼ exp  exp  o GðvÞ 4t 4t Let us define g1 ðtÞ ¼

    1 4t v or2D exp  2v GðvÞ or2D 4t

ð60Þ

ð61Þ

and the second inversion iteration is the same as that in Appendix B Eq. (B4). Let us also define g2 ðt; tD tÞ by Eq. (B4). The full solution should read   lt Z tD exp  g1 ðtÞ o CD ðt; tD tÞ ¼ ð62Þ 0 ð1 þ g2 ðt; tD tÞÞ dt After the same manipulations used in Appendix B, g2 ðt; tD tÞ is by Eq. (B14). Let us define lt o Then, the full solution should read Z tD expðeÞg ðtÞ 1 CD ðt; tD tÞ ¼ ð1 þ g2 ðt; tD tÞÞ dt 0 e¼

ð63Þ

ð64Þ

where 8 pffiffi pffiffiffi 9 expðð e ZÞ2 Þ > > > > > Z lðtD tÞ=ð1oÞ > ffi = < expð2pffiffiffiffi g2 ðt; tD tÞ eZÞ dZ ¼ pffiffiffiffiffi > > expðeÞ I1 ð2 ZeÞ 0 > > > > pffiffiffiffiffiffiffi ; : Z=e

ð65Þ

2.9. Validation of numerical inversion of Laplace transforms Numerical inversion of Laplace transforms have been extensively used to solve equations of fluid flow and solute transport in porous media. Among a large number of techniques, which are reviewed by Davies and Brian [24], two methods have emerged as the most commonly employed techniques of numerical inversion. The first method is based on the Fourier series approximations of the functions and involves Laplace space function evaluations using complex variables. Dubner and Abate [8] have first introduced this method using Fourier Cosine series and the method has been improved by Durbin [25], Crump [26] and DeHoog et al. [27], who all used both Cosine and Sine series

1350

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

approximations together with an acceleration algorithm. While Durbin [25] had obtained limited success, Crump [26] and later DeHoog et al. [27] obtained significant success in the application of acceleration algorithms. A different line of research has also been followed on the computation of inverse Laplace transforms using Fourier series, which involved Fast Fourier Transform (FFT) technique, stated by Ichikawa [28], Inoue et al. [29] and Yonemoto et al. [30]. Since then, the Fourier series approach has been the preferred inversion algorithm of flow and transport problems in porous media [31,10,32]. Note, however, that the fundamental features of Dubner and Abate algorithm such as accuracy and capability or weakness to handle functions with steep slopes is still preserved in the improved algorithms of Durbin [25], Crump [26] and DeHoog [27]. The original Dubner Abate algorithm requires a large number of function evaluations using complex variables; however, it also allows one to specify the relative error of the numerically inverted function value. In other words, the user can control the accuracy of the inverted value. The Stehfest algorithm is based on a simple probabilistic concept, which does not require any complex variables function evaluations. The Stehfest algorithm has worked with outstanding accuracy for fluid flow problems in porous media such as the line source solution [33], finite wellbore solution in an infinite medium and cumulative production solution for a constant pressure inner boundary solution. Its computational efficiency is also exceptional as it involves only a small number of Laplace space function evaluations. These findings have lead to its employment with certain confidence in some other problems whose real space solutions were not developed and hence it has been the most commonly employed method of numerical Laplace inversion up to today for fluid flow problems [33–38]. Nevertheless, the exact limits of applicability of Stehfest algorithm were not determined as it has no mechanism of quantifying either the absolute error or the relative error in the value of numerically inverted function values. Since both methods are still being used for tracer transport problems [39–44,32], a comparison of the their efficiencies would serve as a guide for future applications in similar problems. Figs. 1–3 show the tracer return profiles for single fracture (boundary conditions dominated system) model over a range of longitudinal and transverse Peclet numbers Pe and Pt. Note that transverse Peclet number Pt embedded in l is defined by Eq. (19). While both Stehfest and Dubner and Abate methods yield quite accurate results for low longitudinal Peclet number (i.e. Pe=30, Fig. 1), matching the exact solution perfectly, for high Peclet numbers (i.e. Pe=100, Figs. 2 and 3), Stehfest algorithm differs significantly from exact solution. These results show that one can use the Dubner and Abate algorithm confidently for the boundary condition dominated systems. Figs. 4 and 5 show the instantaneous injection tracer return profiles for dead-end porosity system (fluid mechanic dominated system), where Peclet number is fixed to 100 and l is varied. The flowing fraction of the pore volume is fixed for all cases at o=0.1. Since we have a high Peclet number, Stehfest algorithm yields significantly inaccurate results in both cases. The Dubner and Abate algorithm reproduces the exact solutions with great accuracy except at one single point in Fig. 5. In fact the same result also exists in the results for Fig. 4, but since the inaccurate point is beyond 2 on the tD axis, it is unseen in the figure. We have verified that if the coefficients required for Dubner and Abate algorithm are computed using a greater time value than the actual measurement time, then the inaccurate point totally is displaced from the sampling domain of the return profiles. These results also increase one’s confidence in using the Dubner and Abate algorithm in more complex problems where developing an exact analytical solution is not feasible.

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

1351

2.10. Comparison of limiting forms of solutions For fluid mechanics dominated system, four different approximate solutions can be developed for the limiting values of tD, o, l and Peclet number. The first limiting solution is the classical convection dispersion equation given by Eq. (66), which is obtained for o=1.     1 xD tD xD þ tD pffiffiffiffiffi þ expðxD Þ erfc pffiffiffiffiffi CD ¼ erfc ð66Þ 2 2 tD 2 tD The second limiting solution is given by Eq. (35). For short times, which correspond to large s, or for a very small l, which makes the second part of Eq. (33) negligible, Eq. (33) reduces tof ðsÞ ¼ l, which in turn reduces the solution to Eq. (35). The third approximation is obtained by assuming l to be zero, which reduces the value of y=1, and the solution for any value of o and for all times reduces to Eq. (67) ! 3 2 xD tD =o pffiffiffiffiffiffiffiffiffiffiffi þ 7 6 erfc 7 2 tD =o 16 7 6 ! ð67Þ CD ¼ 6 7 7 26 x þ t =o D D 5 4 expðxD Þerfc pffiffiffiffiffiffiffiffiffiffiffi 2 tD =o Finally, assuming a long flow system, which makes Pe to approach infinity, the solution again reduces to the classical convection dispersion Eq. (66)) with an apparent Peclet number equal to [2]   1 1 1 ð1oÞ2 2 ¼ þ ð1oÞ NDa ¼ 1þ ð68Þ Pa Pe Pe l Note that the apparent Peclet number is smaller than the Peclet number and hence there is a greater dispersion accounted by the apparent Peclet number. These limiting forms of the solutions are indeed quite useful to gain some insight into the transport mechanisms as exhibited by Figs. 6 and 7. Fig. 6 shows that for small l, the exchange of tracer between the interconnected and dead-end pores is of limited quantity

Fig. 1. Single fracture model profiles for Pe=30 and l=100.

1352

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

Fig. 2. Single fracture model profiles for Pe=100 and l=30.

Fig. 3. Single fracture model profiles for Pe=100 and l=100.

Fig. 4. Dead-end pores model profiles for Pe=100, l=0.001 and o=0.9.

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

1353

Fig. 5. Dead-end pores model profiles for Pe=100, l=0.1 and o=0.9.

and slow in time, and hence, the solutions, Eqs. (35) and (44), assuming the transport is taking place within the interconnected pores is matching the exact solution fairly well until late times. Only at late times, the classical solution with an apparent dispersion coefficient Pa (Eq. (66) with Eq. (68)) is matching the exact solution. Fig. 7 shows that for large l values, the complete mixing model, Eq. (43) with Pa, is in perfect agreement with the exact solution. Eq. (44), which assumes l=0 is differing at all times, and Eq. (35) yields a zero value for all times. This means it has no contribution to the full solution for large l values, and hence, exact solution is solely determined by Eq. (37). While differing from the exact solution significantly in Fig. 6, Eq. (43) with Pe matches the exact solution perfectly in Fig. 7. This is due to the fact that Pa=1.01Pe for o=0.9 and l=0.1 or l=1.0 values the latter l value happens to be the case in Fig.7. Even though the limiting solutions are instructive to identify dominant transport mechanisms for certain times, none of them reproduces the complete return profiles at any time. Therefore, the fairly simple exact analytical solution should be used for interpretation and forecasting purposes. 2.11. Insight into the transport mechanisms More insight is gained on the transport mechanisms using the exact solutions for instantaneous injection case. Fig. 8 shows that for large l values (i.e. lZ0.1), a fast complete mixing is taking place between the interconnected and dead-end pores resulting in a symmetrical tracer return profile with a center of gravity arrival time approximately equal to one total pore volume injection. Fig. 8 also shows that for small l values (i.e. lr0.001), a slow and significantly incomplete mixing is taking place between the interconnected and dead-end pores also resulting in a symmetrical tracer return profile with a center of gravity arrival time approximately equal to one interconnected pore volume injection. Note that since o=0.8, the interconnected pore volume is equal to 0.8 total pore volume. For other l values between those limits, we expect a substantially asymmetrical return profile with a breakthrough close to that of insignificant mixing and significant tailing at later times, as exhibited by the curve of l=0.01 in Fig.8. The tailing is absent in both cases of large and small l values.

1354

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

Qualitatively, one can reach the conclusion that the exchange between interconnected and dead-end pore would delay the transport of tracer through the system. However, the advantage of the model is that the model provides us the quantitative values of parameter ranges of limiting conditions, of delay and of tailing in time domain in the return profile and hence a better insight into the transport mechanisms.

2.12. Design and interpretation of tracer tests The field tracer tests usually involve an instantaneous injection of a small mass of tracer into the fluid flux entering the system. Thus, to interpret return profiles of such systems, an instantaneous injection solution should be used, which is simply the time derivative of the continuous injection solutions developed earlier. Usually, an optimization software is employed in the interpretation of the tracer return profiles. The computerized matching of measured tracer data from a field experiment with the model-generated values involves a nonlinear regression analysis or nonlinear curve fitting. The instantaneous injection is not only convenient because it requires only small amount of tracer to be injected, but also useful as its return profile provides better resolution for the influence of heterogeneity. Therefore, it is expected to obtain a faster convergence in nonlinear regression analysis compared to the continuous injection solutions. For the design of tracer test, there are two important issues besides several other concerns: amount of tracer to be injected and sample collection interval in time. The model can run and set the limits of concentration values expected in an observation point. More importantly, the parametric sensitivity analysis shows it is most important to collect samples related to the transition region of return profiles. Therefore, sampling frequency in time is determined and preparations are made to collect samples at the specified time values using results of model-generated return profiles. For the interpretation of tracer tests, a nonlinear curve fitting technique is used to match the field data by the model-generated profiles. From the match, the likely values of the parameters of the system are inferred. Fig. 9 shows a typical nonlinear match result for a three parameter system using a simulated field data.

Fig. 6. Simplified model profiles for Pe=100, and l=0.001 and o=09.

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

1355

Fig. 7. Simplified model profiles for Pe=100, and l=1.0 and o=0.9.

Fig. 8. Exact solution dead-end pores model profiles for Pe=100 and o=0.8.

3. Conclusions A number of novel analytical solutions are presented for modeling tracer transients in heterogeneous porous media. These solutions are developed using the iterated Laplace transform technique. Both boundary conditions dominated systems and fluid mechanics dominated systems have been considered in the modeling. In both cases, the new analytical solutions have been used to serve three important goals, namely: verifying numerical inversion algorithms, gaining insight into tracer transport mechanisms and functioning as design and interpretation tools of actual tracer tests. The exact solutions are compared with the two most commonly used numerical Laplace inversion algorithms, namely Stehfest and Dubner and Abate algorithms. For both systems, while Dubner and Abate algorithm matched quite accurately the exact solutions

1356

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

Fig. 9. Curve fitting of dead-end pores model profiles for Pe=100, l=1.0 and o=0.9.

for all practical ranges of Peclet numbers, Stehfest algorithm matched the exact solution only for highly dispersive systems. The limiting solutions are investigated to specify the applicability range of limiting solutions. Also, and limiting values of parameters are investigated to identify the expected return profile shapes. For design and interpretation purposes, parametric sensitivity analysis and nonlinear regression should be used, respectively. For optimization purposes, instantaneous injection solutions should be used as they facilitate the convergence in nonlinear optimization procedure. One can confidently use either the exact solution or the Dubner and Abate numerical inversion algorithm in design and interpretation procedures, provided that the required time adjustment is made in computing inversion coefficients. Finally, the iterated Laplace transform is a versatile method of developing analytical solutions, and hence, we expect it to find many areas of application in solving fluid flow and transport problems in porous media. Appendix A. Derivation of solutions to single fracture model in linear geometry The Laplace transform application reduces Eqs. (12) and (13) to @2 C D @C D 1=2 @C mD  þl sC D ¼ 0 2 @xD @yD yD ¼0 @xD @2 C mD sC mD ¼ 0 @y2D

ðA1Þ

ðA2Þ

Solving Eq. (A2) with boundary conditions given by Eqs. (14), (16) and (17), one can obtain its derivative at fracture matrix interface pffiffi @C mD ¼  sC D ðA3Þ @yD yD ¼0

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

1357

Substituting Eq. (A3) into Eq. (A1) and solving with initial and boundary conditions of Eqs. (14), (16) and (20), will yield the Laplace space solution given by Eq. (22). Using the following relations [6]:   pffiffi 1 1 L1 expð sÞ ¼ pffiffiffiffiffiffiffi exp  ðA4Þ 4t 2 pt3     1 b t L f ðas þ bÞ ¼ exp  t F a a a 1

application of first iteration of Laplace transform inversion to Eq. (22) will result  2 x  x xD D D pffiffiffiffiffiffiffi exp CD ðt; sÞ ¼ exp 2 2 pt3 4t p ffiffiffiffi ffi 1 expðð1=4 þ lsÞtÞ s Let us define  2 x t  xD xD D pffiffiffiffiffiffiffi exp g1 ðtÞ ¼ exp  3 2 4 2 pt 4t pffiffiffipffiffi 1 g2 ðt; tD tÞ ¼ L1 expðt l sÞ s

ðA5Þ

ðA6Þ

ðA7Þ

ðA8Þ

Using [6] 1

L

  pffiffi 1 k expðk sÞ ¼ erfc pffiffiffiffiffi s 2 tD

ðA9Þ

we obtain

pffiffiffi ! l g2 ðt; tD tÞ ¼ erfc t pffiffiffiffiffiffiffiffiffiffiffi 2 tD t

The final results may be written as Z tD g1 ðtÞg2 ðt; tD tÞ dt Cðt; tD tÞ ¼

ðA10Þ

ðA11Þ

0

Substitution of g1 and g2 will yield the desired solution given by Eq. (23). Appendix B. Derivation of solutions to dead-end porosity model in linear geometry The Laplace transform reduces Eqs. (28) and (29) to @2 C D @C D  osC D ð1oÞsC mD ¼ 0 @xD @x2D

ðB1Þ

ð1oÞslC mD ¼ ðC D C mD Þ

ðB2Þ

Solving Eq. (B2) for CmD and substituting it into Eq. (B1), and then, solving Eq. (B1) with boundary conditions of Eqs. (14) and (15) results in Eq. (32).

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

1358

To develop the real space solution, we use Eqs. (A6) and (A7), to obtain the first iteration of Laplace inversion       pffiffiffiffi o ðxD t=oÞ2 1 lt l2 t=o ðB3Þ exp  CD ðt; sÞ ¼ xD pffiffiffiffiffiffiffi exp exp s o ð1oÞs þ l 4t=o 2 pt3 For the second iteration of Laplace inversion, we define g1 by Eq. (A7) but redefine g2 as     1 l2 t=o exp g2 ðt; tD tÞ ¼ L1 1 ðB4Þ s ð1oÞs þ l The final solution should read 8 9  t Z tD < exp  g1 ðtÞ = lo dt CD ðt; tD tÞ ¼ 0 : ð1 þ g2 ðt; tD tÞÞ ;

ðB5Þ

To facilitate for the reader to follow the derivations, let us define the following:   Z tD lt exp  g1 ðtÞdt CD1 ¼ o 0 G2 ¼ expð Z CD2¼

lt Þg2 ðt; tD tÞ o

ðB6Þ

ðB7Þ

tD

g1 ðtÞG2 dt

ðB8Þ

0

Thus, Eq. (34) will represent the final solution. Substituting Eq. (A7) into Eq. (B6) and performing a change of variables leads to ffi pffiffiffiffi Z pffiffiffi tD o  ðB9Þ CD1 ¼ xD pffiffiffi

pexp au2 þ ub2 du 0 where a¼

  1 1 1þ ; 4o l

b¼

x2D f 4

ðB10Þ

Eq. (B9) is integrated using following relation [6]:   9d 8 b > > > expð2abÞerf au þ þ > > pffiffiffi > Z d = < u p au2 ðb=u2 Þ   e du ¼ b > 4a > c > > > ; : expð2abÞerf au u >

ðB11Þ

c

Carrying out the integration in Eq. (B9) will result in Eq. (35). The following two relations and Eq. (A5) are employed to obtain the function g2 [5,6] 1 pffiffi L1 ðexpð1=sÞ1Þ ¼ pffiffi tI1 ð2 tÞ Z t 1 f ðsÞ L ¼ F ðtÞdt s 0

ðB12Þ

ðB13Þ

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

After some algebraic manipulations we obtain 9 8 > Z ðtD tÞ=ð1oÞl > < expðZÞ pffiffiffiffiffiffiffiffiffiffiffiffi = I1 ð2 Zt=ðolÞÞ dZ g2 ðt; tD tÞ ¼ > > 0 ; : pffiffiffiffiffiffiffiffiffiffiffiffiffi Zlo=t

1359

ðB14Þ

Substituting Eq. (B14) into (B7) yields Eq. (38) and subsequently substitution of Eq. (A7) into Eq. (B8) yields Eq. (37). Appendix C. Derivation of solutions to single fracture model in radial geometry The Laplace transform of Eqs.(43) and (44) becomes   pffiffiffi@C mD @2 C D 2 þ Pe @C D  sC D þ l Þ 2 2rD @rD @yD yD ¼0 @rD @2 C mD sC mD ¼ 0 @y2D

ðC1Þ

ðC2Þ

Solving Eq. (C2) with boundary conditions given by Eqs. (48) and (49) and substituting the result into Eq. (C1) yields   pffiffiffiffiffi @2 C D 2 þ Pe @C D  ðs þ lsÞC D ¼ 0 ðC3Þ 2rD @rD @r2D The solution of Eq. (C3) is given by [17] pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C D ¼ ArvD Kv ð ðs þ fs ÞrD Þ þ BrvD Iv ð ðs þ fs ÞrD Þ

C4Þ

where Pe þ1 4 pffiffiffiffiffi fs ¼ ls v¼

Using the following relation from Tang and Peaceman and Eq. (47) lim pffiffiffi v pffiffiffi  brD Þ Kv ð brD ¼ 2v1 GðvÞ rD -0

ðC5Þ ðC6Þ

ðC7Þ

We determine the second constant A as A¼

bv=2 1 v1 s 2 GðvÞ

ðC8Þ

where b ¼ s þ fs

ðC9Þ

Setting B=0 and substituting A into Eq. (C4), we obtain the solution in Laplace space given by Eq. (53). Eq. (53) is rewritten as follows: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðr2 s þ r2 fs ÞÞv=2 Kv C D ¼ D v1 D r2D s þ r2D fs ðC10Þ s2 GðvÞ

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

1360

Using the following inversion relation:   pffiffi 1 1 1 v=2 L s K v ð sÞ ¼ exp  4t ð2tÞvþ1

ðC11Þ

and the relation given by Eq. (A5), we obtain the first iteration inversion as Eq. (55). Since pffiffiffiffiffi expðfs tÞ expð lstÞ ¼ ðC12Þ s s The inversion of Eq. (11) is given by Eq. (A10) as second iteration of inversion and applying the convolution to full solution after second iteration inversion results in Eq. (56). Appendix D. Derivation of solutions to dead-end porosity model in radial geometry The Laplace transform of Eqs. (56) and (57) and their solutions in Laplace space is given by (D1) pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ðos þ fs ÞrD þ BrvD Iv ðos þ fs ÞrD C D ¼ ArvD Kv ðD1Þ where fs ¼

ð1oÞls sð1oÞ þ l

Applying the initial and the boundary conditions to Eq. (D1) would result in qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðr2 os þ r2 fs Þv=2 Kv C D ¼ D v1 D r2D os þ r2D fs s2 GðvÞ

ðD2Þ

ðD3Þ

Using Eqs. (C11) and (A5), we obtain the first iteration of inversion as Eq. (60). Note Eq. (60) may be written as  v   t1 or2D or2D CD ðt; sÞ ¼ exp  GðvÞ 4t 4t     ðD4Þ 2 1 lt l t=o exp  exp s o ð1oÞs þ l The second iteration of inversion directly follows from Eq. (B4). The final convoluted solution is obtained as Eq. (62). References [1] I.A. Sneddon, in: The Use of integral Transforms, MCGRaw Hill, New York, 1972. [2] F. De Smedt, P.J. Wierenga, A generalized solution for solute flow in soils with mobile and immobile water, Water Resour. Res. 15 (5) (1981) 1137–1141. [3] F. De Smedt, P.J. Wierenga, A. Van der beken, in: Theoretical and Experimental Study of Solute Movement Through Porous Media with Mobile and Immobile Water, Vrije Universtait Brussel, Brussel, 1981. [4] F. De Smedt, W. Brevis, P. Debels, Analytical solution for solute transport resulting from instantaneous injection in streams with transient storage, J. Hydrol. 315 (25-39) (2005). [5] H. Bateman, A. Erdelyi, Tables of Integral Transforms, vol. 1, McGraw Hill, 1954. [6] M. Abramowitz, I.A. Stegun, in: Handbook of Mathematical Functions, 10th ed., Dover Publications, New York, 1972. [7] H. Stehfest, Algorithm 368: numerical inversion of Laplace transform D-5, Commun. ACM 13 (1) (1970) 47–49.

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

1361

[8] H. Dubner, J. Abate, Numerical inversion of Laplace transforms by relating them to the Finite Fourier Cosine Transform, J. ACM 13 (1) (1968) 115–123. [9] I. Kocabas, Modeling tracer flow in oil reservoirs containing high permeability streaks, SPE Paper 81429, Proceedings of MEOS, June 9–12, Bahrain, 2003. [10] A.C. Correa, K.K. Pande, H.J. Ramey Jr., W.E., Brigham, Prediction and interpretation of miscible displacement using a transverse matrix dispersion model, SPE 16704, September 1987. [11] F.T. Lindstrom, M.N.L. Narasimhan, Mathematical theory of a kinetic model for dispersion of previously distributed chemicals in a sorbing porous medium, SIAM J. Appl. Math. 24 (1973) 496–510. [12] C.L. Carnahan, J.S. Remer, Nonequilibrium and equilibrium sorption with a linear sorption isotherm during mass transport through an infinite porous medium: some analytical solutions, J. Hydrol. 73 (1984) 227–258. [13] M.N. Goltz, P.V. Roberts, Three-dimensional solutions for solute transport in an infinite medium with mobile and immobile zones, Water Resour. Res. 22 (7) (1984) 1139–1148. [14] J.J.A. van Kooten, A method to solve the advection dispersion equation with a kinetic adsorption isotherm, Adv. Water Resour. 19 (4) (1995) 193–206. [15] J. Huang, M.N. Goltz, P.V. Roberts, Comment on ‘Analytical solution for solute transport resulting from instantaneous injection in streams with transient storage’ by F. De Smedt, W. Brevis, P. Debels, Journal of Hydrology 315(2005) 25–39, Journal of Hydrology 330(2006)759–760. [16] F. De Smedt, B. Brevis, P. Debels, Reply to comment: Analytical solution for solute transport resulting from instantaneous injection in streams with transient storage by F. De Smedt, W. Brevis, P. Debels, Journal of Hydrology 315(2005) 25–39, Journal of Hydrology 330(2006)761–762. [17] I. Kocabas, M.R. Islam, Concentration and temperature transients in heterogeneous porous media, part I: linear transport, J. Pet. Sci. Eng. 26 (2000) 211–222. [18] S.C. Chen, Solutions for radionuclide transport from an injection well into a single fracture in a porous formation, Water Resour. Res. 22 (4) (1986) 508–518. [19] A.G. Spillette, Heat transfer during hot fluid injection into an oil reservoir, in: D.N. Dietz (Ed.), Thermal Recovery Techniques, 10, SPE Reprint Series, 1965, pp. 21–26. [20] J.P. Sauty, A.C. Gringarten, A. Menjoz, P.A. Landel, Sensible energy storage in aquifers 1. Theoretical study, Water Resour. Res. 18 (2) (1982) 245–252. [21] J.P. Sauty, A. Gringarten, P.A. Landel,. The effect of thermal dispersion on injection of hot water in aquifers. Proceedings of the 2nd Invitational Well Testing Symposium, Division of Geothermal Energy, U.S. Department of Energy, Berkeley, CA, 1978, pp. 122–131. [22] A. Kreft, A. Zuber, On the physical meaning of the dispersion equation and its solutions for different boundary conditions, Chem. Eng. Sci. 33 (1978) 1471–1480. [23] A.J. Valocchi, Effect of radial flow on deviations from local equilibrium during sorbing solute transport through homogeneous soils, Water Resour. Res. 22 (3) (1986) 1673–1701. [24] B. Davies, M. Brian, Numerical inversion of Laplace transforms: a survey and comparison of methods, J. Comput. Phys. 33 (1979) 1–32. [25] F. Durbin, Numerical inversion of Laplace transforms: an efficient improvement to Dunner and Abate’s method, Comput. J. 17 (1974) 371–376. [26] K.S. Crump, Numerical inversion of Laplace transforms using a Fourier series approximation, J. Assoc. Comp. Mach. 23 (1976) 80–96. [27] F.R. DeHoog, J.H. Knight, A.N. Stokes, An improved method for numerical inversion of Laplace transforms, SIAM J. Sci. Stat. Comput. 3 (3) (1982) 357–366. [28] S. Ichikawa, Numerical processing of the two dimensional inverse Laplace transform and its application to equation of heat conduction, Mem. Fac. Eng. Kyoto Univ. 44 (2) (1982) 231–248. [29] H. Inoue, M. Kamibayashi, K. Kishimoto, T. Shibuya, T. Koizumi, Numerical Laplace transformation and inversion using fast Fourier transform, JSME Int. J. Ser. 1: Solid Mech. Strength Mater. 35 (3) (1992)319–324. [30] A. Yonemoto, T. Hisakado, K. Okumura, An improvement of convergence of FFT-based numerical inversion of Laplace transforms, Proc. IEEE Int. Symp. Circuits Syst. 5 (2002) V/769–V/772. [31] J.A. Barker, Laplace transform solutions for solute transport in fissured aquifers, Adv. Water Resour. 5 (1982) 98–104 June. [32] M.R. West, B.H. Kueper, K.S. Novakowski Semi-analytical solutions for solute transport in fractured porous media using a strip source of finite width, Adv. Water Resour. 27 (11), (2004) 1045–1059. [33] B.G. Deruyck, D.P. Bourdet, G. DaPrat, H.J. Ramey, Interpretation of interference tests in reservoirs with double porosity behavior—theory and field examples, SPE paper# 11205, Presented at 57th Annual Fall Technical Conference and Exhibition held in New Orleans, LA, September 26–29, 1982.

1362

I. Kocabas / Journal of the Franklin Institute 348 (2011) 1339–1362

[34] H. Cinco-Ley, V. Samniego, F. Kucuk, The pressure transient behavior for naturally fractured reservoirs with multiple block size, SPE Paper#14168, Presented at 60th Annual Fall Technical Conference and Exhibition held in Las Vegas, NV, September 22–25, 1985. [35] R. Pecher, J.F. Stanislav, C.V. Easwaran, Pressure transient testing of the heterogeneous reservoirs, J. Can. Pet. Technol. 34 (3) (1995) 19–26. [36] G. Zhao, L.G. Thompson, Modeling of complex geometry reservoirs, SPE Reservoir Eng. B (2000) 775–787. [37] M.M. Levitan, G.E. Crawford, General heterogeneous radial and linear models for well—test analysis, SPE J. 7 (2002) 131–1382 7 (2002) 131–138. [38] El-Khatib, Noaman A. F., Transient pressure behavior for a reservoir with continuous permeability distribution in the invaded zone, Proceedings 16th SPE Middle East Oil and Gas Show and Conference, MEOS 2009, vol. 2, pp. 609–624. [39] C.S. Chen., Analytical and approximate solutions to radial dispersion from an injection well to a geological unit with simultaneous fusion into adjacent strata, Water Resour. Res. 21 (8) (1985) 1069–1076. [40] D.H. Tang, D.W. Peaceman, New analytical and numerical solutions for the radial convection-dispersion problem, Soc. Pet. Eng. Reservoir Eng. 2 (3) (1987) 343–359. [41] A.F. Moench, Convergent radial dispersion: a Laplace transform solution for aquifer tracer testing, Water Resour. Res. 25 (3) (1989) 439–447. [42] M. Bai, Z. Shu, J. Coa, M. Zaman, J.C. Roegiers, A semi-analytical solution for a two-dimensional capacitance model in solute transport, JPSE 22 (1999) 275–292. [43] M. Coronado, J. Ram ´ırez-Sabag, Analytical model for tracer transport in reservoirs having a conductive geological fault, J. Pet. Sci. Eng. 62 (2008) 73–79. [44] H. Zhan, Z. Wen, G. Huang, D. Sun, Analytical solution of two-dimensional solute transport in an aquifer-aquitard system, J. Contaminant Hydrol. 107 (2009) 162–174.