Transp Porous Med (2010) 81:219240

DOI 10.1007/s11242-009-9397-y
Compressed Air Flow within Aquifer Reservoirs
of CAES Plants
R. Kushnir A. Ullmann A. Dayan
Received: 18 February 2009 / Accepted: 9 April 2009 / Published online: 20 May 2009
Springer Science+Business Media B.V. 2009
Abstract A model on the air ow within aquifer reservoirs of Compressed Air Energy
Storage (CAES) plants was developed. The design of such CAES plants requires knowledge
of the reservoir air pressure distribution during both the charging and discharging phases.
Also, it must assure air/water interface stability to prevent water suction during discharge. An
approximate analytical solution for the pressure variations within the anisotropic reservoir
porous space was developed, subject to the Darcy equation and for conditions of partially
penetrating wells. Sensitivity analyses were conducted to identify the dominant parameters
affecting the well pressure and the critical ow rate (water suction threshold). It is dem-
onstrated that water coning is a factor that could severely limit the discharge air ow rate.
A significant diminishment of that limitation and reduction of the pressure uctuation can
be achieved by enlargement of the air layer height and discharge period. Likewise, aquifers
with larger horizontal permeability impose less restrictive critical ows. A conclusion on the
preferred screen length could not be merely drawn from technological considerations, but
should also involve important economic aspects.
Keywords Compressed air energy storage (CAES) Porous reservoirs
Partially penetrating well Water coning Critical ow rates
List of Symbols
C Constant, C=e

=1.781072
CD Charging discharging time ratio
f Porosity
F(t

) Dimensionless air mass ow-rate at the well

g Gravitational acceleration
g

Dimensionless hydrostatic pressure,

w
gH/P
0
h Well screen length
R. Kushnir A. Ullmann (B) A. Dayan
Department of Fluid Mechanics and Heat Transfer, School of Mechanical Engineering,
Tel Aviv University, Tel Aviv 69978, Israel
e-mail: ullmann@eng.tau.ac.il
1 3
220 R. Kushnir et al.
h

Dimensionless well screen length, h/H

H Gas layer height
k Permeability
m
c
Air mass ow rate through the compressor
m

c
Dimensionless air mass ow-rate through the compressor,
m
c
ZRT
Hk
r
P
2
0
p Pressure
P
0
Initial air pressure in the reservoir
p

Dimensionless pressure, p/P

0
r Radial coordinate
r

Dimensionless radial coordinate,

_
k
z
k
r
r
H
r
w
Radius of the well
r

w
Dimensionless well radius,
_
k
z
k
r
r
w
H
R Specic gas constant
s Laplace transform parameter
t Time
t
p
Time period of the cycle
t
i
i =1, 2, 3, process duration times, see Fig. 2
t

Dimensionless time, t /t
p
t

Dimensionless elapsed time from the beginning of each stage

T Temperature
v Supercial velocity of the gas
y

Relative interface rise,

H
Hh
z Vertical coordinate
z

Dimensionless vertical coordinate, z/H

Z Gas compressibility factor
Greek Symbol
Pneumatic diffusivity, kP
0
/( f )
Modied pressure, dened in Eq. 21

Dimensionless modied pressure, /P

0
Interface location

Dimensionless interface location, /H

Viscosity of the gas

n
Lag angle
Density of the gas

w
Density of the water

z
Dimensionless cycle time period, t
p

z
/H
2
Subscript
cr Critical
f First cycle
r Radial
s Steady periodic cycle
z Vertical
1 3
Compressed Air Flow of CAES Plants 221
1 Introduction
Typical diurnal electric power consumption exhibits significant variations. It reaches its
peak during daylight and drops to its trough at nighttime. Storage of excess power capac-
ity during off-peak hours is desirable both economically and environmentally. In principle,
off-peak excess electrical energy is stored for subsequent use during hours of peak demand.
In this respect, incorporation of a Compressed Air Energy Storage (CAES) facility is highly
attractive. In a CAES plant, air is compressed into an underground reservoir through the con-
sumption of inexpensive excess electrical power. Later, during peak hours, the compressed
air is heated (red) and then driven to expansion in a gas turbine for electrical power genera-
tion. Three geologically different types of underground reservoirs are feasible: salt caverns,
hard rock caverns, and porous reservoirs (such as aquifers or depleted reservoirs). For the
design of a CAES plant, prediction of the underground reservoir pressure oscillations during
plant operation is required, and specifically the oscillations extreme values. The maximum
pressure constitutes a constraint that the compressor train must meet. The minimum pressure
essentially determines the turbine inlet pressure.
This study addresses the air storage in aquifer reservoirs. To date, there is no commercial
CAES plant that is linked to such reservoirs. However, experimental testing in a single-well
and two-well environments was conducted and conrmed that an aquifer reservoir is indeed a
suitable mediumfor compressed air energy storage (Allen et al. 1984, ANRStorage Company
1986). The experiments also showed that incorporation of a compressor after-cooler reduces
air temperature variations within the reservoir (the gas ow becomes essentially isothermal).
Based on Darcys law, an isothermal transient gas ow in porous media is described by
a single, nonlinear, partial differential equation (Muskat 1937). An approximate analytical
steady periodic solution of that equation was derived for one-dimensional radial ow by
Kushnir et al. (2008). The solution was developed for typical operating conditions, namely,
two periods of constant well mass ow-rate for the charging and discharging phases, and no
ow in between.
In general, the radial ow solution is applicable to studies of dry reservoirs, in which the
well can fully penetrate the air zone. In wet aquifer reservoirs, partially penetrating wells are
used to prevent undesired water suction during the discharge phase. Consequently, the air
ow near the well is two dimensional and axis symmetrical. Furthermore, instead of a xed
impervious bottom boundary condition, a uctuating airwater interface must be addressed.
During the discharge stage, a local pressure drop causes the interface to rise toward the well,
forming a cone shape. Muskat and Wyckoff (1935), who coined the term water coning,
indicated the existence of a critical cone height and a corresponding critical owrate. Beyond
that height, the interface becomes unstable and water would ow into the well. Accordingly,
it is necessary to avoid the condition of interface instability.
The air injection and withdrawal induces a two-phase (airwater) owprocess. Numerical
solutions of the two-phase ow equations for a partially penetrating single well subject to
a weekly cycle were presented by Wiles and McCann (1981). The work made part of an
extensive CAES research effort conducted at the Pacic Northwest Laboratory, U.S.A. (e.g.,
Wiles 1979). The reservoir pressure uctuations were found to be highly sensitive to the well
screen length. Also, it was shown that for certain operating conditions and reservoir char-
acteristics, water can enter the well. Meiri and Karadi (1982) developed a numerical model
for similar conditions (two phase, two dimensional, single well), but subject to a daily cycle.
They examined reservoir permeability effects and found that the airwater displacement is
strongly inuenced by it. The breadth of the transition zone (the zone where both air and
1 3
222 R. Kushnir et al.
water are present in the void space at various degrees of saturation) narrows with increasing
permeability.
If the transition zone is considerably narrower than the air ow domain, it is justied to
assume that a sharp airwater interface would be representative, where a single-phase air ow
exists in the upper domain. Based on that assumption, Braester and Bear (1984) developed
a two-dimensional numerical model for a partially penetrating single well subject to a daily
periodic cycling. The model accounts for variations of the airwater interface location. In the
range of the parameters studied, the airwater interface uctuations were found to be small.
Calculations indicated that the well pressure uctuations were pronounced in cases of short
well penetrations (smaller than 20% of the gas layer height).
Clearly, application of each of the aforementioned numerical solutions is restricted to the
specic range of parameters for which it was developed. As such, those investigations cannot
reveal the general sensitivity of both the well pressure and water coning to the reservoir char-
acteristics and operating conditions. Furthermore, no quantitative data on the critical owrate
have been reported. Additionally, excluding the study of Meiri and Karadi (1982), all other
investigations addressed only the rst cycle solution, to minimize computation cost. There-
fore, there exists a need for explicit formulae that cover a wide range of parameters. In order
to address that need, this study was undertaken with the purpose of providing an analytical
tool to calculate the reservoir pressure distribution, and to determine the interface stability
during a plant operation. An approximate analytical solution for the reservoir pressure distri-
bution was developed for typical operating conditions, namely, two periods of constant well
mass ow-rate for the charging and discharging phases and no ow in between. Based on
the solution, the critical ow rate and interface height were determined in accordance with
the concepts of Muskat and Wyckoff (1935).
2 Formulation of the Problem
Consider a partially penetrating well located in an aquifer reservoir as shown in Fig. 1. Ini-
tially, the porous reservoir contains a stationary air layer of height H, bounded by a water
aquifer below and an impervious caprock above. During a CAES plant operation, air ows
into and out of the reservoir cyclically. Owing to that ow, the airwater interface moves
downward in the charging stage and upward during the withdrawal stage. Two-phase ow
phenomena are expected to periodically affect a relatively small volume of the reservoir,
located near the airwater interface. Such effects are, therefore, likely to be of limited signif-
icance. Hence, it is assumed that the airwater interface is sharp, meaning that the air-bubble
domain does not contain movable water (i.e., only water at irreducible saturation). For such
conditions, the continuity and momentum equations (neglecting gravitational effects) for the
air region, subject to the generalized gas state equation, are as follows:
( f )
t
+
1
r
(rv
r
)
r
+
(v
z
)
z
= 0 (1)
v
r
=
k
r

