Journal of Petroleum Science and Engineering 127 (2015) 179–189

Contents lists available at ScienceDirect

Journal of Petroleum Science and Engineering

journal homepage: www.elsevier.com/locate/petrol

Combustion enhance recovery of shale gas

Grigori Chapiro a,n, Johannes Bruining b
a
Departamento de Matemática, Universidade Federal de Juiz de Fora, 36036-330, Juiz de Fora – MG, Brazil
b
Department of Geoscience & Engineering, Delft University of Technology, Stevinweg 1, 2628CN Delft, The Netherlands

art ic l e i nf o a b s t r a c t

Article history: The amount of extracted natural gas from low permeable organic rich shale has increased rapidly over
Received 28 April 2014 the past decade. The technologies used to improve the permeability of a shale gas reservoir are
Accepted 26 January 2015 horizontal drilling and hydraulic fracturing. The potential environmental impact of hydraulic fracturing
Available online 4 February 2015
motivates research aiming at alternative permeability enhancement methods. In this paper the
Keywords: possibility of in-situ combustion to improve permeability is investigated. Two possibilities for the in-
Gas shales situ fuel source are considered, viz., methane or kerogen. A one dimensional model is considered. Under
in-situ combustion simplifying hypotheses a quasi-analytical solution for the corresponding Riemann problem is obtained.
partial differential equations The solution was analysed with parameters that correspond to resonance conditions, i.e., conditions for
resonance conditions
which the speeds of the heat wave and the combustion wave are equal. We conclude that methane
combustion cannot generate enough heat to enhance the permeability. However kerogen, if present in
sufficient quantities, makes this possible. We present the set of parameters for which the combustion
reaches the optimal temperature to improve the permeability. Finally, the analytical approximations are
validated with direct numerical simulations. The pressure distribution in the reservoir and the
production improvement by the thermal treatment are calculated with the numerical model.
& 2015 Elsevier B.V. All rights reserved.

1. Introduction One possible alternative method is combustion. The literature

on chemical processes occurring during in-situ combustion (ISC)
The production of natural gas from shale gas reservoirs has incr- has many analogies with combustion in gas shales.
eased over the past few years and become an important component The division of the ISC processes from the production to the injec-
of the world's energy supply, (Kuuskraa et al., 2011). A shale gas tion well can be schematically separated into different zones: the
reservoir is characterized by an extremely low matrix permeability, light hydrocarbons zone, the steam zone, the cracking zone, the coke
clusters of mineral-filled “natural” fractures and an organic-rich formation zone and the coke combustion zone, (Gates and Ramey,
deposition (Zhang et al., 2009; Busch et al., 2008). The current state 1958). According to Abu-Khamsin et al. (1988) there are three over-
of the art technology to obtain economical production from a shale lapping processes caused by temperature elevation due to heat tran-
gas reservoir uses horizontal drilling and hydraulic fracturing, which sfer downstream the reaction zone. The first process is distillation
is used to improve the permeability of the reservoir, (Kennedy et al., (evaporation) of light oil components, which occurs at lower tem-
2012). As a result of hydraulic fracturing a slit is created, which can peratures. Next, we separate the “visbreaking” stage that occurs at
significantly extend the effective radius of the well, (Hagoort, 1981; medium temperatures and breaks the initial oil into low volatile
Prats, 1961). The potential environmental impact of hydraulic fractur- components and a “visbroken” oil (with smaller viscosity and specific
ing possibly includes the contamination of ground water, risks to air gravity). Finally, the severe cracking (coking) reaction transforms the
quality, the migration of gases and the migration hydraulic fracturing visbroken oil into coke and some lighter components.
chemicals to the surface, etc, (Howarth et al., 2011; Birol, 2011). These In this text we study the corresponding processes in the com-
hazards are of particular concern in densely populated areas like bustion zone and in the light hydrocarbon zone. Inside the com-
Europe. This motivates research aiming at alternative methods to bustion zone the residual coke reacts with the injected oxygen. The
improve the permeability of tight reservoirs. micro-porous structure of coke and related parameters, such as
reactive area of coke particles and transport of oxygen into the
pores, are important to describe this process, see Bird et al. (2002)
n and Levenspiel (1999). There is abundant literature, for example
Corresponding author. Tel.: þ 55 32 21023308 (Office).
E-mail addresses: [email protected] (G. Chapiro), (Essenhigh, 1981; Smith, 1978), focusing on the reaction kinetics
[email protected] (J. Bruining). for general materials containing carbon. The combustion zone is

http://dx.doi.org/10.1016/j.petrol.2015.01.036
0920-4105/& 2015 Elsevier B.V. All rights reserved.
180 G. Chapiro, J. Bruining / Journal of Petroleum Science and Engineering 127 (2015) 179–189

typically tens of feet wide and the temperature of the combustion of methane present in the matrix or (2) the fuel consists of the
front is usually around 430–540 1C. Condensation of all light oil immobile part of heat treated kerogen, see below.
that has been produced elsewhere occurs inside the light hydro- Gas shales contain methane of three different types, (Glorioso et
carbon zone. In some cases, oxygen can break through the coke al., 2012): methane in the matrix pores, methane in the kerogen
zone and come into contact with the light hydrocarbons where the pores and methane adsorbed at the surfaces of shale and kerogen.
reaction is described by “Low Temperature Oxidation” (LTO) mec- When the temperature rises the adsorbed methane can be liber-
hanism, see Cinar et al. (2009). ated and react with oxygen.
In the case of methane combustion, the separation in zones is A shale gas reservoir is a source rock that over geologic times
much simpler (Bubnovich et al., 2006), as there will be no visbreak- retained hydrocarbons, which comprises the total organic carbon
ing process and no hydrocarbon vaporization. Therefore we study in (TOC) currently stored inside the matrix. The soluble part of TOC is
this paper two types of fuel, viz., (1) the immobile part of the ker- classified as bitumen and the insoluble part as kerogen. The TOC
ogen, which is analogous to coke combustion and (2) methane content is an important property of the reservoir and is the subject
present in the reservoir. of many studies, e.g., (Chelini et al., 2010; Akkutlu and Fathi, 2012).
Shales permeabilities vary between the micro-Darcy range and The TOC content present in the rock can be poor (less than 0.5%),
the nano-Darcy range depending on the depth of burial see, e.g., fair (0.5–1%), good (1–2%), very good (2–4%) and excellent (more
(Soeder, 1988). The permeability of gas shales is affected by the than 4%). In extreme cases it can reach 20% and more, see Chelini
organic rich deposition (kerogen) (Allan and Mavko, 2012; Akkutlu et al. (2010). When subject to intense heat some types of kerogen
and Fathi, 2011, 2012; Fathi and Akkutlu, 2012) and removal by release crude oil or natural gas as shown in Fig. 1, see Seewald
combustion of the kerogen might enhance the permeability. (2003) for details. This fact triggers the interest in in-situ combus-
It is asserted that heat treatment of gas-shale reservoirs can imp- tion, as it can induce the liberation and production of natural gas.
rove the permeability in the same way as for oil shales, (Sanmiguel Similarly to what happens during in-situ combustion (Gates and
et al., 2002; Busch and Amann-Hildenbrand, 2013). Thermal treatment Ramey, 1958) part of the kerogen (about 10%) stays immobile in a
enhances the permeability in tight gas reservoirs by vaporizing the form of heavy residual oil after heat treatment and can be used for
capillarily blocked water, dehydrating the clay-bound water, destroy- combustion (Perry and Lee, 2007).
ing the clay lattice and creating thermally induced microfractures, The main goal of this paper is to understand under which con-
(Jamaluddin et al., 2000). In their paper, Jamaluddin et al. (2000) ditions it is possible to increase the permeability of the gas shales
obtain an increment in permeability from 50% to 675% depending on using thermal treatment by in-situ combustion. We solve the system
the core type and the temperatures to which these cores were of partial differential equations as a Riemann problem obtaining the
exposed. The maximum increment was obtained for a temperature solution in a form of a sequence of waves, viz., contact discontinuities
around 650 1C. In this paper the possibility whether a combustion and a combustion traveling wave. The last wave is studied following
process in a reservoir at resonance conditions can be used to improve the technique used in Bruining et al. (2009), Aldushin et al. (1999)
the permeability is investigated. At resonance conditions the heated and Mailybaev et al. (2011) and cited references therein. This
zone essentially stays of constant size. Thus the permeability is technique uses the separation of the combustion wave into low
improved by heating the reservoir only locally avoiding the use of temperature (pre-reaction) and high temperature (reaction) layers.
excessive energy. This paper is organized as follows. Section 2 presents a math-
In order to produce such a combustion wave an in-situ fuel is ematical model of the in-situ combustion process. Section 3
needed. Two possibilities are considered, viz., (1) the fuel consists analyzes the case of methane combustion. Section 4 studies the

