Thermal conductivity of granular porous media: A pore scale modeling approach

R. Askari, S. Taheri, and S. H. Hejazi

Citation: AIP Advances 5, 097106 (2015); doi: 10.1063/1.4930258

View online: http://dx.doi.org/10.1063/1.4930258
View Table of Contents: http://aip.scitation.org/toc/adv/5/9
Published by the American Institute of Physics
 AIP ADVANCES 5, 097106 (2015)

Thermal conductivity of granular porous media: A pore

scale modeling approach
R. Askari,1 S. Taheri,2 and S. H. Hejazi2
1
Geological and Mining Engineering and Sciences, Michigan Technological University,
Houghton, MI 49931, United States
2
Chemical and Petroleum Engineering Department, University of Calgary,
Calgary, T2N 1N4, Canada
(Received 5 June 2015; accepted 24 August 2015; published online 2 September 2015)

Pore scale modeling method has been widely used in the petrophysical studies to
estimate macroscopic properties (e.g. porosity, permeability, and electrical resistivity)
of porous media with respect to their micro structures. Although there is a sumptuous
literature about the application of the method to study flow in porous media, there
are fewer studies regarding its application to thermal conduction characterization,
and the estimation of effective thermal conductivity, which is a salient parameter in
many engineering surveys (e.g. geothermal resources and heavy oil recovery). By
considering thermal contact resistance, we demonstrate the robustness of the method
for predicting the effective thermal conductivity. According to our results obtained
from Utah oil sand samples simulations, the simulation of thermal contact resistance
is pivotal to grant reliable estimates of effective thermal conductivity. Our estimated
effective thermal conductivities exhibit a better compatibility with the experimental
data in companion with some famous experimental and analytical equations for
the calculation of the effective thermal conductivity. In addition, we reconstruct a
porous medium for an Alberta oil sand sample. By increasing roughness, we observe
the effect of thermal contact resistance in the decrease of the effective thermal
conductivity. However, the roughness effect becomes more noticeable when there
is a higher thermal conductivity of solid to fluid ratio. Moreover, by considering
the thermal resistance in porous media with different grains sizes, we find that the
effective thermal conductivity augments with increased grain size. Our observation is
in a reasonable accordance with experimental results. This demonstrates the useful-
ness of our modeling approach for further computational studies of heat transfer
in porous media. C 2015 Author(s). All article content, except where otherwise
noted, is licensed under a Creative Commons Attribution 3.0 Unported License.
[http://dx.doi.org/10.1063/1.4930258]

I. INTRODUCTION
The effective thermal conductivity is an important factor to accurately estimate oil production
rate for the thermal recovery methods such as Steam Assisted Gravity Drainage (SAGD) process.1
The effective thermal conductivity, or other petorphysical parameters such as permeability, can
be obtained from analytical, experimental and numerical methods. Analytical approaches provide
exact solutions to a temperature field from which the effective thermal conductivity for a specific
physical model is determined (e.g., Refs. 2–5). Analytical methods have the disadvantage of over-
simplified assumptions that restrain their applications for real complicated problems.6 Based on
laboratory measurements, experimental methods (e.g. needle probe,7 parallel plates,8 and optical
scanning9) are used to determine the effective thermal conductivity.10–16 However, logistic or eco-
nomic constraints might confine the practicability of experimental measurements.17 In quest of
addressing the challenges and limitations associated with the analytical and experimental methods,
computational techniques have been widely considered in recent years (e.g. Refs. 18–27) among

2158-3226/2015/5(9)/097106/16 5, 097106-1 © Author(s) 2015

097106-2 Askari, Taheri, and Hejazi AIP Advances 5, 097106 (2015)

which the pore network modeming28 and the pore scale modeling29,19 approaches have demon-
strated promises to predict petrophysical parameters of porous media.
Macroscopic properties of porous media (e.g. porosity, electrical conductivity, permeability,
and thermal conductivity) can be obtained from their microstructures (the geometry and topology of
void spaces and solid volumes). Based on this presumption, Fatt28 introduces the network modeling
technique where the void spaces of a porous medium are represented by a graph of connected
sites. Sites are considered as pores, and are connected by several bonds that are defined as throats.
Network modeling technique qualitatively, within an order of magnitude, predicts some petrophysi-
cal parameters such as permeability for sand packs and sandstones.30–34 However, its application for
the quantitative predictions of petrophysical properties has been challenging because this method is
based on some assumptions that oversimplify a real structure of void spaces. For instance, at each
site, a particular value is assigned to the number of bonds; or throats are assumed to have the same
length.17
Bryant et al.17 introduce the pore scale modeling approach to predict macroscopic petrophys-
ical parameters of porous media by considering real shapes of grains, and other geometrical and
topological features. Based on this procedure, a pixel based pore model is reconstructed by either
a direct imaging techniques such as X-ray microcomputed tomography (micro-CT),35–37 thin sec-
tion analysis38,39 or simulation of random grains packing process.40,41 Using the reconstructed pore
models, different petrophysical parameters are determined. The advantage of pore scale modeling
over pore network modeling lies in the fact that the physics and image of a porous medium is
used directly in simulation, therefore, there is no need for any simplification of pores, throats, and
grain shapes. Furthermore, using computational simulation for grain packing process, we can create
numerous realizations of porous media for a given grain size distribution to provide a reasonable
correlation among different macroscopic parameters (for example, porosity versus effective thermal
conductivity). For these reasons the pore scale modeling approach is greatly utilized in different
petophysical studies such as the predictions of relative permeability, electrical conductivity, and
capillary pressure.42
In this study, we develop a computational algorithm to reconstruct 2D granular porous media
by which the thermal conduction behavior and the effective thermal conductivity can be studied
in pore scales. Our approach creates a granular porous medium for a given grain size distribu-
tion such as Normal, bi-Normal, Log-Normal, Weibull, or any distribution directly derived from
experimental data, in which the thermal contact resistance is simulated by applying roughness to
the grains. The approach allows simulation of random grain packing process by a swelling circle
procedure.43 In this procedure, a set of random points in a 2D plane is created where each point
is a center of a circle. Circles grow (swell) incrementally based on their given sizes. In each step
of swelling, the circles’ centers are slightly rearranged to avoid circles interferences. In addition,
in our methodology, we can increase the thermal contact resistance by the increment of grains’
roughness.
We simulate granular porous media with respect to a real grain size distribution of Utah uncon-
solidated oil sand samples at different porosities and demonstrate the potential of our method to
describe the thermal conduction process. In addition, we reconstruct a 2D granular porous media
based on a real grain size distribution of Alberta unconsolidated oil sand sample and measure its
effective thermal conductivity with respect to different roughness values. Furthermore, we estimate
the effective thermal conductivity with respect to grain size.
Grain size is a textural factor that affects the thermal conductivity.16 Rzhevsky and Novik44
observe that the thermal conductivity increases with augmented grain size. Following Rzhevsky and
Novik,44 other experimental studies confirm this observation (e.g. Refs. 16, 45, and 46). In spite
of these experimental observations, recent modeling studies report an inverse relationship between
the grain size and effective thermal conductivity.47 In this study, we successfully exhibit the grain
size effect by considering the thermal contact resistance for granular porous media. To best of the
authors’ knowledge, it is for the first time that the surface contact roughness is incorporated for a
proper pore level modeling of heat conduction through granular porous media such as unconsoli-
dated oil sands. We hope that our work will be an incentive to initiate the application of pore scale
modeling in those studies.
097106-3 Askari, Taheri, and Hejazi AIP Advances 5, 097106 (2015)

FIG. 1. Swelling circles algorithm.

II. GENERATION OF POROUS MEDIA

We develop a computer program based on the swelling circles algorithm43 (Figure 1) to recon-
struct granular porous media whose grains are circle shape. At the first step, the dimension of a
medium in which a porous media is simulated is determined (the dashed square in Figure 2). Then,
based on the grain size distribution of a porous medium, a set of random points in a 2D plane are
created (Figure 2(a)). Each point corresponds to a circle center to which a radius growth rate ratio
is assigned with respect to its radius. This parameter determines how fast the circle will grow rela-
tive to other circles. Thereafter, circles start growing (Figure 2(b)-2(f)). During the growing stage,
some circles might overlap each other. Therefore, the circle centers are slightly rearranged to remove
overlaps. Figure 2(c) shows a situation where two circles overlap, and the overlap is taken away by
small reconfiguration of the centers in Figure 2(d). In some situations, in particular in the last steps
of swelling, overlaps might not be resolved by scores of reconfigurations (Circle A in Figure 2(e)).
In this case, after N iterations which has been defaulted in the program, some circles will be slightly
shrunk so that they can fit into the packing space (Circle A in Figure 2(f)). The growing and shrinking
of the circles will continue until the circles radii reach their assigned values (Figure 2(f)).

III. CALCULATION OF THE EFFECTIVE THERMAL CONDUCTIVITY

In a porous medium, which is a composition of solid (matrix) and fluid (e.g., gas, water,
and oil), we consider an overall thermal conductivity for the entire medium. This overall thermal
097106-4 Askari, Taheri, and Hejazi AIP Advances 5, 097106 (2015)

FIG. 2. The steps of swelling circles algorithm. (a) Selection of random points, (b) initial growth of circles, (c) an example
of the overlap of two circles, (d) the removal of the overlap in Fig 2(c) by reconfiguring the circles, (e) one of the last steps
of the growth, (f) the final step of swelling after slightly shrinking the circle A.

conductivity is called the effective thermal conductivity and is a function of the thermal conduc-
tivities of matrix and fluid, and the configurations of grains and pores.16 In this study, by solving
the steady state heat conduction equation, we determine temperature distributions for 2D granular
porous media. Using the calculated 2D temperature profile, we estimate the heat transfer rate and
consequently the effective thermal conductivity.
We digitize a 2D reconstructed porous medium into a pixel based 2D image. Each pixel is given
either zero or one value based on its correspondence (pore or grain). Finite volume method48 has
been widely used in numerical heat transfer studies (e.g., Ref. 49). In addition, finite volume method
is independent from the thermal conductivity of matrix to fluid ratio, which provides accurate re-
sults when the thermal conductivity of fluid is low (for example, gas).50 Therefore, we adopt this
technique to solve the steady state 2D heat conduction equation with no heat flow through upper and
lower boundaries and constant temperatures at the right and left boundaries as shown in Figure 3.
We obtain a set of linear equations to compute temperature distribution. The system of equations
can be written as51
Gm = d, (1)
where the coefficient matrix, G, contains the thermal conductivities assigned to the grids, the vector
m is composed of unknowns parameter (here the temperature profile of the grids), and d is the
right constant vector which is determined by the boundary condition (here the temperatures at the
right and left sides of a medium). We use the Bi-Conjugate gradient stabilized method,52 which is
097106-5 Askari, Taheri, and Hejazi AIP Advances 5, 097106 (2015)

FIG. 3. Schematic of the domain of solution and boundary conditions.

a fast and smoothly converging variant of the BiConjugate Gradient method, to solve our system of
equations.
We discretize the domain using uniform grids in both x and y directions. The grid sizes are
chosen to be equal or smaller than the image pixel sizes which ensures the conservation of image
resolution. Hence, surface roughness effects are properly captured. Moreover, this guarantees that
the simulation results are grid independent. To demonstrate the grid intendancy, two simulations
with different grid sizes are performed (Figure 4) for the reconstructed porous medium in Fig-
ure 11(a). With respect to Figure 4, we observe a perfect match for the temperature contour lines

FIG. 4. Temperature contours for a porous medium using two different grid sizes: 5 µm × 5 µm (solid lines) and 2.5 µm
× 2.5 µm (hashed lines).
097106-6 Askari, Taheri, and Hejazi AIP Advances 5, 097106 (2015)

1000 × 1000 grids with the size of 5 µm × 5 µm (solid lines) and 2000 × 2000 grids with the size of
2.5 µm × 2.5 µm (hashed lines). Accordingly, in this study for every simulation, we consider grids
to be small enough so that the calculated effective thermal conductivity becomes independent from
the grid size.

A. Thermal contact resistance simulation

Surface roughness of grains in a porous medium results in the imperfect contacts of the grains
which consequently hinders the heat transfer from one grain to another. This phenomenon known
as thermal contact resistance is analytically modeled by assuming that a temperature jump (∆T ′)
exists at the solid-solid contact surface (Fig. 5). Therefore, the overall thermal contact resistance, R,
stemming from imperfect contact is given by53
R = ∆T ′/q, (2)
′
where q denotes the heat flow rate, and ∆T is the temperature difference between the two adjacent
solid-solid contact surface. The temperature difference, ∆T ′, can be obtained from54
∆T ′ = γcD (1 + λ) (3)
where cD is the displacement constant and has a units of length which depends on the surface
geometry, λ is the ratio of the thermal conductivities of the media in the left and right sides of the
boundary (k1/k2), and γ is the ratio of the flux to the medium thermal conductivity (q/k1 ). ∆T ′
depends upon the geometry of roughness and the thermal conductivity of the fluid filling roughed
area. To explain equation (3) and demonstrate the role of surface roughness on the overall thermal
conductivity, we obtain the temperature gradient, ∆T ′, for a void space in Figure 6 with dimensions
of l 1 and l 2. The void space is centered in the boundary between two media of the same dimensions
L 1/2 and L 2 whose thermal conductivities are k1 and k2 respectively. The heat flux is from the left to

FIG. 5. The micro structure of the contact surface (bottom) between tow grains.
097106-7 Askari, Taheri, and Hejazi AIP Advances 5, 097106 (2015)

FIG. 6. A medium where in the boundaries of k 1 and k 2, there is a void with thermal conductivity of k.

right side of the medium. We assume the void to be filled by a fluid with the thermal conductivity of
k. We divide the medium into three sections of a, b, and c. The overall thermal conductivity of the
system, k o is obtained from
1
ko = 1−ε 2 ε2
(4)
α + (1−ε 1)α+kε 1

where ε 1 = Ll1 , ε 2 = Ll2 and α = k2k+k

1k 2
which is the medium overall thermal conductivity without
1 2 1 2
the void space. Now consider two cases (1) where the system does not have a void space, i.e., no
roughness and (2) there is surface roughness as depicted in Figure 6. In order to have similar heat
flux through these media, an additional temperature drop is resulted for the case of roughed surfaces
as a result of thermal contact resistant. This additional temperature drop, ∆T, is calculated as:
ε 1ε 2 (k − α)
∆T ′ = ∆T (5)
ε 1 (k − α) + α
where ∆T is the temperate difference between the left and right boundaries for the case of no void
space. As expected, Equation (5) shows that the additional temperature drop will disappear for zero
dimension void space ε 1 = 0, ε 2 = 0 or when the void is equally filled with similar materials as k1
and k 2, i.e. k = α
For the above example, the coefficient is expressed by
.
ε 1ε 2(k − α)
CD = L 1 (6)
(ε 1 (k − α) + α)(2k1/α)
Note that if l 1 → 0 or l 2 → 0, CD → 0, and consequently, ∆T ′ → 0.
The surface roughness has complex geometries which makes it impossible to be analytically
incorporated in real porous medium cases similar to what we calculated in Equations (5) and (6).
To express the importance of roughness simulation for the thermal conduction modeling of porous
media, we simulate porous media with respect to an experimental study on Utah oil sand samples.13
In our simulations, we assume that the thermal conductivity of matrix (solid) is 7.4 Wm−1K−1,
and pores are filled by air wimake th the thermal conductivity of 0.0262 Wm−1K−1. The grain
diameter is assumed to be 250 µm with a standard deviation of 30 µm. We consider two cases. In
the first case, no roughness is assumed for the grains, and therefore, the thermal contact resistance
is minimum. We reconstruct porous media for more than 550 realizations for porosities varying
from 0.18 to 0.50. The solid points in Figures 7 are the predicted values of the effective thermal
conductivity, and the empty circles represent experimental data. We fit second degree polynomials
to the predicted data (the black solid lines in Figure 7) to obtain a general trend of the effec-
tive thermal conductivity with respect to porosity, and compare it with the experimental data (the
097106-8 Askari, Taheri, and Hejazi AIP Advances 5, 097106 (2015)

FIG. 7. The estimated effective thermal conductivities of the reconstructed Utah oil sand sample (the solid points) with re-
spect to the experimental data (empty circles represents),13 and other correlations. 1: Geometrical mean value, 2: Receiprocal,
3: EMT model, 4: Maxwell model, and 5: Veerenda and Chaudhary model.

TABLE I. Some well-known analytical and experimental models for the effective thermal conductivity of two-component
materials. In the above equations α = k f /k s .

kf ϕ
( )
1 ke = k s ks Weighed geometric mean55
√
ke 1+( α−1)ϕ
2 ks = √
1+( 1/α−1)ϕ
Reciprocity model55

k −k e
3 −k e
(1 − ϕ) kkes+2k e
+ k ff+2k e = 0 ETM model55

ke 3(α−1)ϕ
4 ks = 1 + (α+2)−(α−1)ϕ Maxwell model55
( )1
ks 3
5 k e = (1 − ϕ)k H + ϕk L − ϑ kf Veerendra and Chaudhary model56
where
k H = k s e βϕ
k L = k f e −β(1−ϕ)
k
β = k sf − 1
and, the coefficient ϑ is adjusted with
respect to the experimental results.

empty circles). Figure 7 demonstrates that our predicted thermal conductivities for the case when
no roughness is presumed are overestimated with regard to the experimental measurements. In
fact, neglecting the surface roughness results in zero thermal contact resistance and consequently
over-predicted effective thermal conductivities. In addition, we compare the experimental data,
with some well-known experimental and analytical models such as weighted geometric mean,55
reciprocal model,55 ETM,55 Maxwell model,55 and Veerendra and Chaudhary model56 (Table I).
None of the models present a strong correlation with the experimental data. However, the weighted
geometric mean exhibits a more acceptable coherency with the experimental data.
Various mathematical and experimental analyses have shown that the surface roughness results
in the thermal contact resistance.57–59 Thermal contact resistance becomes very important in heat
conduction through unconsolidated granular porous media because imperfect contact points among
grains dramatically affect the thermal behavior of a medium. In some studies, this phenomenon has
been taken into account by assuming that the grains are covered by a thin layer of material with a
thermal conductivity lower than that of grains (e.g., Refs. 60 and 61). In our study, i.e. pore level
modeling of heat flow in granular porous media, we capture the influence of imperfect grain-grain
contacts in micro scales by inclusion of roughness to the grains’ surfaces. To this aim, around any
circle, we consider an infinitesimally thin layer to be called “outer layer” in this paper (the gray
layer in Figure 8; however, its thickness is widely shown for visualization purpose). By dividing
097106-9 Askari, Taheri, and Hejazi AIP Advances 5, 097106 (2015)

FIG. 8. Roughness simulation, by randomly removing a fraction of the outer layer (the gray ring of the circle A) which
results in circle B.

FIG. 9. The estimated effective thermal conductivities of the reconstructed Utah oil sand sample (the solid points) after
applying roughness to grains. Empty circles represents experimental data.13

the layer into many small sections and randomly removing a percentage of them, we obtain a circle
whose surface is not smooth (circle B in Fig. 8). Thus, by repeating the process for all circles in a
porous medium, roughness is simulated which accordingly results in thermal contact resistance.
Using the above approach, we include roughness to our simulations of Utah oil sand by
randomly removing 50% of the grains’ outer layers (Figure 9). The solid line in Figure 9 is a second
degree polynomial fitted to the simulations. When the roughness is considered, our simulations
are in excellent agreement with the experimental measurements. Moreover, our simulation results
exhibit a better coherency with the experimental data in comparison with the weighted geometric
mean model. This demonstrates the potential of the pore scale modeling to characterize the ther-
mal conduction mechanism in granular porous media when their micro structural features, such
as roughness, are taken into account. From Figures 7 and 9, it can be concluded that roughness
simulation is an important factor for the pore scale modeling approach to obtain reliable predictions
of the effective thermal conductivities.

B. Alberta oil sand sample reconstruction

We reconstruct an Alberta oil sand sample based on its grain size distribution (Figure 10). Fig-
ure 11(a) shows the reconstructed medium in which 10% of the sections of the grains’ outer layers
have been randomly removed to simulate roughness. The medium dimension is 5 mm × 5 mm and
the porosity is 0.16. To investigate the roughness influence on the effective thermal conductivity, we
calculate the effective thermal conductivity of the medium by assuming that the thermal conduc-
tivity of the matrix is 7.4 Wm−1K−1 and the fluid saturated pores is 0.0262 Wm−1K−1.62 The heat is
transferred from the left to the right side (x direction) where the temperatures for the left and right
sides are 200◦C and 20◦C, respectively, while the upper and lower sides are assumed to be adiabatic.
097106-10 Askari, Taheri, and Hejazi AIP Advances 5, 097106 (2015)

FIG. 10. Grain size distribution of an Alberta, Canada oil sand sample.

FIG. 11. (a) A porous medium reconstructed based upon the grain size distribution in Fig. 10, and (b) its temperature
distribution.
097106-11 Askari, Taheri, and Hejazi AIP Advances 5, 097106 (2015)

FIG. 12. The temperature distributions of the reconstructed Alberta oil sand sample with respect to different roughnesses.
(a) 20% (b) 30% (c) 40% (d) 50%.

Figure 11(b) shows the temperature distribution for the medium (we recall this as the reference
temperature distribution). The estimated effective thermal conductivity is 3.87 Wm−1K−1.
We gradually increment roughness by removing more percentages (20%, 30%, 40%, and 50%)
of the outer layers’ sections. Figures 12(a)-12(d) show the temperature distributions for the reduced
outer layers’ sections, 20%, 30%, 40%, and 50% respectively. The circles in Figure 13 are the
estimated effective thermal conductivities with respect to the percentage of the removed sections.
The effective thermal conductivity decreases (almost linearly) when more sections of the outer
layers are discarded which consequently culminates in increasing roughness and the thermal con-
tact resistance among the grains. To better comprehend the thermal contact resistance effect and
its dependency on roughness, we subtract the reference temperature distribution from the temper-
ature distributions in Figure 12. Figures 14(a)-14(b) show the difference between the reference
temperature distribution and the temperature distributions in Figure 12. We observe the increase
of the thermal contact resistance with roughness as the temperature differences in the boundaries
of the grains augment accordingly. Therein lies the reason why the effective thermal conductivity
decreases with increased roughness as it hampers the heat transfer from one grain to another.
In addition, we estimate the effective thermal conductivity for the reconstructed Alberta oil
sand sample by assuming different k s /k f ratios (1, 4, 12, and 50). When k s /k f is 1, the effective
097106-12 Askari, Taheri, and Hejazi AIP Advances 5, 097106 (2015)

FIG. 13. The estimated effective thermal conductivities of the reconstructed Alberta oil sand sample with respect to different
roughnesses when pores are filled by air (circles), and water (squares).

FIG. 14. (a), (b), (c), and (d) the differences of the temperature distribution in Fig. 11(b) from the temperature distributions
in Fig. 12 respectively.
097106-13 Askari, Taheri, and Hejazi AIP Advances 5, 097106 (2015)

FIG. 15. The effective thermal conductivity with respect to grains diameters (3.75, 7.5, 15, 30, and 60).

thermal conductivity exhibits no dependency on roughness as there is no difference in the rate of

the thermal energy transferred through either solid or fluid parts. However, by increasing the k s /k f
ratio, we observe that the effective thermal conductivity decays more rapidly with roughness. This
decay can be explained by the effect of thermal contact resistance because when there is a higher ra-
tio of k s /k f , thermal energy is mainly transferred through the solid part. In this case, the roughness
can hinder thermal energy in the boundaries among grain, and hence decreases the effective thermal
conductivity in a porous medium.

C. Grain size
To investigate the effect of the grain size on the effective thermal conductivity, we estimate
the effective thermal conductivity with respect to grains’ diameters of 3.75 µm, 7.5 µm, 15 µm,
30 µm, and 60 µm as shown in Figure 15. Since the influence of grains size is explained by the
thermal contact resistance effect, we assume k s /k f to be 1000 where the most of thermal energy is
transferred through the solid part of a porous medium, and hence, the roughness effect and thermal
contact resistance are more observable. We consider two different porosities of 0.36, and 0.41. For
both porosities, we observe the same trend of increment of the effective thermal conductivity (k e )
with grain size. In addition, the trend intensifies with increased grain diameter (in particular after
30 µm). Our simulated trends of the effective thermal conductivity with respect to grain size are in
strong agreement by experimental observations by Midttømme and Roaldset.16

IV. DISCUSSION
We show how to implement micro structures of granular porous media such as grains size
distribution to estimate the effective thermal conductivity as a macroscopic property. Specifically,
we focus on roughness as a decisive parameter to characterize thermal behavior of unconsoli-
dated granular porous media. Roughness causes imperfect surface contacts among grains which
consequently hampers heat transfer in a porous medium. Therefore, to obtain realistic estimates
of granular porous medium properties, a meticulous simulation of roughness is necessary. For
instance, the grain size influence on the effective thermal conductivity is explained by the thermal
contact resistance effect. In a porous medium with small grain size, there are a larger number
of thermal energy transfer hindering zones along the path of heat flow which accordingly culmi-
nates in decreasing thermal conductivity.16 Therefore, the success of any computational modeling
to perceive the grain size effect depends on the thermal contact resistance simulation. Quite the
contrary to the experimental results, Wang et al.47 propose that the effective thermal conductivity
decreases with increased grain size. We believe their inconsistent result is due to not considering the
thermal resistance effect in their simulation.
097106-14 Askari, Taheri, and Hejazi AIP Advances 5, 097106 (2015)

In this study, we utilize a simple method for the simulation of roughness which is based on
randomly removing some outer layers’ sections. However, to grant better quantitative prediction of
the effective thermal conductivity, an accurate simulation of porous media with respect to fractural
geometry of roughness63 and 3D effects are necessary, which is the theme of an undergoing study.

V. CONCLUSION
In real unconsolidated granular porous media, thermal contact resistance is always present.
This explains why the reconstructed Utah oil sand models do not match with the experimental
measurements when roughness is not considered. According to our results from an unconsolidated
Alberta oil sand sample simulations, the effective thermal conductivity decreases with roughness
in the wake of increasing thermal contact resistance. The thermal contact resistance effect can be
seen in the increase of the temperature gradients in the boundaries among the grain. When pores are
filled by air (high ratio of the thermal conductivity of matrix to fluid), the roughness effect is more
observable because in this case, most of thermal energy is transferred through gains (as the matrix)
and imperfect contacts among them dramatically result in significant energy transport hampering.
Moreover, both experimental studies and our simulations reveal that the effective thermal conduc-
tivity increases with augmented grain size due to the thermal contact resistance effect as there are
less heat transfer hindering points along the heat flux.
The consistencies of our computational results in the Utah oil sand simulation and grain size
effect with the experimental studies demonstrate the promises of the pore scale modeling approach
for further numerical studies in the realm of heat transfer in unconsolidated porous media.

ACKNOWLEDGMENT
The authors are grateful to Drs. H. Motahhari and Y. Haseli for fruitful discussions and valu-
able comments. The authors also acknowledge financial support from the Natural Sciences and
Engineering Research Council of Canada (NSERC).
1 R. M. Butler, G. S. McNab, and H. Y. Lo, “Theoretical studies on the gravity drainage of heavy oil during in-situ steam
heating,” Can. J. Chem. Eng. 59, 455–460 (1981).
2 J. M. Beck and A. E. Beck, “Computing thermal conductivities of rocks from chips and conventional specimen,” J. Geophysc.

Res. 70, 5227-5239 (1965).

3 A. E. Beck, “An improved method of computing the thermal conductivity of fluid-filled sedimentary rocks,” Geophysics

41, 133-144 (1976).

4 F. Yu, G. S. Wei, X. X. Zhang, and K. Chen, “Two effective thermal conductivity models for porous media with hollow

spherical agglomerates,” Int. J. Thermophys. 27, 293-303 (2006).

5 F. Tong, L. Jing, and R. W. Zimmerman, “An effective thermal conductivity model of geological porous media for coupled

thermo-hydro-mechanical systems with multiphase flow,” International Journal of Rock Mechanics and Mining Sciences
46, 1358-1369 (2009).
6 C. Clauser and H. Villinger, “Analysis of conductive and convective heat transfer in a sedimentary basin, demonstrated for

the Rheingraben,” Geophys. J. Int. 100, 393-414 (1990).

7 R. Von Herzen and A. E. Maxwell, “The Measuremenot f Thermal Conductivityo f Deep-Sea Sediments by a Needle-Probe

Method,” Journal of Geophysical Research 64, 1557-1562 (1959).

8 A. Beck, “A steady-state method for the rapid measurement of the thermal conductivity of rocks,” J. $ci. Instr. 34, 186-189

(1957).
9 Y. A. Popov, D. F. C. Pribnow, J. H. Sass, C. F. Williams, and H. Burkhardt, “Characterization of rock thermal conductivity

by high-resolution optical scanning. 28 (2),” Geothermics 28(2), 253–276 (1999).

10 W. Woodside and J. H. Messmer, “Thermal conductivity of porous media,” J. Appl. Phys. 32, 1688–1706 (1961).
11 J. Anand, W. H. Somerton, and E. Gomaa, “Predicting thermal conductivities of formations from other known properties,”

Society of Petroleum Engineers Journal 13, 267-273 (1973).

12 R. Goss, J. Combs, and A. Timur, “Prediction of thermal conductivity in rocks from other physical parameters and from

standard geophysical well logs,” in Proceedings of the SPWLA 16th Annual Logging Symposium (1975).
13 W. R. Lindberg, R. R. Thomas, and R. J. Christensen, “Measurements of specific heat, thermal conductivity and thermal

diffusivity of Utah tar sands,” Fuel 64, 80-85 (1985).

14 G. Vasseur, F. Brigaud, and L. Demongodin, “Thermal conductivity estimation in sedimentary basins,” Tectonophysics 224,

167–174 (1995).
15 K. Midttømme, E. Roaldset, and P. Aagaard, “Thermal Conductivities of Argillaceous Sediments,” in Modern Geophysics in

Engineering Geology Special Publications (Geological Society Engineering Geology, London, 1997), Vol. 12, pp. 355-363.
16 K. Midttømme and E. Roaldset, “The effect of grain size on thermal conductivity of quartz sands and silts,” Petroleum

Geoscience 4, 165–172 (1998).

097106-15 Askari, Taheri, and Hejazi AIP Advances 5, 097106 (2015)

17 S. L. Bryant, P. R. King, and D. W. Mellor, “1993a, Network Model Evaluation of Permeability and Spatial Correlation in
a Real Random Sphere Packing,” Transport in Porous Media 11, 53-70 (1993).
18 S. Mohanty, “Effect of multiphase fluid saturation on the thermal conductivity of geologic media,” J. Appl. Phys. 30,

L80–L84 (1997).
19 G. J. Cheng, A. B. Yu, and P. Zulli, “Evaluation of effective thermal conductivity from the structure of a packed bed,”

Chemical Engineering Science 54, 4199–4209 (1999).

20 J. K. Carson, S. J. Lovatt, D. J. Tanner, and A. C. Cleland, “An analysis of the influence of material structure on the effective

thermal conductivity of theoretical porous materials using finite element simulations,” International Journal of Refrigeration
873–880, 26 (2003).
21 M. R. Wang and N. Pan, “Predictions of effective physical properties of complex multiphase materials,” Materials Science

and Engineering R. 1-30, 63 (2008).

22 M. R. Arab, B. Pateyron, M. El Ganaoui, and N. Calve, “Lattice Boltzamnn simulation for thermal conductivity estimation

in heterogeneous materials,” Defect and Diffusion Forums 283-286, 364-369 (2009).

23 J. Kou, F. Wu, H. Lu, Y. Xu, and F. Song, “The effective thermal conductivity of porous media based on statistical

self-similarity,” Physics Letters A 374, 62–65 (2009).

24 M. G. Alishaev, I. M. Abdulagatov, and Z. Z. Abdulagatova, “Effective thermal conductivity of fluid-saturated rocks: Ex-

periment and modeling,” Engineering Geology 135–136, 24–39 (2012).

25 D. P. Do and D. Hoxha, “Temperature and Pressure Dependence of the Effective Thermal Conductivity of Geomaterials:

Numerical Investigation by the Immersed Interface Method,” Journal of Applied Mathematics 2013, Article ID 456931
(2013).
26 M. Wang, J. Wang, N. Pan, and S. Chen, “Mesoscopic Predictions of the Effective Thermal Conductivity of Microscale

Random Porous Media,” Physical Review E 75, 036702 (2007).

27 M. Wang and Z. Li, “Nonideal gas flow and heat transfer in micro- and nanochannels using the direct simulation Monte

Carlo method,” Physical Review E 68, 046704 (2003).

28 I. Fatt, “The network model of porous media (in three parts),” Pet. Trans. AIME 207, 144-181 (1956).
29 S. Bryant and M. Blunt, “Prediction of relative permeability in simple porous media,” Physical Review A 46, 2004-2011

(1992).
30 A. C. Payatakes, C. Tien, and R. M. Turian, “A new model for granular porous media: Part I. Model formulation,” AIChE

J. 19, 58-67 (1973).

31 F. A. Dullien, Porous Media - Fluid Transport and Pore Structure (Academic Press, New York, 1979).
32 A. C. Payatakes and M. A. Neira, “Model of the constricted unit cell type for isotropic granular porous media,” AIChE J.

23, 922-930 (1977).

33 A. C. Payatakes, K. M. Ng, and R. W. Flumerfelt, “Oil ganglion dynamics during immiscible displacement: model formu-

lation,” AIChE J. 26, 430-443 (1980).

34 C. David, Y. Gueguen, and G. Pampoukis, “Effective medium theory and network theory applied to the transport properties

of rock,” J. Geophys. Res. 95(B5), 6993-7005 (1990).

35 J. M. Dvorkin, M. Armbruster, C. Baldwin, Q. Fiang, N. Derzhi, C. Gomez, B. Nur, and A. Nur, “The future of rock physics:

Computational methods vs. lab testing,” First Break 26, 63–68 (2008).
36 M. A. Knackstedt, S. Latham, M. Madadi, A. Sheppard, T. Varslot, and C. Arns, “Digital rock physics: 3D imaging of core

materialand correlations to acoustic and flow properties,” The Leading Edge 28(no. 1), 28–33 (2009).
37 E. H. Saenger, S. M. Schmalholz, M. A. Lambert, T. T. Nguyen, A. Torres, S. Metzger, R. M. Habiger, T. Müller, S. Rentsch,

and E. Méndez-Hernández, “A passive seismic survey over a gas field: Analysis of low-frequency anomalies,” Geophysics
74(no. 2), O29-O40 (2009).
38 R. M. Holt, E. Fjaer, O. Torsaeter, and S. Bakke, “Petrophysical Laboratory Measurements for Basin and Resrervoir Eval-

uation,” Marine and Petroleum Geology 13, 383-391 (1996).

39 P. E. Øren, S. Bakke, and O. J. Arntzen, “Extending predictive capabilities to network models,” SPE J 3, 324-336 (1998).
40 J. Finney, “Random packings and the structure of simple liquids. I. The geometry of random close packing,” Proc. Roy. Soc.

319A, 479–494 (1970).

41 L. M. Schwartz and J. R. Banavar, “Transport properties of disordered continuum systems,” Phys. Rev. B 39, 11965-11970

(1989).
42 P. H. Valvatne, M. Piri, X Lopez, and M. J. Blunt, “Predictive pore-scale modeling of single and multiphase flow,” Transport

Porous Med. 58(1–2), 23–41 (2005).

43 M. A. Mousavi and S. L. Bryant, “Connectivity of Pore Space as a Control on Two-Phase Flow Properties of Tight-Gas

Sandstones,” Transp Porous Med 94, 537–554 (2012).

44 V. Rzhevsky and G. Novik, The Physics of Rocks (Mir Publishers, Moscow, 1971), translated from the Russian by A. K.

Chatterjee.
45 F. X. Alvarez, D. Jou, and A. Sellitto, “Pore-size dependence of the thermal conductivity of porous silicon, A phonon

hydrodynamic approach,” Appl. Phys. Lett. 97, 33103 (2010).

46 J. J. Zhao, “Thermophysical properties and heat transfer mechanisms of microsclae and nanoscale structures in aerogel-based

composite insulators,” PhD Thesis, Tsinghua University, Beijing, China, 2012.

47 M. R. Wang, N. Pan N., J. Wang, and S. Chen, “Mesoscopic simulation pf phase distribution effects on the effective thermal

conductivity of micro-granular porous media,” J. Colliod Interface Sci 311, 562-570 (2007).
48 E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics. A Practical Introduction, 2nd ed. (Springer-Verlag,

Berlin, 1999).
49 M. A .A. Mendes, S. Ray, and D. Trimis, “A simple and efficient method for the evaluation of effective thermal conductivity

of open-cell foam-like structures,” Int. J. Heat Mass Transfer 66, 412–422 (2013).
50 C. Demuth, M. A. A. Mendes, S. Ray, and D. Trimis, “Performance of thermal lattice Boltzmann and finite volume methods

for the solution of heat conduction equation in 2D and 3D composite media with inclined and curved interfaces,” International
Journal of Heat and Mass Transfer 77, 979–994 (2014).
097106-16 Askari, Taheri, and Hejazi AIP Advances 5, 097106 (2015)

51 R. C. Aster, B. Borchers, and C. H. Thurber, Parameter Estimation and Inverse Problems (Academic Press, New York,
2004).
52 H. A. van der Vorst, “Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear

systems,” SIAM Journal on scientific and Statistical Computing 13, 631-644 (1992).
53 J. P. Holman, Heat Transfer (McGraw-Hill Book Company, Toronto, ON, 1986).
54 M.M. Fyrillas and C. Pozrikidis, “Conductive heat transport across rough surfaces and interfaces between two conforming

media,” International Journal of Heat and Mass Transfer 44, 1789-1801 (2001).
55 M. Wang and N. Pan, “RPredictions of effective phy sical prop erties of complex multiphas e mat erials,” Materials Science

and Engineering 63, 1-30 (2008).

56 A. Aurangzeb and Maqsood, “Modeling of effective thermal conductivity of consolidated porous media with different satu-

rations: A test case of Gabbro rocks,” International Journal of Thermodynamic 28, 1371-1386 (2007).
57 I. Savija, J. R. Culham, and M. M. Yovanovich, “Review of Thermal Conductance Models for Joints Incorporating Enhance-

ment Materials,” Journal of Thermophysics and Heat Transfer 17(1), 43-52 (2003).
58 M. M. Fyrillas and C. Pozrikidis, “Conductive heat transport across rough surfaces and interfaces between two conforming

media,” International Journal of Heat and Mass Transfer 44, 1789-1801 (2001).
59 M. A. Lambert and L. S. Fletcher, “Thermal Contact Conductance of Non-Flat, Rough, Metallic Coated Metals,” ASME

Journal of Heat Transfer 124, 405-412 (2002).

60 Y. T. Feng, K. Han, and D. R. J. Owen, “Discrete thermal element modelling of heat conduction in particle systems: Basic

formulations,” journal of computational physics 227, 5072–5089 (2008).

61 Y. T. Feng, K. Han, and D. R. J. Owen, “Discrete thermal element modelling of heat conduction in particle systems:

Pipe-network model and transient analysis,” Powder Technology 193, 248–256 (2009).
62 N. B. Vargaftik, Tables of the Thermophysical Properties of Liquids and Gases, 2nd ed. (Wiley, New York, 1973).
63 A. Majumdar and C. L. Tien, “Fractal characterization and simulation of rough surfaces,” Wear 136, 313–327 (1990).