p
r
, v
z
=
k
z

p
z
(2)
=
p
ZRT
(3)
where v
r
and v
z
are the gas radial and vertical supercial velocities, f is the porosity, and
k
r
and k
z
denote the medium permeability in the radial (horizontal) and vertical directions,
respectively. The remaining symbols are consistent with common notations. The model is
1 3
Compressed Air Flow of CAES Plants 223
H
h
r
w
r
z
Undisturbed interface
p(r,z,t)
Compressed Air
Discharge stage interface
Charging stage interface
Water aquifer
(r,t)
Fig. 1 Schematic of a partially penetrating well in a CAES aquifer reservoir
based on the assumption that the reservoir can adequately be represented as a homogeneous
porous space with both constant effective porosity and permeability.
The air ow is essentially isothermal owing to both air cooling following the compression
stage and the immense thermal inertia of the porous medium(as compared to that of the gas).
Indeed, as aforementioned, the eld test data (Allen et al. 1984) revealed that air temperature
variations in the reservoir are minor when a compressor after cooler is present. Substitution
of Eqs. 2 and 3, for a constant temperature, into Eq. 1 yields the following nonlinear partial
differential equation:
p
t
=
k
r
2f
_

2
p
2
r
2
+
1
r
p
2
r
_
+
k
z
2f

2
p
2
z
2
. (4)
The initial reservoir gas pressure, P
0
, is uniform and equal to the local hydrostatic head
(the reservoir gas static head is of negligible inuence). The cyclical air injection and with-
drawal produce pressure uctuations in the stored air. Given that the pressure uctuations
are smaller than P
0
and subject to the isothermal-ow assumption, it is reasonable to con-
sider the uid viscosity and compressibility factor as constants and, respectively, equal to
= (T, P
0
), Z = Z(T, P
0
).
The pertinent initial and boundary conditions for the indicated coordinate system are
t = 0, p = P
0
(5)
z = 0, 0 < r < ,
p
z
= 0 (6)
0 z h, r 0,
h
_
0
2rv
r
dz = m
c
F(t ) (7)
h < z H, r = 0,
p
r
= 0 (8)
0 z H, r , p P
0
(9)
where h is the well screen length (penetration length), and the product m
c
F(t ) represents
the mass ow rate of the gas at the well during plant operation. m
c
is the air mass-ow rate
through the compressor, and F(t ) is a dimensionless periodic function with a cycle time
t
p
. Figure2 shows the variations of F(t ) of a CAES plant operating with constant air mass
ow-rates through the compressor and the turbine. The indicated time intervals are t
1
for the
charging time, t
2
t
1
for the storage time, and t
3
t
2
for the power generation time, and CD
represents the discharging to charging mass ow ratio (also equal to the ratio of the charging
1 3
224 R. Kushnir et al.
Fig. 2 The dimensionless air
mass ow-rate at the well during
a CAES plant cycle
Discharge
-CD
1
Storage Storage
Charging
t
p
1
t t
2
t
3
t
F(t)
time to the discharging time). The mass ow boundary condition is applied at r 0 rather
than at the well radius r = r
w
. This is justied since for typical operating conditions and
reservoir characteristics, seconds within the operation initiation of each stage, asymptotic
conditions are reached with a reservoir pressure distribution that is independent of the wells
radius (Kushnir et al. 2008). Consequently, the well is represented as a line source instead
of a cylindrical source. The boundary condition (9) addresses a large reservoir whose effec-
tive radius exceeds the uctuating gas penetration radii (the effective region for gas storage
around the well). However, as will be discussed later, the solution can be extended to account
for nite reservoir boundaries and cases of multiple wells.
In order to complete the model, the boundary condition at the airwater interface must
be formulated. As indicated in Fig. 1, the interface coordinates are specied by (r, t ). The
kinematic boundary condition requires that
D( z)
Dt
=

t
+
v
r
f

r