2.0

1.5 T 50ºC
Other products
Increasing time and temperature
Atomic H/C ratio

1.0

Oil

0.5
Natural gas

T 160ºC

0.0
0.00 0.05 0.10 0.15 0.20 0.25
Atomic O/C ratio Relative yield

Fig. 1. (a) Modified Van Krevelen diagram indicating carbon, hydrogen and oxygen ratios during the evolution of the immature kerogen of varying compositions (Types I, II,
III and IV), see Seewald (2003). (b) Organic products generated during progressive burial in sedimentary basins. The form of this figure is constrained by the maturation
trends shown in the Van Krevelen diagram (see (a)).
 G. Chapiro, J. Bruining / Journal of Petroleum Science and Engineering 127 (2015) 179–189 181

combustion of the immobile part of kerogen. Section 5 presents will therefore be liberated slowly. We will use RY ¼ 0 and 0 o RZ o1
some numerical simulations. Section 6 presents the possible effect in Section 3.
of the combustion on the shale gas production. Section 7 ends with
Next dimensionless variables are introduced
some final remarks and conclusions.
T  T res Y Z ρk x t
θ¼ ; Y¼ ; Z¼ ; ρk ¼ ; x¼ ; t¼ ; ð7Þ
Tn Yn Zn ρnk xn tn
2. Modeling
where the thermometric conductivity αs ¼ D=C m with C m ¼
ð1  φÞcm ρm was introduced following Landau and Lifshitz (1987,
An infinite cylindrical gas shale core is considered where one
Section 50). The reference quantities are
dimensional flow is in the axial x direction. Mixture of nitrogen and
oxygen is injected from the left into a homogeneous porous αs xn
xn ¼ ; tn ¼ ; ρnk ¼ ρres
k ; T n ¼ T res ; Y n ¼ Y inj ; Z n ¼ Z res :
medium initially containing methane and an immobile fraction of ug ug
kerogen. ð8Þ
Gas diffusion is disregarded. Local thermal equilibrium between
In dimensionless form (see Chapiro et al., 2012; Landau and
the fluid and the porous medium is considered. In order to permit
Lifshitz, 1987 for similar calculations) the system (1)–(5) becomes
rigorous identification of traveling waves, the model is simplified by
ignoring the dependence of the gas density on temperature, pre- ∂θ ∂θ ∂2 θ
þvθ ¼ 2 þ μθm Φm þ μθk Φk ; ð9Þ
ssure and composition, i.e. the gas density is constant. The model ∂t ∂x ∂x
used here is based on a model proposed in Chapiro et al. (2012). The
∂Y ∂Y
small variations in pressure P due to the flow of gas are disregarded þ vY ¼  μYm Φm  μYk Φk ; ð10Þ
resulting in a constant gas speed ug. Therefore the pressure P is ∂t ∂x
considered constant. Moreover the heat capacity of gas is neglected ∂Z ∂Z
with respect to the heat capacity of the formation. The model consists þ vZ ¼ μm Φm ; ð11Þ
∂t ∂x
of a heat balance (1), an oxygen-, a methane- and an immobile ker-
ogen mass balance (2)–(4). The reaction rate is given by the linear law ∂ρk
¼  μk Φk ; ð12Þ
of mass action, see Bird et al. (2002); Eq. (5) for methane and Eq. (6) ∂t
for immobile kerogen. The temperature dependence for both reac- where dimensionless reaction rates for methane and kerogen are
tions is given by the Arrhenius law. The primary variables are: the given by
temperature T (K), the partial pressures of oxygen Y ðÞ and methane

 Em
Z ðÞ, the scaled immobile fuel concentration ρk ðÞ obtained by Φm ¼ YZ exp ; ð13Þ
θþ1
dividing the current concentration of fuel ρ k by the initial concentra-
tion of fuel ρres
k , where ðÞ indicates a dependence on dimensionless


E k
variables. The combustion products of which the concentration is Φk ¼ Yρk exp : ð14Þ
θþ1
denoted by X satisfies the equation X ¼ 1  Y  Z. The equations read
Here vθ , vY and vZ represent the dimensionless speeds of the
∂T ∂T ∂2 T thermal, oxygen and methane waves; E m and E k are the dimen-
ð1  φÞcm ρm þ cg ρg ug ¼ D 2 þ Q m W m þ Q k W k ; ð1Þ
∂t ∂x ∂x sionless activation energies; μθm , μm and μYm are the dimensionless
pre-exponential coefficients for methane combustion representing
∂Y P ∂Y P the generated heat, consumed fuel and oxygen respectively; μθk , μk
φð1 þ RY Þ þ ug ¼  kYm W m  kYk W k ; ð2Þ
∂t P 0 ∂x P 0 and μYk represents the dimensionless pre-exponential coefficients
for kerogen combustion representing the generated heat, con-
∂Z P ∂Z P sumed fuel and oxygen respectively. These quantities are defined
φð1 þ RZ Þ þ ug ¼ km W m ; ð3Þ
∂t P 0 ∂x P 0 as follows:

∂ρk cg ρg 1 1 Em Ek
¼ kk W k ; ð4Þ vθ ¼ ; vY ¼ ; vZ ¼ ; Em ¼ ; Ek ¼ ;
∂t Cm φð1 þ RY Þ φð1 þ RZ Þ RT n RT n

where reaction rates for methane and kerogen are respectively Q m kpm t n Y n Z n P Q k kpk t n Y n ρnk P

 μθm ¼ ; μθk ¼ ;
P  Em Cm T n P0 Cm T n P0
W m ¼ kpm Y Z exp ; ð5Þ kk kpk Y n t n P
P0 RT km kpm Y n t n P
μm ¼ ;μ ¼ ;
φð1 þ RZ Þ P 0 k φ P0


