Pressure transient behavior in a multilayer reservoir with

Prévia do material em texto

Journal of Petroleum Science and Engineering 208 (2022) 109376
Available online 31 August 2021
0920-4105/© 2021 Elsevier B.V. All rights reserved.
Pressure transient behavior in a multilayer reservoir with 
formation crossflow 
Jing Lu a,*, Md Motiur Rahman b, Erlong Yang a,**, Mohamed Tariq Alhamami b, Huiying Zhong a 
a Key Lab of Enhanced Oil Recovery, Ministry of Education, Northeast Petroleum University, Daqing, China 
b Department of Petroleum Engineering, Khalifa University of Science and Technology, Abu Dhabi, United Arab Emirates 
A R T I C L E I N F O 
Keywords: 
Analytical solution 
Multilayer reservoir 
Formation crossflow 
Pressure behavior 
A B S T R A C T 
This study proposes a new solution to the pressure transient behavior for a uniform-flux with a fully penetrating 
vertical well in a multilayer reservoir with formation crossflow. The general problem of pressure transient in 
multilayer reservoirs, in which any two adjacent layers are crossflowing in the formation, is solved analytically. 
The solution is presented for an infinite-acting system with no wellbore storage or skin factor effects. This 
analytical solution is obtained using Laplace transform, double Fourier transform, and Green’s functions method. 
This solution provides an accurate and fast tool to: evaluate a vertical well performance in a multilayer reservoir, 
and estimate the effects of formation properties on pressure behavior at locations both far away from and around 
the well. Based on the solution, the pressure of the vertical well, which produces with a constant rate in the 
multilayer reservoir, can be examined in detail. An expression of dimensionless crossflow coefficient is obtained 
by solving a symmetric tri-angular matrix equation. Based on the solution, it is verified that the dimensionless 
crossflow coefficient is not constant but behavies as a function of the time and distance away from the wellbore. 
When producing time is sufficiently long, crossflow will cease to exist and crossflow coefficients are equal to 
zero, which means there are no effects of adjacent layers on the pressure drop of each layer. These findings 
indicate that each layer produces independently. 
1. Introduction 
Petroleum reservoirs are seldom made up of a single homogeneous 
layer with constant properties but multiple layers with varying forma-
tion characteristics. Fluid flows into a layer or out of a layer and into an 
adjacent one, i.e., formation crossflow may occur when different layers 
are produced at different pressures and when layers exhibit different 
properties, and these layers affect each other through interlayer cross-
flow in reality. To this end, multilayer reservoirs with formation cross-
flow are often observed in oil field development. 
Whenever there is a pressure difference between two layers, cross-
flow will occur if there is communication between layers. A more 
permeable layer produces rapidly, which causes a more significant 
pressure drop. Thus, at the same distance from the wellbore, the for-
mation pressure becomes higher in the less permeable layer than in the 
more permeable layer, and the reservoir fluid starts flowing from the less 
to the more permeable layer. This crossflow phenomenon has many 
characteristic effects both on the response of the wellbore pressure and 
on the production rate from each layer. In a multilayer reservoir system, 
the pressure of any mid-layer bounded between an upper and a lower 
producing intervals depends on the amount of fluid flux gained or lost at 
the upper and lower interfaces. 
Interests in the behavior of multilayer reservoir systems with for-
mation crossflow have prompted many studies in the last four decades. 
Russell and Prats (1962) were the first to conduct research work on the 
multilayer system with crossflow. They stated that if crossflow occurs, 
the result is a shorter operating life and a higher primary ultimate re-
covery. Gao (1984) proposed that the interface between any two-layer 
could be represented as a single independent layer (semi-permeable 
wall) having specific pressure and vertical permeability values. Bourdet 
(1985) proposed that the two-layer reservoir with crossflow could be 
considered as the dual-porosity-dual-permeability naturally fractured 
reservoir. Ehlig-Economides and Joseph (1987) made a significant 
advance on the problem of the multilayer system with formation 
crossflow. They investigated both early and late time behaviors of the 
production rate of each layer, and they used the earlier semi-permeable 
* Corresponding author. 
** Corresponding author. 
E-mail addresses: nepulj@nepu.edu.cn (J. Lu), yangerlong@nepu.edu.cn (E. Yang). 
Contents lists available at ScienceDirect 
Journal of Petroleum Science and Engineering 
journal homepage: www.elsevier.com/locate/petrol 
https://doi.org/10.1016/j.petrol.2021.109376 
Received 20 January 2021; Received in revised form 10 August 2021; Accepted 12 August 2021 
mailto:nepulj@nepu.edu.cn
mailto:yangerlong@nepu.edu.cn
www.sciencedirect.com/science/journal/09204105
https://www.elsevier.com/locate/petrol
https://doi.org/10.1016/j.petrol.2021.109376
https://doi.org/10.1016/j.petrol.2021.109376
https://doi.org/10.1016/j.petrol.2021.109376
http://crossmark.crossref.org/dialog/?doi=10.1016/j.petrol.2021.109376&domain=pdf
Journal of Petroleum Science and Engineering 208 (2022) 109376
2
wall model presented by Gao (1984) to account for the crossflow 
behavior at the interface. Larsen (1988) illustrated the difficulty in the 
interpretation of the crossflow coefficient with a three-layer example. 
Later, Park and Horne (1989) extended Bourdet’s dual permeability 
model to any number of layers for different specifications and boundary 
conditions. Gomes and Ambastha (1993) summarized the models pub-
lished from 1960 to 1993 for a multilayer commingled reservoir; they 
pointed out the existing research mainly focused on analyzing and 
acquiring of physical parameters of stratification from the pressure and 
production dynamic data. Guo et al. (2002) discussed pressure responses 
from a formation with two communicating layers in which a fully 
penetrated high permeability layer is adjacent to a low-permeability 
layer. Al-Ajmi et al. (2003) used an analytically derived formula to es-
timate the storativity ratio of a layered reservoir with crossflow from 
pressure transient data. Dubost et al. (2015) presented a method to 
assess layer connectivity and determine dynamic reservoir properties for 
use as inputs in a numerical model for a multilayer reservoir, pressure 
transients recorded in different wells from a wireline formation tester in 
small scale drill stem test, or from conventional well tests, were then 
used to constrain the dynamic model properties through a grid-based 
inversion technique. Nooruddin and Rahman (2017) proposed an 
analytical procedure to estimate interlayer crossflow, and they assumed 
that a well completed only in the reservoir layer of interest, other 
adjacent layers communicated hydraulically with and contribute fluid to 
the layer of interest through reservoir cross flow only. Shi et al. (2020) 
et al. established a pressure transient behavior analysis model consid-
ering the vertical inhomogeneous closed boundary radii.Vasquez and 
Adrian (2021) presented the analytical derivation of the 
Palacio-Blasingame type curves to analyze production data of a 
two-layer reservoir model without crossflow or hydraulic communica-
tion between them. 
Multilayer reservoirs with crossflow are different from commingled 
stratified reservoirs, in a sense that their mathematical treatment is more 
complex. The proposed two-layer model of Bourdet (1985) suggests that 
for a dual-porosity, dual permeability system there are two overlapping 
media at any point in the formation; the matrix system and the fissure 
system. This is an incorrect description of the problem since there is only 
one medium at any point in a multilayer systemthere exists only one 
medium. Therefore it is believed that Bourdet’s model does not apply to 
multilayer reservoirs. 
Lu et al. (2019) challenged Bourdet’s dual permeability model and 
provided a new solution to the transient pressure equation of a vertical 
well in a two-layer reservoir with crossflow. They stated that the pres-
sure behavior in each layer of a two-layer system is controlled by the 
hydraulic diffusivity ratios and layer thickness ratios. They proposed a 
parameter to quantitatively describe the strength of the crossflow be-
tween adjacent layers. In this study, we propose an analytical procedure 
to calculate crossflow coefficients between the adjacent layers of a 
multilayer reservoir, which is not available in the literature. The analysis 
techniques offered by Lu et al. (2019) for their two-layer system are the 
precursors for this work. 
In this study, the general problem of n homogenous layers, in which 
any two adjacent flowing layers are crossflowing in the formation, is 
solved analytically. However, a different route is taken away from 
Bourdet’s dual-permeability model and Gao’s semi-permeable wall 
model, and an extension of the two-layer model that was developed by 
Lu et al. (2019) is made to model the n-layer system. The adopted 
methodology in the study serves to model the response of any arbitrary 
number of reservoir layers with formation crossflow. For simplicity, we 
only present the results for a three-layer system, but the formalism can 
easily be generalized to multilayer reservoir systems. Our proposed 
model gives an algorithm to calculate/predict the pressure of each layer 
and crossflow coefficients when a well is produced at a constant flow 
rate and completed in all layers in a multilayer reservoir, which have not 
been modelled previously. Futute development of this study is to pro-
pose a procedure to calculate/predict crossflow coefficients when the 
well is produced at a constant flowing bottom pressure and to extend the 
algorithm in this study to a well completed only in the layer of interest. 
Uncertainty of the measured results will also be presented in our future 
study. 
2. Model development 
The reservoir model considered in this study are the same as in (Gao, 
1984; Lu et al., 2019), which is shown schematically in Fig. 1. The 
reservoir consists of an arbitrary number of horizontal layers each of 
which are permeability isotropic, homogeneous, but different from each 
other. The top and the bottom of the reservoir are sealed by imperme-
able layers. A single well penetrates all the layers, which are initially in 
hydraulic pressure equilibrium. The reservoir fluid, produced at a con-
stant rate from the well head, is slightly compressible single phase fluid 
of constant viscosity. Wellbore storage effects and mechanical skin 
factors effects are ignored. Gravity force is negligible. 
2.1. Mathematical model for multilayer reservoir 
Fig. 1 shows the schematic of a multilayer reservoir model, a fully 
penetrating vertical well is producing in a multilayer reservoir with 
infinite extension in x and y directions. The assumptions for the reservoir 
are the same as stated in Lu et al. (2019). 
The reservoir initial pressure is 
P(r)|t=0 =Pini, (1a) 
and the reservoir pressure at infinity is equal to the reservoir initial 
pressure: 
P(t)|r→∞ =Pini. (1b) 
The center of the well is located at (0,0) in x-y coordinate system. 
The drainage domain is Ω = ( − ∞,∞)× (0, ht). Note that in Fig. 1, z 
direction is downward and 0 = z0 < z1 < z2 < z3 = ht. 
For each layer, there holds (Lu et al., 2019) 
∂2Pi
∂x2 +
∂2Pi
∂y2 +
∂2Pi
∂z2 =
(
μφici
ki
)
∂Pi
∂t
+
μqiδ(x)δ(y)
kihi
, (2) 
where Pi,ki, φi, ci, qi, hi are pressure, effective permeability, porosity, 
formation compressibility, flow rate and layer thickness of layer i, 
respectively. δ(x) and δ(y) are Dirac delta functions. The hydraulic 
diffusivity coefficient of layer i is defined as below (Craft and Hawkins, 
1991; Lee et al., 2003) 
ηi =
ki
μφici
, (i= 1, 2, ..., n), (3) 
and the hydraulic diffusivity coefficient ratio of layer i (Lu et al., 2019) 
σi =
η1
ηi
, (i= 1, 2, ..., n). (4) 
Obviously σ1 = 1. Then Equation (2) becomes 
(
1
ηi
)
∂Pi
∂t
−
(
∂2Pi
∂x2 +
∂2Pi
∂y2 +
∂2Pi
∂z2
)
=
− μqiδ(x)δ(y)
kihi
. (5) 
Note that 
∂P1
∂z
|z=0 =
∂Pn
∂z
|z=ht
= 0. (6) 
At the interface z = zi, (i = 1, 2, ...,n − 1), obviously there hold the 
following conditions (Lu et al., 2019): 
lim
(zi− 1 ,zi)∍
z→ziP(z)= lim
(zi ,zi+1)∍
z→ziP(z), lim
(zi− 1 ,zi)∍
z→zikihi
∂P(z)
∂z
= lim
(zi ,zi+1)∍
z→ziki+1hi+1
∂P(z)
∂z
(7) 
J. Lu et al. 
Journal of Petroleum Science and Engineering 208 (2022) 109376
3
If the flow rate in layer i is qi, then there holds (Gao, 1984) 
qi =
(
2πkihi
μB
)(
r
∂Pi
∂r
)
|r=Rw
, (i= 1, 2, ..., n). (8) 
Using the same definitions in Lu et al. (2019), the following 
dimensionless parameters are given: 
xD =
x
ht
, yD =
y
ht
, zD =
z
ht
, rD =
r
ht
, hD,i =
hi
ht
, tD =
k1t
μφ1c1h2
t
. (9) 
where ht =
∑n
i=1hi denote the height of the whole multilayer reservoir. 
Define the dimensionless pressure as in Lu et al. (2019): 
PD,i =
2πkihi(Pini − Pi)
μqi
, (i= 1, 2, ..., n). (10) 
By the same process in Lu et al. (2019), Equations (9) and (10) lead to 
the following equation 
σi
∂PD,i
∂tD
−
(
∂2PD,i
∂x2
D
+
∂2PD,i
∂y2
D
+
∂2PD,i
∂z2
D
)
= 2πδ(xD)δ(yD). (11) 
The dimensionless initial condition, top and bottom boundary con-
ditions are: 
PD,i|tD=0 = 0,
∂PD,1
∂zD
|zD=0 =
∂PD,n
∂zD
|zD=1 = 0. (12) 
Taking Laplace transform with respect to tD on both sides of Equation 
(11), we obtain (Haberman, 2004; Gradshteyn and Ryzhik, 2007) 
sσiPD,i −
⎛
⎝∂2PD,i
∂x2
D
+
∂2PD,i
∂y2
D
+
∂2PD,i
∂z2
D
⎞
⎠=
(
2π
s
)
δ(xD)δ(yD) (13) 
where s is the Laplace transform variable with respect to tD, PD,i is the 
image of PD,i after taking Laplace transform with respect to tD. 
Taking double Fourier transform with respect to xD and yD on both 
sides of Equation (13), then (Zwillinger, 1996; Tuma, 1997) 
−
∂2PD,i
∂z2
D
+
(
sσi +ω2
1 +ω2
2
)
PD,i =
2π
2πs
=
1
s
, (14) 
where ω1 and ω2 are double Fourier transform variables with respect to 
xD and yD,PD,i is the image of PD,i after taking double Fourier transform 
with respect to xD and yD (Zhang and Zhang, 2007). 
Equation (14) can be presented in a simple way: 
−
∂2PD,i
∂z2
D
+
PD,i
βi
=
1
s
, (15) 
where 
βi =
1
sσi + ω2
1 + ω2
2
, (i= 1, 2, ..., n). (16) 
The corresponding boundary conditions are changed as below: 
∂PD,1
∂zD
|zD=0 =
∂PD,n
∂zD
|zD=1 = 0. (17) 
Recalling Equation (7), Equation (15) must satisfy the following 
interface conditions: 
lim
(zDi− 1 ,zDi)∍
zD→zDiPD,i(zD)= lim
(zDi ,zDi+1)∍
zD→zDiPD,i+1(zD) (18a) 
lim
(zDi− 1 ,zDi)∍
zD→zDiβi
∂PD,i(zD)
∂zD
= lim
(zDi ,zDi+1)∍
zD→zDiβi+1
∂PD,i+1(zD)
∂zD
. (18b) 
Consequently, in each interval (zD,i− 1,zD,i), i = 1, 2, ..., n − 1, there 
hold 
−
∂2PD,i
∂z2
D
+
PD,i
βi
=
1
s
, (19) 
βi
∂PD,i
∂zD
|zD=zD,i− 1
= λ i− 1, βi
∂PD,i
∂zD
|zD=zD,i
= λ i, (i= 1, 2, ..., n − 1), (20) 
Because the top and bottom boundaries are impermeable, there holds 
λ0 = ​ λn = 0. (21) 
The undetermined coefficients λi (i= 1, 2, ..., n − 1) must satisfy 
lim
(zDi− 1 ,zDi)∍
zD→zDiβi
∂PD,i(zD)
∂zD
= lim
(zDi ,zDi+1)∍
zD→zDiβi+1
∂PD,i+1(zD)
∂zD
= λi. (22) 
Recalling Equation (18a), for λi (i= 1,2, ..., n − 1) there holds 
PD,i
(
zD,i
)
=PD,i+1
(
zD,i
)
(i= 1, 2, ..., n − 1). (23) 
And λi represents the magnitude of the crossflow on the interface 
Fig. 1. Schematic of a multilayer reservoir with a fully penetrating vertical well. 
J. Lu et al.Journal of Petroleum Science and Engineering 208 (2022) 109376
4
between the adjacent layers, λi is dimensionless crossflow coefficient in 
integral transform space. 
εi =
(
sσi + ω2
1 + ω2
2
)1/2
=
(
1
βi
)1/2
, (i= 1, 2, ..., n − 1). (24) 
By the same process in Lu et al. (2019), we obtain (Beck, 1992; 
Stakgold, 1998) 
PD,i
(
zD,i− 1
)
= −
[
λi− 1εi cosh
(
εihD,i
)
sinh
(
εihD,i
)
]
+
λiεi
sinh
(
εihD,i
)+
βi
s
, (i= 1, 2, ..., n),
(25) 
PD,i
(
zD,i
)
= −
[
λi− 1εi
sinh
(
εihD,i
)
]
+
λiεi cosh
(
εihD,i
)
sinh
(
εihD,i
) +
βi
s
, (i= 1, 2, ..., n). (26) 
Let i = 1 in Equation (26), and note that λ0 = 0, consequently, 
PD,1
(
zD,1
)
= −
[
λ0ε1
sinh
(
ε1hD,1
)
]
+
λ1ε1 cosh
(
ε1hD,1
)
sinh
(
ε1hD,1
) +
β1
s
=
λ1ε1 cosh
(
ε1hD,1
)
sinh
(
ε1hD,1
)
+
β1
s
.
(27) 
Recalling Equation (23), use Equations (25) and (26), for i = 1,2,...,
n − 1, there holds 
−
[
λi− 1εi
sinh
(
εihD,i
)
]
+
λiεi cosh
(
εihD,i
)
sinh
(
εihD,i
) +
βi
s
= −
[
λiεi+1 cosh
(
εi+1hD,i+1
)
sinh
(
εi+1hD,i+1
)
]
+
λi+1εi+1
sinh
(
εi+1hD,i+1
)+
βi+1
s
(28) 
Let i = n in Equation (25), and note that λn = 0, then 
PD,n
(
zD,n− 1
)
= −
[
λn− 1εn cosh
(
εnhD,n
)
sinh
(
εnhD,n
)
]
+
βn
s
. (29) 
Rearranging Equation (28), we obtain 
2.2. Dimensionless crossflow coefficient 
Equation (30) can be expressed as follows: 
⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎝
a11 a12 0 0 ... 0
a21 a22 a23 0 ... 0
0 a32 a33 a34 ... 0
0 0 ... ... ... 0
0 0 0 ... ... ...
0 0 0 ... an− 1,n− 2 an− 1,n− 1
⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎠
⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎝
λ1
λ2
λ3
...
...
λn− 1
⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎠
=
⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎝
g1
g2
g3
...
...
gn− 1
⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎠
(31) 
where 
ai,i =
εi cosh
(
εihD,i
)
sinh
(
εihD,i
) +
εi+1 cosh
(
εi+1hD,i+1
)
sinh
(
εi+1hD,i+1
) (i= 1, 2, 3, ..., n − 1) (32a) 
ai− 1,i = −
[
εi
sinh
(
εihD,i
)
]
= ai,i− 1(i= 2, 3, ..., n − 1) (32b) 
gi =
(
1
s
)
(βi+1 − βi), (i= 1, 2, ..., n − 1) (32c) 
The matrix on the left side of Equation (31) is a symmetric tri-angular 
matrix, and Equation (31) can be solved by using Cramer’s rule (Zwil-
linger, 1996; Tuma, 1997). 
For a four-layer reservoir, we have 
⎛
⎝
a11 a12 0
a21 a22 a23
0 a32 a33
⎞
⎠
⎛
⎝
λ1
λ2
λ3
⎞
⎠=
⎛
⎝
g1
g2
g3
⎞
⎠ (33) 
And a11, a22, a33 can be calculated by Equation (32a), a12, a21 can be 
calculated by Equation (32b), b1, b2, b3 can be calculated by Equation 
(32c). 
By using Cramer’s rule, we have 
λ1 =
g1a22a33 − g1a23a32 − g2a12a33 + g3a12a23
a11a22a33 − a11a23a32 − a12a21a33
​ (34a) 
λ2 =
g2a11a33 − g1a21a33 − g3a11a23
a11a22a33 − a11a23a32 − a12a21a33
​ (34b) 
λ3 =
g1a21a32 − g2a11a32 + g3a11a22 − g3a12a21
a11a22a33 − a11a23a32 − a12a21a33
​ (34c) 
For a three-layer reservoir, we have 
(
a11 a12
a21 a22
)(
λ1
λ2
)
=
(
g1
g2
)
(35) 
and 
λ1 = ​ (g1a22 − g2a12)/(a11a22 − a12a21) (36a) 
λ2 = ​ (g2a11 − g1a21)/(a11a22 − a12a21) (36b) 
The dimensionless crossflow coefficient in the real space can be 
calculated by (Zakian, 1970; Stehfest, 1979) 
ΛD,i
(
tD, rD, hD,i
)
= F− 1L− 1{λi(s,ω1,ω2)}, (i= 1, 2, ..., n − 1). (37) 
where L− 1 is inverse Laplace transform operator, F− 1 is inverse double 
Fourier transform operator. 
2.3. Dimensionless average pressure in each layer 
By a similar process in Lu et al. (2019), for i = 1, 2, ..., n, there holds 
(Beck, 1992; Stakgold, 1998) 
PD,i(s,ω1,ω1, zD)= −
λi− 1 cosh
[
εi
(
zD,i − zD
)]
εiβi sinh
(
εihD,i
) +
λi cosh
[
εi
(
zD − zD,i− 1
)]
εiβi sinh
(
εihD,i
)
+
βi
s
,
(
zD,i− 1 ≤ zD ≤ zD,i
)
(38) 
Note that λ0 = 0, zD,0 = 0, then 
PD,1(s,ω1,ω2, zD)=
λ 1ε1 cosh(ε1zD)
sinh
(
ε1hD,1
) +
β1
s
,
(
0< zD ≤ zD,1
)
(39) 
−
[
εi
sinh
(
εihD,i
)
]
λi− 1 +
[εi cosh
(
εihD,i
)
sinh
(
εihD,i
) +
εi+1 cosh
(
εi+1hD,i+1
)
sinh
(
εi+1hD,i+1
)
]
λi −
[
εi+1
sinh
(
εi+1hD,i+1
)
]
λi+1
=
1
s
(βi+1 − βi), (i = 1, 2, ..., n − 1)
(30) 
J. Lu et al. 
Journal of Petroleum Science and Engineering 208 (2022) 109376
5
Using λn = 0, zD,n = hD,t, consequently 
PD,n(s,ω1,ω2, zD)= −
λn− 1 cosh
[
εn
(
hD,t − zD
)]
εnβn sinh
(
εnhD,n
) +
βn
s
,
(
zD,n− 1 ≤ zD ≤ zD,n
)
(40) 
Equation (38) shows that PD,i(s,ω1,ω2, zD) is a function of zD, which 
indicates that the dimensionless pressure PD,i is a function of pressure 
measurement point, zD. And the average value of PD,i(s,ω, zD) can be 
calculated from Equation (38): 
Pave,D,i(s,ω1,ω2)=
(
1
hD,i
)∫ zD,i
zD,i− 1
PD,i(s,ω1,ω2, zD)dzD = −
(
λi− 1
hD,i
)
+
λi
hD,i
+
βi
s
, (i= 1, 2, ..., n)
(41) 
Then the dimensionless average pressure in Layer i can be calculated 
by: 
Pave,D,i(tD, rD)=F− 1L− 1
{
Pave,D,i(s,ω1,ω2)
}
. (42) 
where rD =
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
x2
D + y2
D
√
, 
Taking inverse double Fourier transform in Equation (41), we obtain 
(Haberman, 2004) 
Here PD is the image of PD after taking inverse double Fourier 
transform, and J0(x) is the Bessel function of the first kind of order zero. 
The Gauss-Laguerre quadrature formula with m points is utilized in 
Equation (43), (Zwillinger, 1996; Tuma, 1997), then 
Pave,D,i(s, rD) ≈
∑m
k=1
γ(m)
k ζ(m)
k Pave,D,i
(
s, ζ(m)
k
)
J0
(
rDζ(m)
k
)
. (44) 
Zakian’s method (Zakian, 1970, 1971) of numerical inverse Laplace 
transform is untilized in Equation (44), and the dimensionless pressure 
at tD and rD can be calculated by: 
Pave,D,i(tD, rD) ≈
∑m
k=1
γ(m)
k ζ(m)
k Pave,D,i
(
tD, ζ(m)
k
)
J0
(
rDζ(m)
k
)
. (45) 
3. Results and discussion 
The pressure transient model has been applied to a three-layer 
reservoir and the pressure behavior is examined. Fig. 2 shows the 
schematic of a three-layer reservoir model. The dimensionless formation 
crossflow coeffients are calculated. According to oil field practice, 
pressure is measured at the wellbore and we only calculate the dimen-
sionless pressure at sandface which is the physical interface between the 
formation and the wellbore. At sandface, the radial distance from 
wellbore is R = Rw. 
For a three-layer reservoir, λ1, λ2 can be obtained by Equations (36a) 
and (36b), respectively; and 
a11 =
ε1 cosh(ε1hD1)
sinh(ε1hD1)
+
ε2 cosh(ε2hD2)
sinh(ε2hD2)
(46a) 
a22 =
ε2 cosh(ε2hD2)
sinh(ε2hD2)
+
ε3 cosh(ε3hD3)
sinh(ε3hD3)
(46b) 
a12 = a21 = −
[
ε2
sinh(ε2hD2)
]
(46c) 
g1 =
(
1
s
)
(β2 − β1) (47a) 
g2 =
(
1
s
)
(β3 − β2) (47b) 
Recalling Equation (41), for layer 1, 
Pave,D,1(s,ω1,ω2)=
λ1
hD1
+
β1
s
(48a) 
For layer 2, 
Pave,D,2(s,ω1,ω2)= −
(
λ1
hD2
)
+
λ2
hD2
+
β2
s
(48b) 
Fig. 2. Schematic of a three-layer reservoir with crossflow. 
Pave,D,i
⎛
⎝s, xD, yD
⎞
⎠ =
⎛
⎝ 1
2π
⎞
⎠
∫∞
− ∞
∫∞
− ∞
exp[ − ι(ω1xD + ω2yD
⎞
⎠
⎤
⎦Pave,D,i
⎛
⎝s,ω
⎞
⎠dω1dω2
​ =
∫∞
0
ωPave,D,i
⎛
⎝s,ω
⎞
⎠J0
⎛
⎝rDω
⎞
⎠dω
​ = Pave,D,i
(
s, rD
)
,
(
i = 1, 2, ..., n; rD =
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
x2
D + y2
D
√
, ω =
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
ω2
1 + ω2
2
√
, ι2 = − 1
)
(43) 
J. Lu et al. 
Journal of Petroleum Science and Engineering 208 (2022) 109376
6
For layer 3, 
Pave,D,3(s,ω1,ω2)= −
λ2
hD3
+
β3
s
(48c) 
Note that when λ1 = λ2 = 0, then 
Pave,D,i(s,ω1,ω2)=
βi
s
=
(
1
s
)(
1
sσi + ω2
1 + ω2
2
)
, (i= 1, 2, 3) (49) 
Thus if there is no crossflow between the layers, by taking inverse 
double Fourier transform and inverse Laplace transform, we have 
PD,i(tD, rD)= −
1
2
Ei
(
−
σir2
D
4tD
)
, (i= 1, 2, 3) (50a) 
Substituting the dimensionless definitions, i.e., Equations (9) and 
(10) into Equation (50a), and rearranging the resultant expression, we 
have 
Pini − Pi(t, r)=
(
μqi
4πkihi)[
− Ei
(
r2
4ηit
)]
, (i= 1, 2, 3) (50b) 
Equation (50b) shows that if there is no crossflow between the layers, 
λ1 = λ2 = 0, no effects of adjacent layers on the pressure drop of each 
layer, i.e., each layer produces independently. 
Recalling Equation (4), hydraulic diffusivity coefficient ratio σ1 = 1,
and there hold 
σ2 =
η1
η2
=
(
k1
φ1c1
)/(
k2
φ2c2
)
, σ3 =
η1
η3
=
(
k1
φ1c1
)/(
k3
φ3c3
)
(51a) 
σ2
σ3
=
η3
η2
=
(
k3
φ3c3
)/(
k2
φ2c2
)
. (51b) 
The data of hydraulic diffusivity coefficient ratio and layer dimen-
sionless thickness in Table 1 are utilized to study the pressure behavior 
and crossflow for a three-layer reservoir in the following subsections, 
and we denote Layer 1 in Fig. 1 as top layer, Layer 2 as interlayer, and 
Layer 3 as bottom layer. 
3.1. Dimensionless average sandface pressure and formation crossflow 
coefficients 
Hydraulic diffusivity coefficient η is a measure of how fast the 
pressure drop propagates in porous media. At a given time, the bigger 
hydraulic diffusivity coefficient η, the bigger pressure drop ΔP in the 
reservoir is (Craft and Hawkins, 1991; Lee et al., 2003). Recalling 
Equation (4), hydraulic diffusivity coefficient ratio σi is inversely pro-
portional to hydraulic diffusivity coefficient ηi, and recalling Equation 
(10), dimensionless pressure PD,i is directly proportional to pressure 
drop ΔP = Pini − Pi So, the following conclusion can be reached: If there 
is no formation crossflow, at a given time, the bigger value of σi, the 
smaller value of PD,i is.In the following subsections and figures, PD,i 
means dimensionless average pressure of layer i, which is directly pro-
portional to pressure drop Pini − Pave,i. And Pave,i is average pressure of 
layer i. 
In this subsection, the data of Case 1 in Table 1 are utilized to study 
dimensionless average pressure in each layer at sandface and calculate 
dimensionless formation crossflow coefficients at the upper interface 
and the lower interface. The pressure behavior of a three-layer reservoir 
with a low permeability interlayer is demonstrated in Fig. 3, which 
shows dimensionless average pressure in each layer at sandface versus 
dimensionless time. Fig. 4 shows dimensionless crossflow coefficients 
Table 1 
Layer properties data for the three-layer reservoir in this study. 
σ1 σ2 σ3 hD1 hD2 hD3 
Case 1 1 32 0.2 1/3 1/3 1/3 
Case 2 1 2 0.2 1/3 1/3 1/3 
8 
32 
Case 3 1 100 0.2 0.45 0.1 0.45 
0.48 0.04 0.48 
0.49 0.02 0.49 
Fig. 3. Dimensionless average pressure in each layer at sandface. 
J. Lu et al. 
Journal of Petroleum Science and Engineering 208 (2022) 109376
7
versus dimensionless time. 
Fig. 4 also shows that when dimensionless time is small, the absolute 
values of dimensionless crossflow coefficients ΛD1,ΛD2 are small, which 
indicates the crossflow is weak. Because σ2 > σ1 > σ3, consequently 
when tD < 1, there holds PD3 > PD1 > PD2 in Fig. 3. Note 
thatσ2 = 32 > 1,η1 > η2, consequently, (Pini-Pave,1) in top layer is bigger 
than (Pini-Pave,2) in interlayer at a given time. ThusPave,1 − Pave,2 = (Pini −
Pave,2) − (Pini − Pave,1) < 0. Since the average pressure in top layer is 
Fig. 4. Dimensionless crossflow coefficients at sandface. 
Fig. 5. Dimensionless average pressure vs. dimensionless distance from wellbore. 
J. Lu et al. 
Journal of Petroleum Science and Engineering 208 (2022) 109376
8
smaller than that in interlayer, the direction of the crossflow is pointing 
from interlayer to top layer, which is the negative direction of z axial. So 
the dimensionless crossflow coefficient ΛD1 is negative as shown in 
Fig. 4. 
Note that σ3 = 0.2 < 1, and recalling Equation (51b), σ2/ σ3 = η3/
η2 > 1, η3 > η2, Consequently, (Pini-Pave,3) in bottom layer is bigger than 
(Pini-Pave,2) in intertlayer at a given time. So there holds. 
Pave,2 − Pave,3 = (Pini − Pave,3) − (Pini − Pave,2) > 0. The average 
pressure of interlayer is bigger than that of bottom layer, and ΛD2 is 
positive as shown in Fig. 4. 
Fig. 6. Dimensionless crossflow coefficients vs. dimensionless distance from wellbore. 
Fig. 7. Effects of hydraulic diffusivity coefficient ratio on dimensionless pressure. 
J. Lu et al. 
Journal of Petroleum Science and Engineering 208 (2022) 109376
9
Formation crossflow is caused by property differences between 
adjacent layers. When production time is small, the crossflow effect is 
weak compared with the influence of the property differences and each 
layer produces independently without any effect of crossflow. The ab-
solute values of crossflow coefficients increase with time until they 
reach the peak points in Fig. 4, then the effect of crossflow becomes 
strong. When time is big (tD > 105), the decrease of pressure difference 
between the adjacent layers will cause the decrease of formation 
crossflow. Eventually, no crossflow exists and ΛD = 0 in Fig. 4. Conse-
quently, very small pressure difference between any two pressure curves 
when time is big, and the segregations among the dimenionless pressure 
curves are small in Fig. 3. 
3.2. Effects of radial distance from wellbore 
In this subsection, the data of Case 1 in Table 1 are utilized to discuss 
how the pressure and crossflow coefficients change with the distance 
from wellbore. The effects of radial distance from wellbore on pressure 
and crossflow are demonstrated in Figs. 5 and 6. Fig. 5 shows the 
dimensionless average pressure of each layer versus dimensionless dis-
tance from wellbore when tD = 0.1. The radius of investigation, which is 
defined as the point in the formation beyond which the pressure drop is 
negligible, is a measure of how far a pressure transient has moved into a 
formation following any flow rate change in a well and physically rep-
resents the depth to which formation properties are being investigated at 
any time in a well test (Craft and Hawkins, 1991; Lee et al., 2003).Fig. 5 
also shows the smaller dimensionless pressure farther away from the 
wellbore at the given time.When tD = 0.1,RD > 1.0,PD = 0 pressure 
drop is zero, referring the definition of RD in Equation (9); so when tD =
0.1, the radius of investigation is equal to the total formation thickness 
ht. Since σ2 > σ1 > σ3, consequently at a given dimensionless distance, 
there holdsPD3 > PD1 > PD2. Fig. 6 shows the dimensionless crossflow 
coefficients versus dimensionless distance. Fig. 6 also shows that when 
tD = 0.1,RD < 1.0, dimensionless crossflow coefficients ΛD1 and ΛD2are 
almost constant, but when RD > 1.0,ΛD1 and ΛD2are quickly approaching 
to zero. Thus, beyond the radius of investigation, no crossflow exists. 
From the above analyses, the following conclusions can be reached: 
at a given time, pressure drop will decrease with the increasing distance 
from wellbore, and pressure drop is negligible at the location beyond the 
radius of investigation. Dimensionless formation crossflow coefficient is 
a function of the time and distance away from wellbore. When pro-
duction time is sufficiently long, crossflow will cease to exist; at a given 
time, crossflow will also cease to exist at the location which is far away 
from wellbore. 
3.3. Effects of hydraulic diffusivity coefficient ratio 
In this subsection, the data of Case 2 in Table 1 are utilized to 
demonstratethe effects of hydraulic diffusivity coefficient ratio on 
pressure behavior and formation corssflow coefficients. Fig. 7 presents 
the dimensionless average pressure of interlayer at sandface vs. 
dimensionless time. In Fig. 7, the hydraulic diffusivity coefficient ratio of 
interlayer has three different values of σ2 = 2, 8, 32. The other hy-
draulic diffusivity coefficient ratios areσ1 = 1, σ3 = 0.2.When other 
reserevoir parameters are constants, at a given time when the hydraulic 
diffusivity coefficient ratio of interlayer σ2 increases, the hydraulic 
diffusivity coefficient η2 decreases, and the average pressure drop (Pini- 
Pave,2) decreases; consequently, PD2 (the dimensionless average pressure 
of interlayer) decreases. So, Fig. 7 shows that at a given time, the value 
of PD2 corresponding to σ2 = 2 is the biggest, and the value of PD2 cor-
responding to σ2 = 32 is the smallest. 
When production time is sufficiently long, (tD > 106, ) Figs. 8 and 9 
show that diemsnionless crossflow coefficients are nearly equal to zero 
and the segregations among the dimenionless pressure curves are very 
small in Fig. 7. This is because crossflow does not exist. 
Fig. 8 shows ΛD1 (the dimensionless crossflow coefficient at sandface 
Fig. 8. Effects of hydraulic diffusivity coefficient ratio on ΛD1. 
J. Lu et al. 
Journal of Petroleum Science and Engineering 208 (2022) 109376
10
Fig. 9. Effects of hydraulic diffusivity coefficient ratio on ΛD2. 
Fig. 10. Effects of layer thickness on dimensionless pressure. 
J. Lu et al. 
Journal of Petroleum Science and Engineering 208 (2022) 109376
11
between top layer and interlayer) versus dimensionless time.The three 
curves in Fig. 8 are corresponding to σ2 = 2, 8, 32, respectively. 
Because σ2 = η1/η2 > 1, thus η1 > η2, so at a given time, the pressure 
drop (Pini-Pave,1) in top layer is bigger than the pressure drop (Pini-Pave,2) 
in intertlayer. So, there holds Pave,1 − Pave,2 = (Pini − Pave,2) − (Pini −
Pave,1) < 0. Consequently, the dimensionless crossflow coefficient ΛD1 is 
negative as shown in Fig. 8. A higher value of σ2 indicates a large 
property differences between top layer and interlayer. Fig. 8 also shows 
that a big value of σ2 will result in a big absolute value of dimensionless 
crossflow coefficient. 
Fig. 9 shows ΛD2 (the dimensionless crossflow coefficient at sandface 
between interlayer and bottom layer) vs. dimensionless time. The three 
curves in Fig. 9 are corresponding to σ2 = 2, 8, 32, respectively. Note 
that σ3 = 0.2, and recalling Equation (51b), there holds σ2/ σ3 = η3/
η2 >> 1, so at a given time, the pressure drop (Pini-Pave,3) in bottom layer 
is bigger than the pressure drop (Pini-Pave,2) in intertlayer. Thus Pave,2−
Pave,3 = (Pini − Pave,3) − (Pini − Pave,2) > 0 the direction of the crossflow is 
pointing from interlayer to bottom layer, which is the positive direction 
of z axis. Consequently, in Fig. 9, ΛD2is positive, and a large value of σ2 
will result in a larger value of dimensionless crossflow coefficient. 
3.4. Effects of layer thickness 
The data of Case 3 in Table 1 are utilized to demonstrate the effects of 
layer thickness on pressure behavior and formation corssflow coefficients. 
Fig. 10 shows the dimensionless average pressure of interlayer at sandface 
vs. dimensionless time. In Fig. 10, interlayer is very thin, the dimensionless 
thickness of interlayer has three different values: hD2 = 0.1, 0.04, 0.02,
and hD1 = hD3 = (1.0 − hD2)/2,σ1 = 1, σ2 = 100, σ3 = 0.2, thus inter-
layer has very small hydraulic diffusivity coefficient η2. 
Referring to Equation (8), in this study each layer is producing at 
different constant flow rates, the flow rate of interlayer q2 is a constant. 
When the thickness of interlayer decreases, in order to keep q2 as a 
constant, the pressure drop (Pini-Pave,2) in interlayer has to increase. Thus 
in Fig. 10, at a given time, with the smaller value of hD2, there is bigger 
value of PD2. When the producing time is small, the formation crossflow 
is weak, the segregations among the curves corresponding to different 
values of dimensionless thickness of interlayer are larger. When time 
increases, the effect of crossflow becomes strong, the influence of the 
thickness is gradually cancelled out by crossflow effect. Consequently, 
the segregations among the curves become small, which can be clearly 
observed in Fig. 10. 
Fig. 11 shows ΛD1(the dimensionless crossflow coefficient at sand-
face between top layer and interlayer) versus dimensionless time. The 
three curves in Fig. 11 are corresponding to hD2 = 0.1, 0.04, 0.02,
respectively. Note that when σ2 = η1/η2 = 100, andη1≫η2and when 
production time is small, the crossflow is weak, thus the pressure drop 
(Pini-Pave,1) in top layer is larger than the pressure drop (Pini-Pave,2) in 
interlayer. Consequently, Pave,1 − Pave,2 = (Pini − Pave,2) − (Pini −
Pave,1) < 0.So when time is small, (tD < 0.1) ΛD1 is negative in Fig. 11. 
Referring to Equation (7), fluid flux across interface must be 
continuous. In order to keep fluid flux across interface continuous, (Pini- 
Pave,2) in interlayer must increase due to very small thickness of inter-
layer (hD2 ≤ 0.1). When production time increases, the effect of cross-
flow becomes strong. Consequently, when tD > 10, Pave,1 − Pave,2 =
(Pini − Pave,2) − (Pini − Pave,1) > 0, thus ΛD1 is positive. When the thick-
ness of interlayer hD2 increases, referring to Equation (8), in order to 
keep q2as constant. 
(Pini-Pave,2) must decrease, which will lead a smaller value of Pave,1 −
Pave,2 = (Pini − Pave,2) − (Pini − Pave,1). Consequently, when 
10 < tD < 105, at a given time, the larger value of hD2and the smaller 
value of ΛD1can be clearly observed in Fig. 11. When production time is 
very large, (tD > 105), the effects of crossflow will disappear and η1 >>
η2, then (Pini − Pave,2) < (Pini − Pave,1), and ΛD1 becomes negative again as 
Fig. 11. Effects of interlayer thickness on ΛD1. 
J. Lu et al. 
Journal of Petroleum Science and Engineering 208 (2022) 109376
12
shown in Fig. 11. 
Fig. 12 shows ΛD2(the dimensionless crossflow coefficient at sand-
face between interlayer and bottom layer) vs. dimensionless time. The 
three curves in Fig. 12 are corresponding to three different values of 
dimensionless thickness of interlayer. Note that when σ3 = 0.2, and 
referring to Equation (51b), there holds σ2/σ3 = η3/η2 = 200 > > 1.So 
at a given time, the pressure drop (Pini-Pave,3) in bottom layer is bigger 
than the pressure drop (Pini-Pave,2) in intertlayer. Thus Pave,2− Pave,3 =
(Pini − Pave,3) − (Pini − Pave,2) > 0. Consequently, ΛD2 is always positive as 
shown in Fig. 12. 
When the thickness of interlayer hD2 increases, (Pini-Pave,2) must 
decrease to keep q2as a constant in Equations (7) and (8). Consequently, 
the values of Pave,2 − Pave,3 = (Pini − Pave,3) − (Pini − Pave,2)and ΛD2will 
increase. So when the crossflow is established, (tD > 1.0, ) at a given 
time, the larger values of hD2 and ΛD2 can be clearly observed in Fig. 12. 
4. Summary and conclusions 
Based on this study, the following conclusions are drawn: 
(1) A new solution of the pressure transient behavior is proposed for 
a vertical well in a multilayer reservoir with formation crossflow. 
The general problem ofpressure transient in multilayer reser-
voirs, in which any two adjacent layers are crossflowing in the 
formation, is solved analytically. 
(2) We obtained an expression of dimensionless crossflow coefficient 
by solving a symmetric tri-angular matrix equation. When pro-
ducing time is sufficient long, crossflow will cease to exist, 
crossflow coefficients are equal to zero. 
(3) The dimensionless crossflow coefficient is a function of the time 
and distance away from wellbore. At a given location, the 
crossflow coefficient becomes zero when production time is large. 
At a given time, the crossflow coefficient becomes zero when the 
distance is far away from wellbore. 
(4) The big property difference between two adjacent layers will lead 
to a big absolute value of dimensionless crossflow coefficient. 
Author contributions 
Jing Lu established the model, derived the equations and wrote the 
manuscript. Md. Motiur Rahman improved and revised the manuscript. 
Erlong Yang obtained financial support from the National Natural Sci-
ence Foundation of China and supervised this research project. 
Mohamed Tariq Alhamami and Huiying Zhong wrote MATLAB codes 
and conducted pressure data analysis. Jing Lu: the first author and 
corresponding author, established the model, derived the equations and 
wrote the manuscript. Md. Motiur Rahman: the second author, 
improved and revised the manuscript. Erlong Yang: the third author and 
corresponding author, obtained financial support from the National 
Natural Science Foundation of China and supervised this research 
project. Mohamed Tariq Alhamami: the fourth author, wrote MATLAB 
codes and conducted pressure data analysis. Huiying Zhong: the fifth 
author, wrote MATLAB codes and conducted pressure data analysis. 
Declaration of competing interest 
The authors declare that they have no known competing financial 
interests or personal relationships that could have appeared to influence 
the work reported in this paper. 
Acknowledgement 
This work was financially supported by the National Natural Science 
Foundation of China (51574085), the National Natural Science Foun-
dation of Heilongjiang province (E2016008). 
Fig. 12. Effects of interlayer thickness on ΛD2. 
J. Lu et al. 
Journal of Petroleum Science and Engineering 208 (2022) 109376
13
References 
Al-Ajmi, N.M., Kazemi, H., Ozkan, E., 2003. Estimation of storativity ratio in a layered 
reservoir with crossflow. In: SPE 84294, SPE Annual Technical Conference and 
Exhibition, Denver, Colorado, USA, 5-8 October. 
Beck, J.V., 1992. Heat Conduction Using Green’s Functions. Hemisphere Publishing 
Company, New York City, USA. 
Bourdet, D., 1985. Pressure behavior of layered reservoirs with crossflow. In: SPE 13628, 
SPE California Regional Meeting, Bakersfield, California, USA, 27–29 March. 
Craft, B.C., Hawkins, M., 1991. Applied Petroleum Reservoir Engineering. Prentice Hall 
Inc., Upper Saddle River, New Jersey, USA. 
Dubost, F.X., Rocha, M., Lavin, L.M., Fuentes, R., Ramirez, J.R., Guzman, D., 2015. 
Designing a development plan with a constrained multi-layered reservoir model 
using pressure transients and downhole fluid analysis. In: SPE 177256, SPE Latin 
American and Caribbean Petroleum Engineering Conference, Quito, Ecuador, 18–20 
November. 
Ehlig-Economides, C., Joseph, J., 1987. A new test for determination of individual layer 
properties in a multilayered reservoir. In: SPE Formation Evaluation, September, 
pp. 261–283. 
Gao, C.T., 1984. Single-phase fluid flow in a stratified porous medium with crossflow. 
SPE J. 97–106. February. 
Gomes, E., Ambastha, A., 1993. An analytical pressure transient model for multilayered, 
composite reservoirs with pseudo-steady state formation crossflow. In: SPE 26049, 
SPE Western Regional Meeting, Anchorage, Alaska, USA, 26-28 May. 
Guo, B., Stewart, G., Toro, M., 2002. Linearly supported radial flow - a flow regime in 
layered reservoirs. SPE Reservoir Eval. Eng. 103–110. April. 
Gradshteyn, I.S., Ryzhik, I.M., 2007. Table of Integrals, Series, and Products. Academic 
Press Inc., San Diego, California, USA. 
Haberman, R., 2004. Applied Partial Differential Equations with Fourier Series and 
Boundary Value Problems. Prentice Hall Inc., New York City, USA. 
Larsen, L., 1988. Similarities and differences in methods currently used to analyze 
pressure-transient data from layered reservoirs. In: SPE 18122, SPE Annual 
Technical Conference and Exhibition, Houston, Texas, USA, 2-5 October. 
Lee, J., Rollins, J.B., Spivey, J.P., 2003. Pressure Transient Testing. Society of Petroleum 
Engineers, Richardson, Texas, USA. 
Lu, J., Zhou, B., Rahman, M.M., He, X., 2019. New solution to the pressure transient 
equation in a two-layer reservoir with crossflow”. J. Comput. Appl. Math. 362, 
680–693. 
Nooruddin, H.A., Rahman, N.M.A., 2017. A new analytical procedure to estimate 
interlayer cross-flow rates in layered-reservoir systems using pressure-transient data. 
In: SPE 183689, SPE Middle East Oil & Gas Show and Conference, Manama, Bahrain, 
6–9 March. 
Park, H., Horne, R.N., 1989. Well test analysis of a multilayered reservoir with formation 
crossflow. In: SPE 19800, SPE Annual Technical Conference and Exhibition, San 
Antonio, Texas, USA, 8-11 October. 
Russell, D.G., Prats, M., 1962. The practical aspects of interlayer crossflow”. J. Petrol. 
Technol. 589–594. June. 
Shi, W., Cheng, S., Meng, L., Gao, M., Zhang, J., Shi, Z., Wang, W., Duan, L., 2020. 
Pressure transient behavior of layered commingled reservoir with vertical 
inhomogeneous closed boundary. J. Petrol. Sci. Eng. 189, 1–17. 
Stakgold, I., 1998. Green’s Functions and Boundary Value Problems. Academic Press 
Inc., San Diego, California, USA. 
Stehfest, H., 1979. Numerical Inversion of Laplace Transforms Algorithm. Association for 
Computing Machinery Journal, pp. 47–49. 
Tuma, J.J., 1997. Engineering Mathematics Handbook. McGraw-Hill Professional, New 
York City, USA. 
Vasquez, M.B., Adrian, P.M., 2021. Rate transient analysis in two-layered reservoir 
without crossflow. In: SPE 200906, SPE Trinidad and Tobago Section Energy 
Resources Conference Held Virtually 28 - 30 June. 
Zhang, S., Zhang, M., 2007. Solution of nonlinear dynamic differential equations based 
on numerical Laplace transform inversion. Appl. Math. Comput. 189, 79–86. 
Zakian, V., 1970. Optimization of numerical inversion of Laplace transforms. Electron. 
Lett. 6 (No. 21), 677–679. 
Zakian, V., 1971. Least-squares optimization of numerical inversion of Laplace 
transforms. Electron. Lett. 7, 71–72. 
Zwillinger, D., 1996. Standard Mathematical Tables and Formulae. CRC Press LLC, New 
York City, USA. 
J. Lu et al. 
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref1
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref1
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref1
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref2
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref2
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref3
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref3
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref4
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref4
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref5
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref5
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref5
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref5
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref5
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref7
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref7http://refhub.elsevier.com/S0920-4105(21)01024-X/sref7
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref8
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref8
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref9
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref9
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref9
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref10
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref10
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref11
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref11
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref12
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref12
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref13
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref13
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref13
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref14
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref14
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref15
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref15
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref15
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref16
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref16
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref16
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref16
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref17
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref17
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref17
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref18
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref18
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref19
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref19
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref19
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref20
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref20
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref21
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref21
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref22
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref22
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref23
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref23
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref23
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref24
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref24
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref26
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref26
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref27
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref27
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref28
http://refhub.elsevier.com/S0920-4105(21)01024-X/sref28
	Pressure transient behavior in a multilayer reservoir with formation crossflow
	1 Introduction
	2 Model development
	2.1 Mathematical model for multilayer reservoir
	2.2 Dimensionless crossflow coefficient
	2.3 Dimensionless average pressure in each layer
	3 Results and discussion
	3.1 Dimensionless average sandface pressure and formation crossflow coefficients
	3.2 Effects of radial distance from wellbore
	3.3 Effects of hydraulic diffusivity coefficient ratio
	3.4 Effects of layer thickness
	4 Summary and conclusions
	Author contributions
	Declaration of competing interest
	Acknowledgement
	References