v
z
f
= 0 on z = . (10)
Neglecting the capillary pressure, at equilibrium condition, the pressures on both sides of the
interface are identical. As seen later, excluding short transition periods at the beginning of
each stage, the pressure variation rate near the well is moderate. Consequently, the interface
movement rate is slow. It is, therefore, reasonable to assume that the water layer beneath the
air ow zone is quasi-static. Hence, the position of the interface is described by
p = P
0

w
g (H z) on z = . (11)
The left-hand side of Eq. 11 represents the air pressure at the interface, and the right-hand
side (RHS) stands for the water pressure at the interface under hydrostatic conditions, with

w
being the water density. In order to obtain the kinematic boundary condition in terms of
the air pressure, the derivatives of Eq. 10 (with respect to ) are derived from Eq. 11.
Choosing P
0
and t
p
as the pressure and time scales, respectively, and H as the length scale
in the vertical direction, the dimensionless form of Eqs. 411 becomes
p

=

z
2
_

2
p
2
r
2
+
1
r

p
2
r

+

2
p
2
z
2
_
(12)
t

= 0, p

= 1 (13)
z

= 0, 0 < r

< ,
p

= 0 (14)
0 z

, r

0, r

_
0
p
2
r

dz

= m

c
F(t

) (15)
h

< z

1, r

= 0,
p

= 0 (16)
1 3
Compressed Air Flow of CAES Plants 225
0 z

1, r

, p

1 (17)
z

,
p

z
_
_
p

_
2
+
_
p

_
2
_
+ g

z
p

= 0, (18)
where the location of the interface is
z

, p

= 1 g

_
1 z

_
(19)
and
p

p
P
0
, t

t
t
p
, r

_
k
z
k
r
r
H
, z

z
H
, h

h
H
,

z

t
p
k
z
P
0
H
2
f

t
p

z
H
2
, m

c

m
c
ZRT
Hk
r
P
2
0
, g

w
gH
P
0
(20)
Equations 1219 could be reduced to a simpler formthrough the introduction of a modied
pressure , according to
p
2
= P
2
0
+P
0
, i.e., p
2
= 1 +

. (21)
The substitution of Eq. 21 into Eqs. 1219 yields

= p

z
_

r
2
+
1
r

+

2

z
2
_
(22)
t

= 0,

= 0 (23)
z

= 0, 0 < r

< ,

= 0 (24)
0 z

, r

0, h

= m

c
F(t

) (25)
h

< z

1, r

= 0,

= 0 (26)
0 z

1, r

0 (27)
z

z
2p

_
_

_
2
+
_

_
2
_
+ g

= 0 (28)
and the location of the interface is
z

,
_
1 +

_
1/2
= 1 g

_
1 z

_
. (29)
The integrated well ux (Eq. 15) is represented by a constant average ux (Eq. 25). This
representation entails a negligible error for homogenous reservoirs (Muskat 1937; Ruud and
Kabala 1997), and is clearly of advantageous simplicity.
In terms of the modied pressure, the nonlinearity of the differential equation is repre-
sented solely by the coefcient p

z
. For efcient operation, namely, one that is characterized
by small pressure uctuations relative to P
0
, Eq. 22 can be linearized by substituting p

= 1.
Indeed, for typical operating conditions and reservoir characteristics, the nonlinearity effects
on Eq. 22 are small (subsequently discussed). Hence, for practical conditions, the solution
of the linearized equation is sufciently accurate to represent the full nonlinear equation
solution.
Nonlinearity exists in the airwater interface boundary condition (Eq. 28), which further-
more refers to the unknown surface z

. Nonetheless, for small pressure uctuations,

1 3
226 R. Kushnir et al.
the airwater interface level would exhibit corresponding small variations. As a rst approx-
imation, one may neglect the inuence of the interface shape on the pressure distribution by
assuming a xed airwater interface (represented as an impervious boundary to the air ow
domain). With that simplication, Eq. 28 is replaced by
z

= 1, 0 < r

< ,

= 0. (30)
In fact, Eq. 30 can also be derived from Eq. 28, following perturbation solution method tech-
niques (based on the rst expansion term, while neglecting the interface movement).
From the pressure distribution solution, the shape of the interface z

(r

, t

) during
discharge can be determined from Eq. 29. However, under certain conditions of ow in the
air zone, the interface becomes unstable and water could be sucked into the well. In order
to avoid such occurrence, the conditions that ascertain the airwater interface stability are
subsequently explored.
3 Pressure Distribution
Evidently, from the imposed periodic boundary condition, the solution for the modied
pressure can be constructed from two separate parts (transient and steady periodic). The
transient part, which is affected by the initial condition, vanishes practically within a few
cycles. Therefore, the conditions within the reservoir are determined by the steady periodic
solution.
3.1 Fourier Series Method
A straightforward method is used to obtain the steady periodic solution of Eqs. 2227 and
30, through a Fourier series representation of F(t

). First, the equations are solved for the

mass ow rate expressed as a harmonic function of time, namely,
F(t

) = sin(2nt

+
n
). (31)
It is convenient to introduce the steady periodic modied pressure distribution for an
unbounded reservoir caused by a periodic source element distribution at z

= along the
well (Carslaw and Jaeger 1959, p. 263)
d

s
=
m

c
d
2h

_
n

z
(r
2
+(z

)
2
)
_
r
2
+(z

)
2
sin
_
2nt

+
n

_
n

z
_
r
2
+(z

)
2
_
_
, (32)
where m

c
d/2h

is the element strength amplitude with d being the element dimensionless

length. In order to satisfy the boundary conditions (24) and (30), a corresponding set of source
elements along the z

axis at points (0, 2m ), m = 0, 1, 2, . . . must be accounted for.

The solution can then be obtained by integrating the derived expression for all s from 0 to
h

. Consequently, the steady periodic solution of Eqs. 2227 and 30 (with p

= 1), subject
to F(t

) as expressed in (31), is

s
=
m

c
2h

m=
z

2m+h

_
z

2mh

_
n

z
(r
2
+
2
)
_
r
2
+
2
sin
_
2nt

+
n

_
n

z
_
r
2
+
2
_
_
d. (33)
1 3
Compressed Air Flow of CAES Plants 227
In the derivation, for simplicity, the variable transformation z

2m = was incor-
porated. Generally, the solution is a fast converging series. This is especially valid for smaller

z
where the contributions of remote elements become negligible.
The solution for a harmonic mass ow rate is a building block for the construction of any
periodic mass ow function. This is accomplished with a sine Fourier series of F(t

), which
effectively is a mere superposition of harmonic forcing functions for which the solution (33)
applies. The solution of the actual mass ow rate (see Fig. 2), obtained through the Fourier
series application, is therefore

s
=
m

c
h

n=1
C
n
n

m=
z

2m+h

_
z

2mh

_
n

z
(r
2
+
2
)
_
r
2
+
2
sin
_
2nt

+
n

_
n