P  Ek kYm kpm Z n t n kYk kpk ρnk t n
W k ¼ kpk ρk Y exp : ð6Þ μYm ¼ ; μYk ¼ : ð15Þ
P0 RT φð1 þ RY Þ φð1 þ RY Þ
All parameters are given in Table C1. For the parameter values presented in Table C1, the nondimen-
sional constants given in (15) are
Remark. We prefer to give a highly simplified model of the sor-
ption behavior of the gases present, i.e., oxygen, methane and vθ ¼ 9:9  10  4 ; vY ¼ 5; vZ ¼ 2:5; E m ¼ 52:1; E k ¼ 23:3;
combustion products. The free oxygen Y and free methane Z are μθm ¼ 4:95  109 ; μθk ¼ 3:57  106 ; μm ¼ 7:42  1011 ;
convected with the gas stream, react with each other and sorb on
μk ¼ 2:85  106 ; μYm ¼ 1:41  1013 ; μYk ¼ 3:40  108 : ð16Þ
the medium proportional to their respective concentrations, i.e.,
Y s ¼ RY Y and Z s ¼ RZ Z . These coefficients enter in the accumulation The system (9)–(14) has to be solved on the interval ½0; 1Þ with left
terms in Eqs. (2) and (3). The parameters RY and RZ are considered and right Dirichlet boundary conditions given respectively by
to be independent. We always assume that RY o RZ , otherwise θð0; tÞ ¼ 0; Yð0; tÞ ¼ 1; Zð0; tÞ ¼ 0; ρk ð0; tÞ ¼ 0;
there is no overlap between oxygen and methane and therefore no
θð1; tÞ ¼ 0; Yð1; tÞ ¼ 0; Zð1; tÞ ¼ 1; ρk ð1; tÞ ¼ 1; ð17Þ
reaction occurs. It is also expected that RY o RZ because methane
has penetrated the nanopores of shales during geologic times and and initial conditions:
182 G. Chapiro, J. Bruining / Journal of Petroleum Science and Engineering 127 (2015) 179–189