z
(r
2
+
2
)
_
d, (34)
where
C
n
=
_
A
2
n
+ B
2
n
,
n
= arg(B
n
+ A
n
i ),
A
n
= cos
_
nt

1
_
sin
_
nt

1
_
CDcos
_
n(t

3
+t

2
)
_
sin
_
n(t

3
t

2
)
_
, and (35)
B
n
= sin
2
_
nt

1
_
CDsin
_
n(t

3
+t

2
)
_
sin
_
n(t

3
t

2
)
_
Equation34 is an exact steady periodic solution to the linearized problem. However, in the
time vicinity of a owstep change (see Fig. 2), the solution series (for n) converges extremely
slowly for small values of r

. Hence, the Fourier series solution does not reveal the maximum
and minimum well pressures, occurring at t

1
and t

3
, respectively. Nevertheless, the solution
is useful for all other times, and it also serves as a validation tool of the Laplace transform
solution, subsequently presented.
Although the plant operating conditions are determined by the steady periodic solution,
published numerical studies usually focused only on the rst cycle solution. Hence, it is of
interest to obtain the rst cycle solution and compare the limiting conditions of both solu-
tions. The solution for the charging stage of the rst cycle (when F(t

) = 1) can be obtained
in a similar way as applied to the steady periodic solution method, where the steady periodic
point source is replaced by a continuous source. The modied pressure for 0 t

1
is
therefore

f
=
m

c
2h

m=
z

2m+h

_
z

2mh

1
_
r
2
+
2
erfc
_
_
_
_
r
2
+
2
_
4
z
t

_
_
d 0 t

1
. (36)
The values of

f
in the subsequent times can be readily found by superposition

f 2
(r

, z

, t

) =

f
(r

, z

, t

f
(r

, z

, t

1
) t

1
< t

2
(37a)

f 3
(r

, z

, t

) =

f 2
(r

, z

, t

) CD

f
(r

, z

, t

2
) t

2
< t

3
(37b)

f 4
(r

, z

, t

) =

f 3
(r

, z

, t

) +CD

f
(r

, z

, t

3
) t

3
< t

1. (37c)
Equations 36and37represent anexact rst cycle solutiontothe linearizedproblem. The series
solution converges rapidly and is useful for calculating the rst cycle pressure
variation.
1 3
228 R. Kushnir et al.
3.2 Laplace Transform Method
In order to circumvent the slow convergence of the Fourier series, an alternative solution
based on the Laplace transform is incorporated. The method provides a steady periodic
pressure solution that can be easily evaluated for all values of r

. Applying, rst, the Laplace

transform to Eqs. 2227 and 30 (with p

= 1), and subsequently, the nite Fourier cosine

transform, yields the solution

(r

, z

, s) = m

c

F (s) K
0
__
s

z
r

_
+
2 m

c

F (s)
h

m=1
K
0
__
(m)
2
+
s

z
r

_
sin(mh

) cos(mz

)
m
, (38)
where

(r

, z

, s) and

F(s)denotes the Laplace transform of

(r

, z

, t

) and F(t

),
respectively, and K
0
is the modied Bessel function of the second kind of order zero. The
rst term on the RHS of (38), represents the solution of a fully penetrating well radial ow
(h

= 1). The second term describes the vertical pressure variations around a partially
penetrating well.
The rst cycle compression solution can be found by the substitution

F (s) = 1/s (i.e.,
F(t

) = 1) in Eq. 38. The expression obtained is then integrated, according to the inversion
theorem. Consequently, the modied pressure for 0 t

1
is

f
=
m

c
2
Ei
_
r
2
4
z
t

_
+
2 m

c
h

m=1
_
_
K
0
(mr

)
e
m
2

z
t

_
0
J
0
(r

)e

z
t

2
+(m)
2
d
_
_
sin(mh

)cos(mz

)
m
, 0 t

1
(39)
where Ei, is the exponential integral, and J
0
is the Bessel function of the rst kind of order
zero. The values of

f
in subsequent times can be found by Eqs. 37ac. Equation39, as
it stands, is merely an equivalent form of Eq. 36, and has no advantage on the latter. Actu-
ally, it carries a disadvantage, since as r

0, the series containing K

0
(mr

) requires an
extremely large number of terms before converging (and is singular at r

= 0). As shown sub-

sequently, that problem can be avoided by replacing the Bessel series with its corresponding
power series representation in terms of r

.
The superiority of the Laplace transform method is clearly apparent when considering the
steady periodic solution. The procedure used to obtain the steady periodic solution (follow-
ing the inversion theorem incorporation) is described by Kushnir et al. (2008). Applying that
procedure on Eq. 38 yields

s
=

f
+ m

c
I
0
+
2 m

c
h

m=1
I
m
e
m
2

z
t
sin(mh

) cos(mz

)
m
0 t

1, (40)
where
I
m
=

_
0
1 e
(
2
+m
2

z
)t

1
+CD
_
e
(
2
+m
2

z
)t

3
e
(
2
+m
2

z
)t

2
_
(1 e
(
2
+m
2

z
)
)(
2
+m
2

z
)
J
0
_
r

1/2
z
_
e

2
t

d
m = 0, 1, 2, ... (41)
1 3
Compressed Air Flow of CAES Plants 229
Equation40 is an exact steady periodic solution of the linearized partially penetrating well
representation. For all relevant periods, the solution is a fast converging series, and the
integrals I
m
can be easily calculated. Consequently, the steady periodic modied pressure
expressed by Eq. 40 is far more suitable for pressure computations at small values of r

than
its equivalent Fourier series form (Eq. 34). Furthermore, inspection of the Laplace transform
solutions reveals that for
2

z
t

>> 1 (where t

is the dimensionless elapsed time fromthe

beginning of each stage), transient effects in the vertical direction are practically negligible.
For such times, the solutions could be significantly simplied by neglecting those transient
effects. Accordingly, the steady periodic solution is reduced to

s
=

s,r
+
2 m

c
F (t

)
h

m=1
K
0
_
mr

_
sin(mh

) cos(mz

)
m
, (42)
where

s,r
is the steady periodic solution of a fully penetrating well given by
2

s,r
m

c
= Ei
_
r
2
4
z
t

_
+2I
0
0 < t

1
(43a)
2

s,r
m

c
= Ei
_
r
2
4
z
t

_
+ Ei
_
r
2
4
z
(t

1
)
_
+2I
0
t

1
< t

2
(43b)
2

s,r
m

c
= Ei
_
r
2
4
z
t

_
+ Ei
_
r
2
4
z
(t

1
)
_
+CDEi
_
r
2
4
z
(t

2
)
_
+2I
0
t

2
< t

3
(43c)
2

s,r
m

c
= Ei
_
r
2
4
z
t

_
+ Ei
_
r
2
4
z
(t

1
)
_
+CDEi
_
r
2
4
z
(t

2
)
_
CDEi
_
r
2
4
z
(t

3
)
_
+2I
0
t

3
< t

1. (43d)
At the vicinity of the well (small values of r

), asymptotic representations of the expo-

nential integral and of I
0
can be incorporated into the pressure expressions, which yield
2

s,r
m

c
= ln
4
z
t

Cr
2
+ I
0a
0 < t

1
(44a)
2

s,r
m

c
= ln
t

1
+ I
0a
t

1
< t

2
(44b)
2

s,r
m

c
= ln
t

1
CDln
4
z
(t

2
)
Cr
2
+ I
0a
t

2
< t

3
(44c)
2

s,r
m

c
= ln
t

1
+CDln
t

3
t

2
+ I
0a
t

3
< t

1, (44d)
where I
0a
is the asymptotic representation of 2I
0
I
0a
= ln
(1 +t

1
)
(1 +t

)
+CDln
(1 +t

2
)
(1 +t

3
)
. (45)
is the Gamma function and C = e

, with being the Eulers constant (C=1.781072).

Equation 42 is valid when
2

z
t

>> 1 for 0 < t

1
, and
2

z
(t

1
) >> 1 for
t

1
< t

2
, etc. In effect, these conditions state that the elapsed time from the begin-
ning of each stage should be substantially larger than the characteristic time H
2
/(
z

2
).
1 3
230 R. Kushnir et al.
For such times, pressure transient effects in the vertical direction are negligible, and the
vertical pressure changes are of a steady periodic nature. Notice that the rst cycle solution
can be calculated by Eq. 42, simply from the omission of the 2I
0
(or I
0a
) term. Thus, the
Laplace transform method also exposes the difference between the rst cycle and the steady
periodic cycle solutions.
As seen from Eq. 42, the modied pressure is practically independent of z

for r

s closer
or larger than 2. Effects of the partial penetration of the well are therefore conned to a
dimensionless radial distance smaller than 2. Beyond that distance, the ow is essentially
radial, and the modied pressure can be calculated by

s,r
. Up to a distance of 2, the vertical
ow is prominent, as obtained from the summation term of (42). For smaller r

s, more
terms are needed to calculate the Bessel series. In particular, as r

approaches zero, the series

requires an extremely large number of terms to calculate. For small values of r

, it is more
advantageous to address the two equivalent solutions, Eqs. 36 and 39. On comparing those
solutions, subject to
z
t

, and by expressing (r
2
+
2
)
1/2
in terms of a power series
of r

, the following representation emerges

2
h

m=1
K
0
(mr

)
sin(mh

) cos(mz

)
m
= ln
r

4

1
2h

_
ln
_

2
((z

+h

) /2) sin ( (z

+h

) /2)

2
((z

) /2) sin ( (z

) /2)
_
ln
_
z

+h

+h

+
_
r
2
+(z

+h

)
2
z

+
_
r
2
+(z

)
2
_
+
r
2
16
_

_
2, 1
z

2

h

2
_

_
2, 1
z

2
+
h

2
_
+
_
2,
z

h
2

_
2,
z

+h
2

_
+
_
2
z

+h

_
2

_
2
z

_
2
_

3r
4
512
_

_
4, 1
z

2

h

2
_

_
4, 1
z

2
+
h

2
_
+
_
4,
z

h
2

_
4,
z

+h
2

_
+
_
2
z

+h

_
4

_
2
z

_
4
_
+
_
, (46)
where is the Hurwitz Zeta function. Excluding extreme cases, the value of r

at the well
is small enough so that only the logarithmic terms, on the RHS of (46), must be retained to
calculate the well pressure. Furthermore, the term ln(r

/4) is the only singular function at

r

= 0. As expected, it cancels out when combined with the singular term of

s,r
(see Eqs.
44a, 44c). Consequently, a nite value of the pressure beneath the well (r

= 0, z

> h

)
can be calculated. The pressure variations beneath the well are of particular significance, by
providing the basis for the calculation of the highest interface rise, which occurs at r

= 0.
As aforementioned, the maximumand minimumreservoir pressures are of practical inter-
est. These pressures occur at the top of the well (r

= r

w
, z

= 0) at times t

1
and t

3
,
1 3
Compressed Air Flow of CAES Plants 231
(a) (b)
Fig. 3 First period air pressure variations at the well (r

= r

w
, z

= 0): comparison between Eq. 42 and the

exact solution for m

c
= 0.005, r

w
= 0.014,
z
= 100, h

= 0.5, t

1
= 7/24, t

2
= 14/24, and t

3
= 18/24.
a Charging and subsequent storage and b discharge and subsequent storage
respectively. Upon substituting the indicated location and times in Eq. 42 and 46 one obtains,
2

s max
m

c
= ln

z
t

1
4C
+ I
0a
+
1
h

_
_
ln
_
_
r
2
w
+h
2
+h

2
(h

/2)
_
_
r
2
w
+h
2
h

2
(h

/2)
+ O(r
2
w
)
_
_
2

s min
m

c
= ln
t

3
t

3
t

1
CDln

z
_
t

3
t

2
_
4C
+I
0a

CD
h

_
_
ln
_
_
r
2
w
+h
2
+h

2
(h

/2)
_
_
r
2
w
+h
2
h

2
(h

/2)
+ O(r
2
w
)
_
_
, (47)
where I
0a
is calculated at the appropriate time (t

1
for maximum and t

3
for minimum). The
rst cycle extreme pressures can be calculated simply from the omission of I
0a
. As seen, the
Laplace transform method provides simple and useful expressions of the well pressure.
The above expressions were obtained for conditions of negligible transient effects in the
vertical direction. In order to assess if those transient effects are indeed negligible based
on their duration, well pressure variations during the rst cycle were calculated for typical
operating conditions and reservoir characteristics, as seen in Fig. 3a and b. The gures pres-
ent a comparison between the exact solution (of the linear problem, Eqs. 36 or 39) and the
simplied solution (Eq. 42, omitting I
0a
). All the stages were shifted in time to represent
them on a single time scale (t

). As it turns out, already at

2

z
t

= 2, the solutions are

essentially identical (with proximity smaller than 0.02%). That is, for a cycle time of 24h
and
z
= 100, 3min after the beginning of each stage, vertical pressure transients terminate.
In order to estimate the error entailed by the linearization of Eq. 22 (obtained by setting
p

= 1), the analysis described in the appendix of Kushnir et al. (2008) was incorporated.
The results conrm that the nonlinearity effects of Eq. 22 are indeed negligible. It is worth
noting that the linearized modied pressure representation allows the extension of the sin-
gle well solution to solutions of multiple well elds, simply by superposition. Likewise, for
nite reservoir boundaries, solutions can be derived from the superposition of appropriate
well images.
1 3
232 R. Kushnir et al.
4 Interface Location and Stability
In principle, knowledge of the reservoir pressure distribution enables the calculation of the
interface curve z

(r

, t

) during discharge based on Eq. 29. The highest interface rise

occurs right beneath the well, approximately at the end of the withdrawal stage (r

= 0, t

=
t

3
). Based on the pressure distribution solution, the modied pressure beneath the well at
that time is
2

s
m

= 0
t

= t

3
= ln
t

3
t

3
t

1
CDln

z
_
t

3
t

2
_
4C
+ I
0a

CD
h

ln

2
((z

) /2) sin ( (z

) /2)

2
((z

+h

) /2) sin ( (z

+h

) /2)
. (48)
The substitution of this expression into Eq. 29 yields an implicit equation from which the
extent of the highest interface rise can be determined (subject to the set of pertinent param-
eters: m

c
, h

, g

,
z
, t

1
, t

2
, t

3
).
In order to demonstrate how the location and stability of the interface can be determined,
the air pressure distribution beneath the well at the end of the rst withdrawal stage is plotted
in Fig. 4a and b, for different values of m

c
s and penetration depths, respectively. The linear
dotted line in the gure refers to the water hydrostatic pressure distribution beneath the well.
The intersection points of the dotted line with the air pressure curves represent the solutions
of Eq. 29. As seen from Fig. 4a, for m

c
= 0.006, there are two solutions (points A and B).
At point B, the pressure gradient in the air zone exceeds the hydrostatic gradient in the water
zone. Hence, the height at B is unstable, and water could be drawn into the well. Conversely,
point A pertains to a physically stable height, as the air pressure gradient is smaller than that
of the water zone. As m

c
is increased, the interface height rises up to a critical height (point
C). At that point, the pressure gradient at both sides of the interface is identical, and any
increase in m

c
will drive the water into the well. Therefore, the interface height at point C is
dened as critical, and the corresponding value of m

c
(0.01255) is the critical suction rate.
Beyond that critical rate, a stable interface height cannot exist, as seen from the inspection of
the curve m

c
= 0.024. Similarly, for any given value of m

c
, there exists a critical penetration
depth, as seen in Fig. 4b.
Following the above, for any given compressed air pressure distribution, the airwater
interface stability can be determined. Evidently, the critical mass ow rate and interface rise
are dened by two conditions stating that both the pressure and the vertical pressure gradi-
ent on each side of the interface are identical. Hence, the equations determining the critical
conditions are
r

= 0, z

, t

= t

3
_
1 +

_
1/2
= 1 g

_
1 z

_
(49)
r

= 0, z

, t

= t

= 2g

_
1 g

_
1 z

__
. (50)
Given a set of parameters (h

, g

,
z
, t

1
, t

2
, t

3
), the critical dimensionless mass ow rate,
m

c cr
, and the critical interface rise can thereby be determined.
In reality, critical conditions can involve appreciable interface rises, whereas the derived
solutions are for relatively small interface uctuations. Nonetheless, Muskat (1937) showed
that such an approximation will only affect the quantitative details of the results. It does
not invalidate the general conclusions relating to the existence of critical conditions and
the dependence of the critical ow on its pertinent parameters. Indeed, comparison of the
1 3
Compressed Air Flow of CAES Plants 233
(a) (b)
Fig. 4 Air and water pressure distributions beneath the well at the termination of the rst discharge stage
(r

= 0, t

= t

3
) for
z
= 78.75, g

= 0.03924, t

1
= 10/24, t

2
= 12/24, and t

3
= 22/24. a At different
m

c
s and b at different h

s
Muskat and Wyckoff (1935) approximate solution with solutions that addresses the interface
shape (Wheatley 1985; Hoyland et al. 1989), reveals that the approximation entails somewhat
higher critical ow rates but represents quite accurately the inuence of the relevant param-
eters (well screen length, layer height, etc.) on the critical conditions. Furthermore, from
a quantitative point of view, the approximate critical ow rate is sufcient to estimate the
number of wells required for any desired application, and thereby determines the economic
feasibility.
5 Results and Discussion
5.1 General Considerations
Table1 provides representative ranges of the studied reservoir characteristics, operating con-
ditions, and their corresponding dimensionless parameters. As seen in the table, the values
of r

w
are indeed small enough to justify the omission of terms containing r
2
w
, and terms of
higher power, in Eq. 47. These terms can be significant only if k
z
is much larger than k
r
(the
maximumvalue of r

w
in the table is based on a permeability ratio of unity). In addition, it was
shown in Fig. 3 that vertical pressure transients terminate when the time elapsed at each stage,
from its initiation, is approximately twice the vertical characteristic time, H
2
/(
z

2
). That
is, even for the maximum vertical characteristic time of 2,700s (Table1), vertical pressure
transients terminate way before the end of each stage.
In order to closely inspect the nature of the pressure distribution, calculated steady peri-
odic pressure distributions for a cycle period at different radii and elevations are illustrated
in Fig. 5a and b, respectively (for the indicated set of operating conditions). As seen, the well
pressure undergoes rapid changes at the beginning of each time interval (owing to the ow
variations). Following those sharp responses, the rate of change of the pressure moderates, as
it approaches stable conditions. The well pressure reaches its crest at the end of the charging
stage and drops to its trough at the termination of the power generation stage. Figure5a reveals
a diminution of the pressure amplitude with a slightly progressive phase lag as r

increases.
As observed, the penetration of a well pressure uctuation into the reservoir does not exceed
1 3
234 R. Kushnir et al.
Table 1 Representative ranges of reservoir characteristics, operating conditions, and their corresponding
dimensionless parameters (The maximal value of r

w
is based on the maximal permeability ratio of 1, and
separately, the specic gas constant is R = 0.287kJ/kg K)
Variable Definitions Minimum value Maximum value Units
f Porosity 0.05 0.35 unitless
k
r
Radial permeability 100 5,000 md
k
z
Vertical permeability 100 5,000 md
H Layer height 5 25 m
r
w
Well radius 0.05 0.6 m
P
0
Initial air pressure 2 7 MPa
m
c
Compressor ow rate 1 50 kg/s
T Air temperature 300 400 K
Air viscosity 1.8 10
5
2.4 10
5
kg/m s

w
Water density 930 1,000 kg/m
3
t
p
Time period 24 24 h
Z Air compressibility factor 0.99 1.01 unitless

r
,
z
Pneumatic diffusivities 0.02 40 m
2
/s
H
2
/(
z

2
) Vertical characteristic time 0.07 2.7 10
3
s

z
t
p

z
/H
2
3 1.3 10
5
r

w
(k
z
/k
r
)
1/2
r
w
/H 3 10
4
0.12
m

c
m
c
ZRT/(H k
r
P
2
0
) 8 10
5
20
g

w
gH/P
0
7 10
3
0.12
t

1
t
1
/t
p
6/24 12/24
t

2
t

1
(t
2
t
1
)/t
p
2/24 8/24
t

3
t

2
(t
3
t
2
)/t
p
2/24 10/24
a certain distance. Beyond that distance, the pressure oscillations vanish. Hence, the innite
reservoir approach is valid when the effective reservoir radius exceeds the uctuating gas
penetration radius. For a fully penetrating well, the penetration radius is approximately equal
to (t
p

r
)
1/2
(Kushnir et al. 2008). For a partially penetrating well, the ratio of the maximum
pressure (at any radius) versus the maximum well pressure is likely to drop within a shorter
radius than the fully penetrating well (owning to a higher maximum well pressure). As seen
in Fig. 5b, proceeding vertically downward, the pressure amplitude decreases along the well,
and the phase lag is practically negligible. Based on the latter, calculations of the critical
conditions can be performed at the end of each withdrawal stage, representing the time of
minimal pressure in the gas layer below the well.
Since previous investigations addressed only the rst cycle solution, it is of essence to
compare the extreme conditions of the rst cycle and periodically steady-state solutions.
For that purpose, the relative error ( p

f
p

s
)/p

s
was calculated and plotted (including the
pressure) against m

c
, for the maximum and minimum pressure values, as seen in Fig. 6 a
and b, respectively. For the particular calculations shown in Fig. 6, the relative error is indeed
negligible for moderate pressure uctuations. It becomes perceptible only in cases of imprac-
tical high pressure uctuations. It is important to note that applicable ranges of operating
conditions are to avoid interface instability, and thus are bounded by m

c cr
(here, in particular,
m

c cr
= 0.0163 for g

= 0.12). Similar calculations for the entire range of the parameters

appearing in Table1 lead to the same conclusion. Therefore, it seems that for all practical
conditions, subject to the critical conditions ( m

c
m

c cr
)and within the limitations of mod-
erate pressure uctuations ( p

s min
0.75), the relative error can be considered negligible. In
conclusion, for normal operating conditions, the rst cycle solution adequately represents the
steady periodic solution. Hence, to obtain a cycle period pressure distribution, it is sufcient
1 3
Compressed Air Flow of CAES Plants 235
(a) (b)
Fig. 5 Steady periodic dimensionless pressure oscillations for t

1
= 7/24, t

2
= 14/24, and t

3
= 18/24.
a At different dimensionless radii and b at different dimensionless depths
to know the solution for a constant well mass ow rate. Any rst cycle pressure distribution
(for a daily or a weekly cycle) can be constructed from that solution.
5.2 Parametric Study
For the design of a CAES plant, predictions of the air pressure oscillations in the reservoir
during operation are of interest, in particular, as it pertains to the oscillations extreme values.
In general, an increase in pressure oscillations entails larger compression and lower turbine
output, both, reducing the plant efciency. Clearly, to mitigate pressure uctuations and
associated losses, the well screen length must be extended as much as possible. On the other
hand, larger penetrations run the risk of undesired water suction. In order to demonstrate the
effects of the penetration of the well on pressure uctuations, the range of pressure oscillation
p

s
= p

s max
p

s min
is plotted as a function of h

in Fig. 7ac, for several values of m

c
, r

w
,
and
z
. It is seen that p

s
is highly dependent on h

for a range of small well penetrations.

This range widens along with larger p

s
, as m

c
increases or r

w
decreases. At longer well
penetrations, the p

s
variation moderates.
Expectedly, Fig. 7ac reveals that p

s
increases with both m

c
and
z
, and decreases with
r

w
. From these parameters, m

c
is the most inuencing. Its dominance can be observed in
Eq. 47 as well, since the modied pressure is directly proportional to m

c
. The parameter
z
appears in a logarithmic term and therefore affects mildly p

s
. As seen from, both, Fig. 7b
and Eq. 47, the inuence of r

w
on p

s
is strongly dependent on the well penetration. At
short penetrations, p

s
is strongly effected by r

w
, while at deep penetrations, p

s
becomes
independent of r

w
. This assertion stems from the fact that at high well penetrations, the ow
is radial and thus depends on
z
/r
2
w
(
r
), while at short well penetration, vertical ow is
appreciable and varies with r

w
. Notice that the plots do not address instability limitations
of the ow, and as such are representative only for the corresponding ranges of stable ow
conditions.
Beyondthe considerationof pressure uctuations, the CAESplant must be designedwithin
the limitation of critical ow rates. As aforementioned, the critical dimensionless mass ow
rate and interface rise are calculated by Eqs. 49 and 50. Results of such calculations are seen in
Fig. 8, showing the critical conditions dependence on the relative well penetration for various
values of g

(8a) and
z
(8b). The predicted trends agree with the results of previous inves-
tigations of steady oil wells critical conditions (e.g., Wheatley 1985; Hoyland et al. 1989).
1 3
236 R. Kushnir et al.
(a)
(b)
Fig. 6 Relative error (between the steady periodic and rst period pressures) as a function of m

c
for different
values of r

w
. a At the cycle maximum pressure and b at the cycle minimum pressure
Evidently, critical ow rates are smaller for deeper well penetrations, with a markedly higher
sensitivity at deep penetrations. On the other hand, the relative interface critical rise is larger
at deeper well penetrations (the relative interface rise is dened as y

= (H )/(H h)).
It is seen that m

c cr
increases with g

and decreases with

z
. As it turns out, the critical
ow rate strongly depends on the gas layer height, where a mild extension of H entails a
significant increase of m
c cr
. In addition, the critical ow rate is directly proportional to the
horizontal permeability and approximately inversely proportional to the temperature (effects
of viscosity and compressibility owing to temperature changes are negligible). This is not a
consequence of the linearization and applies also to the exact solution of Eqs. 12 through 18
(since k
r
and T appear only in m

c
). The vertical permeability and the porosity, which appear
solely in
z
, affect little the critical ow rate.
In general, the variables that can mitigate the pressure uctuations also stabilize the inter-
face. Excluded from such variables are the well screen length and the vertical permeability.
1 3
Compressed Air Flow of CAES Plants 237
(a) (b)
(c)
Fig. 7 Ranges of dimensionless pressure variations versus relative well penetrations. a At different m

c
s,
b at different r

w
s, and c at different
z
s
(a) (b)
Fig. 8 Critical dimensionless ow rates and relative interface rises versus relative well penetrations. a At
different g

s and b at different
z
s
This is manifested in Figs. 7 and 8 where larger well penetrations reduce both the pressure
uctuations and the critical ow rate (notice that an increase in vertical permeability can be
interpreted as a well penetration elongation). A decrease in pressure oscillations is desirable
because it reduces the compression work and enlarges the turbine output. However, lower
critical ow rates would require more wells. In order to elucidate these opposing effects, the
1 3
238 R. Kushnir et al.
Fig. 9 Effects of the well penetration on the critical dimensionless ow rate and the corresponding range of
pressure variations
critical ow rate and the corresponding range of pressure oscillations (calculated at m

c cr
)
are plotted in Fig. 9. Obviously, the choice of the preferred screen length must be strongly
based on economics.
It is warranted to extend the discussion and address quantitatively real critical ow rates.
Referring to the Huntorf CAES plant, where the compressor works 12h at 108kg/s and the
turbine works 3h (producing 290MW) with an inlet pressure of 43 bar (Crotogino et al.
2001). Though the plant operates with salt cavern reservoirs, the same operating conditions
are subsequently used to determine how many wells are required to perform similarly with
an aquifer reservoir. The maximum allowable ow rate, per well, is therefore calculated
for t

1
= 12/24, t

2
= 18/24, t

3
= 21/24, and P
0
= 45 bar. Assuming an aquifer with
f = 0.2, k
r
= 500 md, k
z
= 400 md, and T=310K (Z=0.9946, = 1.964 10
5
Pas,

w
= 995 kg/m
3
), and then for H = 25m, the highest critical owrate (h

= 0) is 2.34kg/s.
It means that for a total of 108kg/s, at least 46 wells are needed (excluding well interference).
Accordingly, for H = 15 m, m
c cr
= 0.724 kg/s (149 wells), and for H = 5 m, m
c cr
=
0.0627 kg/s (1,722 wells). The calculation demonstrates that water coning could constitute
a severe limitation, especially for large scale plants. Consequently, it is apparent that CAES
designs must give preference to the smallest practical m

c
(e.g., larger H, and k
r
or smaller
T) for the reduction of both number of wells and pressure uctuations.
All the above results are based on xed time intervals (charging, storage, and power gener-
ation). Realistic bounds for these intervals are shown in Table1. Basically, the time intervals
are determined so as to meet the local power demand and production capacity. Effects of
the duration period of the power generation on the critical ow rate and the corresponding
range of minimum pressure are shown in Fig. 10. As seen, prolonging the duration of power
generation can significantly increase the well storage capacity. Moreover, the corresponding
minimum pressure is even slightly higher owing to smaller discharge ow rates. As for the
compression phase, referring to the bounds shown in Table1, prolongation of its duration
period would not affect the well storage capacity (i.e., m

c cr
t

1
const), but would reduce
p

max
owing to a corresponding smaller ow rate. Consequently, it is desirable to expand the
duration periods of compression and power generation, as much as feasible.
1 3
Compressed Air Flow of CAES Plants 239
Fig. 10 Dimensionless critical ow rates and corresponding minimum pressures versus well penetration for
different power-generation-duration periods
6 Conclusions
A theoretical investigation of a compressible gas ow within CAES plant aquifer reservoirs
was conducted. The analysis provides simple and useful expressions for the periodic air
pressure distribution. Based on these expressions both the well pressure and the stability
conditions of the airwater interface are determined. The following conclusions were drawn
from the investigation:
For normal operating conditions (characterized by moderate pressure uctuations and
smaller than critical mass ow rates), the rst cycle solution adequately represents the
steady periodic solution.
Water coning could impose severe limitations (especially for large scale plants) on the
discharge ow rates.
CAES plant designs must give preference to the smallest practical m

c
(e.g., larger H, and
k
r
or smaller T) for the reduction of both the number of wells and the pressure uctuations.
It is desirable to expand the compression and power generation duration periods, as much
as feasible. Prolongation of the duration period of the power generation, in particular, can
significantly increase the well storage capacity.
The choice of the well screen length must be based on economics (accounting for the
plant energy capacity and the number of required wells).
The analytical solution developed in this study can also be used to construct a solution for
multiple well systems. It provides an important tool that could eventually support optimiza-
tion analyses of compressed air storage.
References
Allen, R.D., Doherty, T.J., Schainker, R.B., Istvan, J.A., Pereira, J.C.: Preliminary results from the pittseld
aquifer eld test applicable to commercialization of CAES technology. Intersociety Energy Conversion
Engineering Conference, San Francisco, USA, pp. 10811090 (1984)
ANRStorage Company: compressed air energy storage in porous media, EPRI Report 2488-10, March (1986)
1 3
240 R. Kushnir et al.
Braester, C., Bear, J.: Some hydrodynamics aspects of compressed-air energy storage in aquifers. J. Hydrol.
(Amst.) 73, 201225 (1984)
Carslaw, H.S., Jaeger, J.C.: Conduction of Heat in Solids, 2nd edn. Oxford University Press, Oxford (1959)
Crotogino, F., Mohmeyer, K.U., Scharf, R.: Huntorf CAES: more than 20 years of successful operation. SMRI
Spring Meeting, Orlando, USA, pp. 351357 (2001)
Hoyland, L.A., Papatzacos, P., Skjaeveland, S.M.: Critical rate for water coning: correlation and analytical
solution. SPE Reserv. Eng. 4(4), 495502 (1989). doi:10.2118/15855-PA
Kushnir, R., Ullmann, A., Dayan, A.: Steady periodic gas ow around a well of a CAES plant. Transp. Porous
Media 73(1), 120 (2008). doi:10.1007/s11242-007-9156-x
Meiri, D., Karadi, G.M.: Simulation of air storage aquifer by nite element model. Int. J. Numer. Anal. Meth.
Geomech. 6(3), 339351 (1982). doi:10.1002/nag.1610060306
Muskat, M.: The Flowof Homogeneous Fluid through Porous Media. 1st edn. McGraw-Hill, NewYork (1937)
Muskat, M., Wyckoff, R.D.: An approximate theory of water coning in oil production. Trans. AIME 114, 144
163 (1935)
Ruud, N.C., Kabala, Z.J.: Response of a partially penetrating well in a heterogeneous aquifer: integrated well-
face ux vs. uniform well-face ux boundary conditions. J. Hydrol. (Amst.) 194, 7694 (1997). doi:10.
1016/S0022-1694(96)03217-9
Wheatley, M.J.: An approximate theory of oil/water coning. Paper SPE 14210, SPE 60th annual technical
conference and exhibition, Las Vegas, USA (1985)
Wiles, L.E.: Numerical analysis of temperature and ow effects in a dry, two-dimensional, porous media res-
ervoir used for compressed air energy storage, technical report PNL-3047, Pacic Northwest Laboratory
(1979)
Wiles, L.E., McCann, R.A.: Water coning in porous media reservoirs for compressed air energy storage,
technical report PNL-3470, Pacic Northwest Laboratory (1981)
1 3