( (
1; 0 ox oxc ; 0; 0 ox oxb ; Table 1
Yðx; 0Þ ¼ Zðx; 0Þ ¼ Dimensional parameters used in recovery rate computation and their typical
0; xc o x; 1; xb ox;
( values.
0; 0 o x o xc ;
ρk ðx; 0Þ ¼ ð18Þ Symbol Physical quantity Value Unit
1; xc o x:
 
ct ρg Spec. total compressibility ct ρg ¼ M g =RT res 0.0112 g/J
The initial condition for the temperature is given by Pn
8
Reference pressure (1 atm) 1  105 Pa

< 0; 0 o x o xc ;
> Pres Initial pressure inside p.m. (100 atm) 1  107 Pa

θðx; 0Þ ¼ 1; xc o x o xb ; ð19Þ Phole Pressure inside the hole (10 atm) 1  106 Pa
>
: 0; x o x: Lres Reservoir length 1000 m
b
Rres Reservoir maximum radius 10 m
Mg Average gas (methane) molar weight 16 g/mol
In order to study the principal contribution to the heat genera-
kini Initial reservoir permeability (100 nD) 10  19 m2
tion in the following sections two situations are analyzed: (1) the μg Average gas viscosity 1:85  10  5 Pa s
total amount of immobile kerogen is negligible and (2) the com-
bustion is dominated by the immobile part of kerogen. In this text
the attention is focused on the combustion wave, because the other impossible to produce a high temperature traveling wave. Conse-
waves are simple contact discontinuities (Chapiro et al., 2012). quently the thermal permeability enhancement in the absence of
kerogen involves heating the entire reservoir and is prohibitively
expensive.
3. Methane combustion
The combustion wave speed is proportional to the gas velocity
as the displacement front between injected gas and methane is
Let us assume that the amount of immobile kerogen is negli-
moving along with the gas speed. The thermal wave speed is
gible. Thus the system (9)–(14) can be reduced to (9), (10), (11),
always lower than the gas velocity by a factor given by the ratio of
(13) with zero combustion rate of kerogen, Φk  0.
the gas heat capacity and the rock heat capacity, which is
The combustion process occurring in the reaction layer and the
practically smaller than one. As the pressure increases the heat
pre-reaction layer is analyzed by following the method used in
capacity of the gas increases and thus the ratio between gas heat
Bruining et al. (2009) as shown in Fig. 2.
capacity and rock heat capacity increases. However for all practical
All the calculation are given in Appendix A. We look for
values of the gas heat capacity the ratio remains smaller than one.
solutions of (9)–(11), (13) in the form of traveling waves moving
with velocity V, given by
μm vY þ μYm vZ 4. Kerogen combustion
V¼ : ð20Þ
μYm þ μm
Notice that by assuming μY 4 0, μZ 4 0, it follows that vZ o In this section it is considered that the reaction is only due to
vY ⟺vZ o V o vY . The temperature of the combustion wave θb in the immobile part of kerogen and the methane content of the
this case is given by reservoir is neglected. Thus the system (9)–(14) can be reduced to

 (9), (10), (12), (14) with zero methane oxidation rate Φm  0.
 Em We use the same method as in the previous section and analyze
expð  E m Þ  E m Ei ðE m Þ  ðθb þ 1Þ exp b

  θ þ1  the combustion process occurring in the reaction layer and the pre-
Em μ  μ V  vθ 
E m Ei b ¼ 2θ ðV  vZ Þ2 ln1 θb m : ð21Þ reaction layer. All the calculation are given in Appendix B. We look for
θ þ1 μm μθm V  vZ  solutions of (9), (10), (12), (14) in the form of traveling waves moving
The resonance conditions happen when the combustion wave with velocity V, given by
speed coincides with the thermal wave speed (Aldushin et al., μk vY
V¼ : ð23Þ
1999; Mailybaev et al., 2011). As observed above the combustion μYk þ μk
wave speed is limited by vZ oV o vY ; the only option to reach the Notice that assuming μYk 4 0, μk 4 0, yields that the combustion wave
resonance condition is to increase the speed vθ . Using (15) one seed is smaller than the oxygen wave speed V o vY . The tem-
realizes that this can only be achieved by increasing the gas density perature of the combustion wave θb and the kerogen concentration
ρg. ρck in this case are given by
cg ρg μ vY þ μYm vZ C m μm vY þμYm vZ


¼ vθ ¼ V ¼ m ⟹ρg ¼ : ð22Þ E E
Cm μYm þ μm cg μYm þ μm expð  E k Þ EEi ðE k Þ ðθb þ 1Þ exp b k þ EEi b k
θ þ1 θ þ1
μ
For the values from Table 1 the density from (22) correspond to ¼ θk2 V 2 lnj 1  ρck j ð24Þ
ρg ¼ 1:4  105 mol=m3 . Using the ideal gas law, this density corre- ðμk Þ
spond to the pressures over 3500 bar . Unfortunately, this is not and
technically feasible for all possible values of RZ A ½0; 1. It is technically
μk V  vθ
ρck ¼ θb : ð25Þ
μθk V
Solving the Eqs. (24) and (25) the temperature θb and the kerogen
concentration ρck are obtained. Notice that the dimensionless quan-
tities involved depend on the injection speed uinj, the initial concen-
tration of the immobile part of kerogen ρres k and the injected oxygen
concentration Yinj. Fig. 3 (left) plots the combustion temperature θb,
which is a three valued function of the injection speed (i.e., low,
medium and high temperature branches). Notice that the approxima-
tion used in Section 4 is only valid for high temperatures, and not for
Fig. 2. Schematic representation of in-situ combustion reaction zone separated in low and medium temperatures, because in this case the separation in
two layers. a two layered system is not accurate (Bruining et al., 2009). As will be
 G. Chapiro, J. Bruining / Journal of Petroleum Science and Engineering 127 (2015) 179–189 183

700

600

θb [K]
500

400

300
0 10 20 30 40 50
u [m/s]
inj

Fig. 3. Left figure: possible temperatures of the combustion front θb as function of the injection speed. Right figure: the high temperature branch of the combustion front as
function of injection gas speed and oxygen concentration of the injected gas. Both figures correspond to an immobile kerogen concentration of ρres
k ¼ 25 kg=m . The bold line
3

on the left figure corresponds to a cross section of the right figure at Y inj ¼ 0:21 indicated by the bold line in both figures.

Fig. 4. This figure shows the regions in the parameter space ρkres and Yinj that correspond to temperatures higher and lower than 650 1C (left) and 400 1C (right). The two curves
correspond to injection speeds uinj ¼ 20 m/day (left) and uinj ¼ 50 m/day (right). It shows the amount of oxygen to be used as a function of immobile fuel present in a reservoir.

shown in Section 5 the numerical simulations confirmed that indeed 5.1. Methane combustion
only the upper branch is stable. The combustion temperature depends
on the injection speed uinj, the initial fuel concentration ρres
k and the Fig. 5 shows the simulation results for the system (9)–(14) for
oxygen concentration of the injected gas Yinj. On the right side of Fig. 3 methane combustion (with kerogen reaction disregarded) corre-
only the high temperature branch that corresponds to sponding to the parameter values from Table C1. Notice that
ρres
k ¼ 25 kg=m is plotted.
3
methane is driven by the bank of combustion gases away from
Fig. 4 uses (24) and (25) and plots the regions in the parameter the high temperature zone and therefore combustion will stop.
space uinj  Y inj for which the combustion front attains tempera-
tures higher than 650 1C or higher than 400 1C on the left and right 5.2. Kerogen combustion
sides respectively.
Figs. 3–6 were plotted using atmospheric pressure, however, Fig. 6 shows the simulation results for the system (9)–(14) for
when high pressures are used the simulations do not change much kerogen combustion (with the methane reaction disregarded) cor-
because of the simplified model. responding to the parameter values in Table C1. Notice that the
dimensionless combustion temperature is θb ¼ 1:3, which corre-
Remark. Resonance conditions happen when the combustion sponds to a temperature of T bsim ¼ 690 1C. The approximated ana-
wave speed coincides with the thermal wave speed, (Mailybaev lytical formula gives T ban ¼ 680 1C, see Fig. 4, which corresponds to
et al., 2011). In practice, resonance conditions are achieved by a difference of less that 1.5%. We can estimate the combustion wave
changing the injection oxygen fraction, which can vary inside the speed from Fig. 4, obtaining V ksim ¼ 0:0417, which is within the
interval Y n A ½0; 1. Equating the combustion wave speed (23) with roundoff accuracy the same as V kan ¼ 0:0417 obtained from the
thermal wave speed (15) the resonance condition is obtained: analytical formula (23).

1 In Section 6.2 we discuss the influence of the temperature
kYk ρnk 1
Yn ¼  1  RY : increment on the production, as shown in Fig. 7.
kk φvθ

6. Production enhancement
5. Numerical solutions
6.1. Permeability enhancement
In this section the solutions obtained by analytical approxima-
tions in Sections 3 and 4 are validated. For this purpose the Crank– Gas shales and clays have permeability of the order of nanodarcy
Nicolson finite difference scheme with an adaptive time step is (Glorioso et al., 2012). The permeability of such shales includes modi-
used. In both simulations high initial temperatures were used to fication due to the slip effect. Much research is devoted to the correct
faster reach a steady state profile. modeling of permeability in tight reservoirs. Different equations exists
184 G. Chapiro, J. Bruining / Journal of Petroleum Science and Engineering 127 (2015) 179–189

5 5
θ θ
4 ρ 4 ρ
Y Y
3 3

2 2

1 1

0 0

0 200 400 600 0 200 400 600

ξ ξ
Fig. 5. Simulation results for methane combustion. Initial conditions are plotted on the left and solution profile at time 1870 on the right.

5 5
θ θ
4 ρ 4 ρ
Y Y
3 3

2 2

1 1

0 0

0 200 400 600 0 200 400 600

ξ ξ
Fig. 6. Simulation results for kerogen combustion. Initial conditions are plotted on the left and the solution profile at time 10,000 on the right.

100 100
P [bar]

P [bar]

50 50

0 0
10 10
1000 1000
5 5
r [m 500 r [m 500
] ]
0 x [m] 0 x [m]

Fig. 7. Pressure profile after 500 days. Left figure corresponds to k ¼100 nD and the right figure corresponds to k ¼ 500 nD.

that incorporate these effects in the Darcy flow equation, e.g., (Civan reservoir temperature
et al., 2011; Fathi et al., 2012). In order to obtain the permeability as 
 
T  T res
function of temperature k(T), the simplified approach based on Jam- kðT; tÞ ¼ max min exp lnð5Þ kini ; 5 kini ; ð26Þ
tZ0 350
aluddin et al. (2000) is used.
In Jamaluddin et al. (2000) the mechanism of the formation heat where kini is the initial permeability.
treatment is described. This technique is based on the vaporization of Other mechanisms like thermally induced fracking and depen-
the blocked water, dehydration of the clay-bound water, destruction dence on the geological origin of shales must be kept for
of the clay lattices and creation of microfractures due to induced future work.
thermal stresses.
Jamaluddin et al. (2000) indicate that the permeability increased
by approximately 50% to 675% depending on the temperature and 6.2. Recovery rate computation
time of heating. Maximum values were obtained at a temperature of
650 1C and heating times of around one hour. In Jamaluddin et al. We consider that characteristic time of the changes in the pressure
(2000) different types of clays were analyzed and here it is assumed field (Dake, 1983, Chapter 5) in response to the permeability modifica-
that similar behavior occurs for the gas shales discussed. For simplicity tion is considerably longer than the characteristic time of the change in
it is considered that the permeability kðT; tÞ grows exponentially by a the permeability field which is the typical time of the combustion
maximum factor of 5 when the temperature grows to 350 1C over the process. This allows us to solve the equations sequentially. After the
 G. Chapiro, J. Bruining / Journal of Petroleum Science and Engineering 127 (2015) 179–189 185

combustion process is performed the final permeability field is The possibility of using this method for recovery of shale gas needs
obtained. to consider a number of other mechanisms such as gas production
In order to illustrate the effect of the permeability enhancement from kerogen, the effect of the initial permeability of the reservoir,
on the gas recovery efficiency we perform two simulations with compression costs and coupling it to existing techniques, which are
constant permeabilities: the first one uses the original reservoir not addressed in this work.
permeability given in Table C1; the second one corresponds to the
situation when the combustion wave already passed through the
reservoir enhancing the permeability to the maximum value 7. Discussion and conclusions
described by Eq. (26). The computation uses a cylindrical reservoir
with length Lres ¼ 1000 m, radius Rres ¼ 10 m, initially filled with  We investigated in-situ combustion for improving the perme-
methane at a pressure P res ¼ 107 Pa. In the middle of the reservoir ability in a gas shale reservoir, which is possibly a supplemen-
there is a hole of small size representing a well that has a constant tary technique for hydraulic fracturing.
pressure of P hole ¼ 106 Pa. This allows us to calculate the methane  A one dimensional combustion model was considered. Under
production profile as explained below. simplifying hypotheses, a quasianalytical solution of the combus-
The temperature dependence on the pressure P is neglected and tion traveling wave was obtained. We consider two possibilities for
it is assumed that the pressure inside the formation is described by the in-situ fuel composition, namely, methane or kerogen.
the diffusion equation:  Methane combustion did not produce a high temperature tra-
! veling wave. Thus it cannot be used for heating the reservoir.
∂P
φct ρg ¼ ∇ðr xÞ 
kðT; tÞ
ρg ∇ðr xÞ P ; ð27Þ  However to obtain a high temperature combustion wave resu-
∂t μg lting from the oxidation of kerogen is possible. The solution was
analyzed with parameters that correspond to resonance condi-
where the constants are given in Table 1, kðT; tÞ is the permeability
tions. We present the set of parameters for which the combustion
of the porous medium, which depends on the temperature history.
reaches the optimal temperature to improve the permeability.
A gas density ρg is described by the ideal gas law at the reference
 In the case that reservoir is initially filled with a mixture of kerogen
temperature, (i.e., ρg ¼ M g P =ðRT n Þ) yielding
! and methane, it can be expected that the combustion will be
n
∂P kðT ; tÞ M g P described only by kerogen oxidation as methane will be displaced
φct ρg ¼ ∇ðr xÞ  ∇ ðr xÞ P : ð28Þ
∂t μg RT n by the combustion products out of reach of the injected oxygen.
 Our theoretical calculations were compared with numerical
In cylindrical coordinates (r,x) with radial symmetry, (28) becomes simulations using finite difference schemes. Both results show
∂P kðT ; tÞM g 1 good agreement with each other.
φct ρg ¼ ∂r ðr P ∂r P Þ þ ∂x ðP ∂x P Þ : ð29Þ  Finally, a simple cylindrical reservoir model with radial sym-
∂t μg RT n r
metry was used to illustrate numerically the pressure changes
Using (7), Eq. (29) is rewritten in dimensionless form: in the reservoir and ensuing production improvement resulting
∂P 1 from permeability enhancement. As explained in the Appendix
¼ ∂r ðer r∂r PÞ þ ∂x ðex ∂x PÞ; ð30Þ the recovery efficiency improves from 1.19% to 4.25% of the gas
∂t r
initially in place for a production of 500 days and a domain
where R  L ¼10 m  1000 m.
Mg P n tn e kðT; tÞ e kðT; tÞ  In 2-D or 3-D the combustion process may or may not be stable.
e¼ ; er ðTÞ ¼ ; ex ðTÞ ¼ : ð31Þ
μg RT n  φct ρg ðr n Þ2 ðxn Þ2 This can be analyzed using a normal mode analysis (Armento
and Miller, 1977) or probabilistic simulation (Biezen et al., 1995;
It is assumed that an initial pressure inside the formation is Pres King and Scher, 1985; King, 1987). Such an analysis, however,
and the pressure inside the hole Phole is maintained constant at all must be kept for future work.
times. Thus the problem can be solved in a cylindrical domain  The possibility of using this method for recovery of shale gas
½1; Rres   ½0; Lres  with boundary conditions: needs to consider a number of mechanisms such as gas pro-
∂P duction from kerogen, initial permeability of the reservoir,
Pðx; r ¼ 1; tÞ ¼ P hole ; ðx; r ¼ Rres ; tÞ ¼ 0;
∂x compression costs and coupling it to existing techniques. This
∂P ∂P has to be kept for future work.
ðx ¼ 0; r; tÞ ¼ 0; ðx ¼ Lres ; r; tÞ ¼ 0: ð32Þ
∂r ∂r  The aim of the present paper is to show that it is possible to
The initial condition is sustain a temperature wave that can reach 650 1C to increase the
permeability according to Jamaluddin et al. (2000) without heating
Pðx; r; t ¼ 0Þ ¼ P res : ð33Þ the full reservoir. Even if this is an important issue, for the full
The initial total volume of a reservoir is given by feasibility and validity of our model it is necessary to experimen-
tally prove that permeability enhancement also occurs when the
V res ¼ φπR2res Lres : ð34Þ formation is heated by combustion. Another issue consists in
Thus, at each time step, the methane content inside the cylindrical studying the effect of the kerogen combustion on the permeability
domain is calculated by integrating the pressure and applying ideal enhancement. However, distillation and condensation can lead to
gas law deterioration of the permeability. We are aware that this aspects
Z Z require further experimental and theoretical studies.
PV 2π Rres Lres
n ¼ rPðr; xÞ dx dr: ð35Þ
RT RT 0 0

Fig. 7 plots the pressure profiles at 500 days. The cumulative

production (the difference between the initial methane content Acknowledgments
and the methane content at 500 days) for a permeability of 0:1 μD
(left figure) is 1.19% of total available gas whereas for a perme- G.C. was supported in part by Coordenação de Aperfeiçoamento de
ability of 0:5 μD (right figure) it is 4.25%; the resulting recovery Pessoal de Nível Superior (CAPES) under Grants 17710-12-0 and 024/
increment is 350%. 2011.
186 G. Chapiro, J. Bruining / Journal of Petroleum Science and Engineering 127 (2015) 179–189

The authors thank Prof. Dr. A. Mailybaev, Prof. Dr. D. Marchesin transfer is dominated by both conduction and flow. Thus the
and Prof. Dr. K.H. Wolf for helpful discussions. The authors thank combustion rate (13) can be neglected in the energy balance and
TU Delft for its hospitality. The authors thank W.S. Pereira for can be approximated with θ¼0 for the other equations, i.e.,
helping to find some typos and an anonymous referee for valuable
Φapprox
m ¼ YZ expð  E m Þ: ðA:11Þ
comments and suggestions.
Using traveling wave coordinates ξ ¼ x  Vt the system (9)–(11)
becomes
Appendix A. Methane combustion
ðvθ V Þdξ θ ¼ dξξ θ; ðA:12Þ
In this section we present the analytical calculations that allow
us obtain the combustion wave speed (20) and the combustion ðvY  V Þdξ Y ¼  μYm Φapprox
m ; ðA:13Þ
wave temperature (21) corresponding to the system (9)–(11), (13).
In traveling coordinates ðξ ¼ x  VtÞ the system (9)–(11), becomes ðvZ  VÞdξ Z ¼  μm Φapprox
m : ðA:14Þ
ðvθ  VÞdξ θ ¼ dξξ θ þ μθm Φm ; ðA:1Þ c
Substituting (A.14) into (A.13) and integrating from ξ to ξ and next
substituting the boundary condition at ξ ¼ ξ0 , i.e., Y 0 ¼ 0 and
ðvY  VÞdξ Y ¼  μYm Φm ; ðA:2Þ Z 0 ¼ 1, it follows that

ðvZ  VÞdξ Z ¼  μm Φm : ðA:3Þ vY  V vZ V vY  V c vZ  V c vY  V 0 vZ V 0

Y Z¼ Y  Z ¼ Y  Z :
μYm μm μYm μm μYm μm
We simplify Eqs. (A.1)–(A.3) for the reaction layer and pre-
ðA:15Þ
reaction layer in the next subsections.
0
Substituting values at ξ , it follows that
A.1. Reaction layer
Y ¼ αðZ  1Þ; Y c ¼ αðZ c  1Þ; ðA:16Þ

The reaction layer is the first layer from upstream and corre- where α is defined in (A.10). Substituting (A.16) into (A.14) and
sponds to ξb rξ rξc . We assume that this layer is thin and thus using (A.11) leads to
that the heat propagation is dominated by conduction. We neglect μY
the first derivative terms in (A.1). Consequently, it is plausible that dξ Z ¼ Zð1  ZÞ expð  E Þ ¼ KZð1  ZÞ; ðA:17Þ
vY  V
the temperature variation inside this layer is small, yielding θb 
θ  θc ; in other words the temperature variation is zero, i.e., where K ¼ μYm =vY  V expð  E Þ: Eq. (A.17) can be solved, i.e.,
dξ θb ¼ 0, dξ Y b ¼ dξ Z b ¼ 0. Thus the oxygen and methane partial Z Z Z ξ
dZ K c eKðξ  ξc Þ
pressures at ξ ¼ ξb are equal to the injection conditions Yb ¼1, ¼ Kdξ⟹Z ¼ ; ðA:18Þ
c Zð1  ZÞ 1 þ K c eKðξ  ξc Þ
Z ξc
Zb ¼0. Combining all of the above Eq. (A.1) becomes
0 ¼ dξξ θ þ μθm Φm ; ðA:4Þ where K c ¼ Z c =ð1 Z c Þ. We solve (A.12) integrating from ξc to ξ and
next substituting the boundary conditions at ξ ¼ ξ0 yielding
Substituting (A.3) into (A.4), integrating from ξb to ξ and next
substituting ξc leads to dξ θ þ ðV  vθ Þθ ¼ dξ θc þ ðV  vθ Þθc ¼ 0: ðA:19Þ

μm dξ θ μθm ðvZ  VÞZ ¼ μm dξ θb  μθm ðvZ  VÞZ b Solving (A.19) the wave profile in the pre-reaction zone is obtained
¼ μm dξ θc  μθm ðvZ V ÞZ c : ðA:5Þ θ ¼ θc exp ðvθ  VÞðξ ξc Þ ; dξ θc ¼ θc ðvθ  VÞ: ðA:20Þ
b
Substituting (A.3) into (A.2), integrating from ξ to ξ and next
substituting ξc leads to
A.3. Connecting the reaction and pre-reaction layers
vY  V vZ  V vY  V b vZ  V b vY  V c vZ  V c
Y Z¼ Y  Z ¼ Y  Z : ðA:6Þ
μYm μm μYm μm μYm μm Basically, the wave speed has to be the same inside both layers
A wave profile in the reaction layer is described by and the heat flow is continuous, i.e., the derivative of θ is equal at
μ both sides of the layer boundary. Thus the relations for the
dξ θ ¼ θm ðvz V ÞZ; ðA:7Þ combustion speed, the combustion temperature and the methane
μm
partial pressure at point ξ ¼ ξc are obtained.
Y ¼ 1 þαZ; ðA:8Þ Wave speed: Equating Yc from (A.10) and (A.16) it follows that
α ¼  1 and that the combustion wave speed is given by Eq. (20).
μm
dξ Z ¼ Φm ; ðA:9Þ Heat conservation: Equating the term dξ θc from (A.10) and (A.20)
V  vZ yields
where at ξ ¼ ξb the boundary conditions Z b ¼ ∂ξ θb ¼ 0, Yb ¼ 1 and μθm
Zb ¼0 were used. For ξ ¼ ξc it follows: ðvZ V ÞZ c ¼ θc ðvθ V Þ: ðA:21Þ
μm
μθm μYm V  vZ
dξ θ c ¼ ðvz  VÞZ c ; Y c ¼ 1 þ αZ c ; where α¼ : Considering θb ¼ θc it follows:
μm μm V  vY
ðA:10Þ V  vZ μθm c μm V  vθ
θb ¼ Z or Z c ¼ θb : ðA:22Þ
V  vθ μm μθm V  vZ
The solution of the system (A.7)–(A.9) is tedious and involves
continuity conditions, which appear in Section A.3. We solve this Combustion temperature inside reaction layer: Dividing (A.7) by
system in Section A.3. (A.9) and substituting (13), (A.8) results in the following ordinary
differential equation (ODE):
A.2. Pre-reaction layer

dθ μ 1 Em
¼  θm ðV  v Z Þ 2
exp ; ðA:23Þ
dZ μ2m 1Z θþ1
The pre-reaction layer corresponds to ξc r ξ o 1, see Fig. 2. We
assume that inside this layer the reaction is slow, thus the heat Separating the variables and integrating inside the reaction layer
 G. Chapiro, J. Bruining / Journal of Petroleum Science and Engineering 127 (2015) 179–189 187

yields  Vdξ ρk ¼  μk Φk : ðB:3Þ

Z θc
 Z Zc
 Em μ dZ
exp dθ ¼  θm ðV  vZ Þ 2
: ðA:24Þ
θb θþ1 μ2m Zb 1  Z B.1. Reaction layer
First the left side of (A.24) is evaluated
The reaction layer is the first layer from upstream and corre-
Z θc
 Z 0

E m  Em sponds to ξb r ξ r ξc . Under the same hypotheses of Section A.1 the
exp dθ  exp dθ
θb θþ1 θb θþ1 boundary conditions are dξ θb ¼ 0, dξ Y b ¼ dξ ρbk ¼ 0, Yb ¼1, ρbk ¼ 0 and


 0 the system becomes
 Em Em
¼ exp ðθ þ 1Þ  E m Ei ; ðA:25Þ
θþ1 θþ1 θb 0 ¼ dξξ θ þ μθk Φk ; ðB:4Þ

where we use the Exponential Integral Function Ei(x), see ðvY  V Þdξ Y ¼  μYk Φk ; ðB:5Þ
Abramowitz and Stegun (1964). In Eq. (A.25) we neglect the
following integral with respect to other terms:  Vdξ ρk ¼  μk Φk : ðB:6Þ
Z θc

 Em After analogous calculations as in Section A.1 the wave profile in
exp dθ: ðA:26Þ
0 θþ1 the reaction zone is described by equations corresponding to (A.7)–
This approach is valid for high activation energies, see Bruining (A.10):
et al. (2009). To obtain simpler expressions it is useful to follow  μθk V
dξ θ ¼ ρk ; ðB:7Þ
Bruining et al. (2009) and expand the fraction inside the exponen- μk
tial around θb:
Y ¼ 1  βρk ; ðB:8Þ
E m  Em  Em E m δθ
¼  þ : ðA:27Þ
θ þ 1 θb þ 1 þδθ θb þ 1 ðθb þ 1Þ2 μk
dξ ρk ¼ Φ : ðB:9Þ
V k
In this case one can approximate the integral in Eq. (A.25) by
Z 0
 Z 1 ! and
 Em E m E m δθ
exp dθ  exp b þ b dðδθÞ  μθk V c μYk V
θb θþ1 0 θ þ 1 ðθ þ 1Þ2 dξ θc ¼ ρk ; Y c ¼ 1  βρck ; where β¼ : ðB:10Þ
" !#  1 μk μk vY  V


 E m ðθb þ 1Þ2 E m δθ
¼ exp b exp Notice that for kerogen the quantity corresponding to the speed vZ
θ þ1 Em ðθb þ 1Þ2 0 is zero.

 b
 E m ðθ þ 1Þ2
¼  exp b : ðA:28Þ
θ þ1 Em B.2. Pre-reaction layer

In the integral calculation, as before, we use the fact that the strong
As before, the pre-reaction layer corresponds to ξc r ξ o 1.
nonlinearity in the exponential function allows us to change the
Under the same hypotheses as in Section A.2, the combustion rate
integral limits from ½0;  θb  to ½0; 1 introducing a small error.
(14) can be approximated for the other equations by
This argument is similar to the one used for the integration of Eq.
(A.25). Now the right side of (A.24) is evaluated. Φapprox
k
¼ ρk Y expð  E k Þ:
Z Zc Z Zc
dZ dZ c Using traveling wave coordinates the system (9), (10), (12)
 ¼  lnj 1  Z j Z0 ¼  lnj 1  Z c j : ðA:29Þ
Z b 1 Z 0 1Z
becomes

Substituting (A.25) and (A.29) back into (A.24) a compete relation ðvθ V Þdξ θ ¼ dξξ θ; ðB:11Þ
between Zc and θb is obtained

 ðvY  V Þdξ Y ¼  μYk Φapprox ; ðB:12Þ
 Em k
expð  E m Þ  E m Ei ðE m Þ  ðθb þ 1Þ exp b
θ þ1  Vdξ ρk ¼  μk Φapprox : ðB:13Þ

 k
Em μθ
þE m Ei b ¼ 2 ðV  vZ Þ lnj 1 Z c j :
2
ðA:30Þ After similar calculations as in Section A.2 a wave profile in the pre-
θ þ1 μm
reaction zone is described by
To obtain a simpler expression one substitutes (A.28), (A.29) into
(A.24) leading to the approximate formula: Y ¼ βð1  ρk Þ; Y c ¼ βð1  ρck Þ; ðB:14Þ

 where β is defined in (B.10) and
 E m ðθb þ 1Þ2 μθm
exp b ¼ 2 ðV  vZ Þ2 lnj 1 Z c j : ðA:31Þ
θ þ1 Em μm K c eKðξ  ξc Þ
ρk ¼ ; θ ¼ θc exp ðvθ  VÞðξ  ξc Þ ; dξ θc ¼ θc ðvθ  VÞ:
1 þ K c eKðξ  ξc Þ
ðB:15Þ
Appendix B. Kerogen combustion

B.3. Connecting the reaction and pre-reaction layers

In this section we present the analytical calculations that allow
us obtain the combustion wave speed, temperature and fuel
Like in Section A.3 the relations for the combustion speed, the
concentration (23), (24) and (25) corresponding to the system
combustion temperature and the methane partial pressure at point
(9), (10), (12) and (14). In traveling coordinates the system (9), (10),
ξ ¼ ξc are obtained in a next section.
(12), becomes
Wave speed: Comparing (B.10) and (B.14) it follows that β¼1 and
ðvθ  VÞdξ θ ¼ dξξ θ þ μθk Φk ; ðB:1Þ the combustion wave speed is given by (23).
Heat conservation: Following Section A.3 and considering θb ¼ θc
ðvY  VÞdξ Y ¼  μYk Φk ; ðB:2Þ results in a relation between the combustion temperature and the
188 G. Chapiro, J. Bruining / Journal of Petroleum Science and Engineering 127 (2015) 179–189

amount of kerogen present in the interface: Reaction rate parameters: The activation energy Ek and pre-
exponential coefficient kpk for kerogen were taken the same as for
V μθk c μ V  vθ
θb ¼ ρ or ρck ¼ θb k : ðB:16Þ coke combustion, see Chapiro et al. (2012). The activation energy
V  vθ μk k μθk V
Em and pre-exponential coefficient kpm for methane combustion
Combustion temperature inside reaction layer: Following Section inside porous media were taken from Bubnovich et al. (2006) and
A.3 a complete formula is obtained Yoshio et al. (1988).


 Porosity: The porosity changes if one assumes kerogen combus-
E E
expð  E k Þ  EEi ðE k Þ  ðθb þ 1Þ exp b k þ EEi b k tion. However, there are other effects that may influence the
θ þ1 θ þ1
μθk 2 porosity such as the evaporation of trapped water, geochemical
c
¼ V lnj 1  ρk j ; ðB:17Þ reactions, etc. To obtain an accurate porosity model requires a
ðμk Þ2
separate study which is outside the scope of this paper. To keep the
and also an approximated formula: text tractable we presently disregard all these effects.

 b
E ðθ þ 1Þ2 μ
exp b k ¼ θk2 V 2 lnj 1  ρck j : ðB:18Þ References
θ þ1 Ek ðμk Þ

Abramowitz, M., Stegun, I., 1964. Handbook of Mathematical Functions with

Formulas, Graphs, and Mathematical Tables. Dover Publications.
Appendix C. Model parameters Abu-Khamsin, S., Brigham, E., Ramey Jr., H., 1988. Reaction kinetics of fuel formation
for in-situ combustion. SPE Reserv. Eng. 3 (4), 1308–1316.
Akkutlu, I.Y., Fathi, E., 2011. Gas transport in shales with local kerogen heterogene-
Average heat capacity of gas: As the gas heat capacity is small ities. In: SPE Annual Technical Conference and Exhibition, SPE.
compared to the rock heat capacity, its dependence on composition is Akkutlu, I.Y., Fathi, E., 2012. Multiscale gas transport in shales with local kerogen
heterogeneities. SPE J. 17 (04), 1002–1011.
neglected and considered to be equal to the methane heat capacity. Aldushin, A., Rumanov, I., Matkowsky, B., 1999. Maximal energy accumulation in a
Consider the heat capacity of methane cg ¼ 35.69 J/mole K and the superadiabatic filtration combustion wave. J. Combust. Flame 118, 76–90.
average gas density ρg ¼ 45 mole=m3 , thus C g ¼ cg ρg  1580 J=m3 K. Allan, A.M., Mavko, G., 2012. The effect of adsorption and diffusion on the gas
permeability of kerogen. In: 2012 SEG Annual Meeting, SEG.
Enthalpy of methane: The methane combustion molar enthalpy Armento, M.E., Miller, C.A., 1977. Stability of moving combustion fronts in porous
at standard temperature is U m ¼ 8:9  105 J/mole. Thus the amount media. Soc. Pet. Eng. J. 17 (06), 423–430.
of energy per cubic meter of porous media is Q m ¼ ρg φU m ¼ Biezen, E., Bruining, J., Molenaar, J., 1995. An integrated 3d model for underground
coal gasification. In: Society of Petroleum Engineers. Annual Technical Con-
45 mole=m3  0:2ðÞ 8:9  105 J=mole ¼ 8:0  106 J/m3.
ference, pp. 929–940.
Enthalpy of kerogen: The coke combustion molar enthalpy is Bird, R., Stewart, W., Lightfoot, E., 2002. Transport Phenomena, 2nd Edition Wiley,
taken equal to carbon molar enthalpy, i.e., U k ¼ 4:0  105 J/mole. New York.
Therefore the amount of energy per cubic meter of porous media is Birol, F., 2011. World Energy Outlook. International Energy Agency, Paris, France.
Bruining, J., Mailybaev, A., Marchesin, D., 2009. Filtration combustion in wet porous
Q k ¼ ρk ð1  φÞU k ¼ 372 mole=m3  0:8 ½  4:0  105 J=mole medium. SIAM J. Appl. Math. 70, 1157–1177.
¼ 1:2  108 J/m3. Bubnovich, V., Zhdanok, S., Dobrego, K., 2006. Analytical study of the combustion
Stoichiometric coefficients: We consider the following chemical waves propagation under filtration of methane-air mixture in a packed bed. Int.
J. Heat Mass Transf. 49 (15), 2578–2586.
reactions for methane CH4 þ2O2 -CO2 þ2H2 O. The residue of kero- Busch, A., Amann-Hildenbrand, A., 2013. Predicting capillarity of mudrocks. Mar.
gen ðCH1:5 O0:05 N0:023 S0:005 Þ, see Yen and Chilingarian (1976), after Petrol. Geol. 45, 208–223.
pyrolysis and cracking, is considered to be pure carbon (C). We Busch, A., Alles, S., Gensterblum, Y., Prinz, D., Dewhurst, D., Raven, M., Stanjek, H.,
Krooss, B., 2008. Carbon dioxide storage potential of shales. Int. J. Greenh. Gas
consider immobile part of kerogen to behave like coke and react as Control 2 (3), 297–308.
C þ O2 -CO2 . Thus the stoichiometric coefficients are kYm ¼ 2, kYk ¼ 1, Chapiro, G., Mailybaev, A.A., Souza, A., Marchesin, D., Bruining, J., 2012. Asymptotic
km ¼ 1 and kk ¼ 1. approximation of long-time solution for low-temperature filtration combustion.
Comput. Geosci. 16, 799–808.
Chelini, V., Muttoni, A., Mele, M., Rossi, E., Galimberti, R., Ortenzi, A., 2010. Gas-shale
Table C1 reservoir characterization: a North Africa case. In: SPE Annual Technical
Dimensional parameters and their typical values. Conference and Exhibition.
Cinar, M., Castanier, L.M., Kovscek, A.R., 2009. Isoconversional kinetic analysis of the
Symbol Physical quantity Value Unit combustion of heavy hydrocarbons. Energy Fuels 23 (8), 4003–4015.
Civan, F., Rai, C., Sondergeld, C., 2011. Shale-gas permeability and diffusivity inferred
φ Porosity 0.2 () by improved formulation of relevant retention and transport mechanisms.
c m ρm Heat capacity of the reservoir J/m3K Transp. Porous Media 86 (3), 925–944.
2  106
Dake, L., 1983. Fundamentals of Reservoir Engineering. Elsevier.
cg ρg Average heat capacity of gas 1580 J/m3K
Essenhigh, H., 1981. Fundamentals of coal combustion. In: Chemistry of Coal
ug Darcy velocity of total gas (1 m/day) 5 m/s
1:0  10 Utilization. Secondary supplementary volume. Wiley-Interscience Publication,
RY Sorbtion coef. of oxygen 0 () New York.
RZ Sorbtion coef. of methane ðRZ A ½0; 1Þ 1 () Fathi, E., Akkutlu, I.Y., 2012. Lattice Boltzmann method for simulation of shale gas
D Thermal conductivity of porous medium 0.87 J/(m s K) transport in kerogen. SPE J. 18, 27–37.
Qm Heat content of fuel (methane) 8  106 J/m3 Fathi, E., Tinni, A., Akkutlu, I., 2012. Shale gas correction to Klinkenberg slip theory.
Qk Heat content of fuel (kerogen) 1.2  108 J/m3 In: SPE Americas Unconventional Resources Conference.
R Gas constant 8.314 J/mol K Gates, C.F., Ramey Jr., H.J., 1958. Field results of South Belridge thermal recovery
Tres Initial reservoir temperature 300 K experiment. AIME Pet. Trans. 213, 236–244 SPE 1179-G.
Glorioso, J., Rattia, A., et al., 2012. Unconventional reservoirs: basic petrophysical
P0 Reference pressure (1 atm) 105 Pa
concepts for shale gas. In: SPE/EAGE European Unconventional Resources
P Considered pressure (1 atm) 105 Pa
Conference & Exhibition-From Potential to Production.
kpm Pre-exponential coef. of methane 2.6  108 1/s Hagoort, J., 1981. Waterflood-Induced Hydraulic Fracturing (Ph.D. thesis), Tech-
kpk Pre-exponential coef. of kerogen 500 1/s nische Hogeschool Delft.
kYm Stoichiometric coef. (O2 with methane) 2 () Howarth, R., Santoro, R., Ingraffea, A., 2011. Methane and the greenhouse-gas
kYk Stoichiometric coef. (O2 with kerogen) 1 () footprint of natural gas from shale formations. Clim. Change 106 (4), 679–690.
km Stoichiometric coef. (methane) 1 () Jamaluddin, A., Bennion, D., Thomas, F., Ma, T., 2000. Application of heat treatment
kk Stoichiometric coef. (kerogen) 1 () to enhance permeability in tight gas reservoirs. J. Can. Pet. Technol. 39 (11),
Em Activation energy of methane 13,0000 J/mol 1–16.
Ek Activation energy of kerogen 58,000 J/mol Kennedy, R., Knecht, W., Georgi, D., 2012. Comparisons and contrasts of shale gas
Yinj Injected oxygen concentration 0.21 mol/mol and tight gas developments, North American experience and trends. In: SPE
Zres Reservoir methane concentration 1.0 mol/mol Saudi Arabia Section Technical Symposium and Exhibition.
ρres Initial immobile kerogen ρres kg/m3 King, M., 1987. Viscous fingering utilizing probabilistic simulation. In: SPE Annual
k k A ½1; 50 25
Technical Conference and Exhibition. Society of Petroleum Engineers.
 G. Chapiro, J. Bruining / Journal of Petroleum Science and Engineering 127 (2015) 179–189 189

King, M., Scher, H., 1985. Probabilistic stability analysis of multiphase flow in porous Seewald, J., 2003. Organic–inorganic interactions in petroleum-producing sedimen-
media. In: SPE Annual Technical Conference and Exhibition. Society of Petro- tary basins. Nature 426 (6964), 327–333.
leum Engineers. Smith, I., 1978. The intrinsic reactivity of carbons to oxygen. Fuel 57 (7), 409–414.
Kuuskraa, V., Stevens, S., Van Leeuwen, T., Moodhe, K., 2011. World Shale Gas Soeder, D., 1988. Porosity and permeability of eastern devonian gas shale. SPE Form.
Resources: An Initial Assessment of 14 Regions Outside the United States. US Eval. 3 (01), 116–124.
Department of Energy. Yen, T., Chilingarian, G., 1976. Oil Shale. Elsevier Scientific Publishing Co.,
Landau, L., Lifshitz, E., 1987. Fluid Mechanics, vol. 6 (Course Of Theoretical Physics). Amsterdam.
Pergamon Press Oxford. Yoshio, Y., Kiyoshi, S., Ryozo, E., 1988. Analytical study of the structure of radiation
Levenspiel, O., 1999. Chemical Reaction Engineering, 3rd Edition Wiley, New York. controlled flame. Int. J. Heat Mass Transfer 31 (2), 311–319.
Mailybaev, A., Marchesin, D., Bruining, J., 2011. Resonance in low-temperature Zhang, X., Changan, D., Deimbacher, F., Crick, M., Harikesavanallur, A., 2009.
oxidation waves for porous media. SIAM J. Math. Anal. 43, 2230–2252. Sensitivity studies of horizontal wells with hydraulic fractures in shale gas
Perry, K., Lee, J., 2007. Report: hard truths-facing the hard truths about energy. Topic reservoirs. In: International Petroleum Technology Conference.
Paper No. 27: Unconventional Gas Reservoirs—Tight Gas, Coal Seams, and
Shales. National Petroleum Council.
Prats, M., 1961. Effect of vertical fractures on reservoir behavior-incompressible fluid
case. SPE J. 1 (2), 105–118.
Sanmiguel, J., Mallory, D., Mehta, S., Moore, R., 2002. Formation heat treatment
process by combustion of gases around the wellbore. J. Can. Pet. Technol. 41 (8),
70–76.