Supporting Information

Four Nanoscale-Extended Equations of State: Phase Behaviour of Confined Fluids

in Shale Reservoirs

by

Kaiqiang Zhang1, Paitoon Tontiwachwuthikul2, Na Jia1*, and Songyan Li3

1
Petroleum Systems Engineering
Faculty of Engineering and Applied Science
University of Regina
Regina, Saskatchewan S4S 0A2
Canada

2
Clean Energy Technologies Research Institute (CETRI)
Faculty of Engineering and Applied Science
University of Regina
Regina, Saskatchewan S4S 0A2
Canada

3
College of Petroleum Engineering
China University of Petroleum (East China)
Qingdao 266580
China

*Corresponding author. Tel: 1-306-337-3287; Fax: 1-306-585-4855.

E-mail: Na.Jia@Uregina.Ca (N. Jia)

 In the Supporting Information, the detailed derivations for the generalized equation of state,

which includes the confinement effects (i.e., pore radius effect) and intermolecular interactions

(i.e., moleculemolecule and moleculewall interactions) are stated. Fig. S1 shows the schematic

diagrams of the nanopore network model, a single nanopore, and the configuration energy in a

single nanopore.

The canonical partition function from the statistical thermodynamics is shown as follows

[1],

1 3 N N
Q( N ,V , T )   e  Ei ( N ,V ) / kT   qint Z ( N ,V , T ) (S1)
i N!

where N is the number of molecules; V is the total volume; T is the temperature; E is the

overall energy state; k is the Boltzmann constant;  is the de Broglie wavelength,

h2
( ) 0.5 , h is the Planck’s constant, m is the molecular mass; qint is the internal
2mkT

partition function; Z is the configuration partition function, which is expressed as,

Z ( N ,V , T )   e  U(r1 ,r2 ,...,rN )/kT dr1dr2 ...drN (S2)

V

where U is the potential energy of entire system of N number of molecules which positions are

described by ri , i = 1,2,…N, and ri is the distance of separation between molecules.

Since the configurational energy ( E conf ) is expressed as [2],

 ln Z
E conf ( N ,V , T )  E ( N ,V , T )  E ( N ,V , T  )  kT 2 ( ) N ,V (S3)
T

So,

T E conf ( N ,V , T )
ln Z ( N ,V , T )  ln Z ( N ,V , T  )   dT (S4)
 kT 2

where only the hard-core repulsive forces between molecules Z HC (  ) are important when the

2
configuration integral at infinite temperature, Z ( N ,V , T  )  Z HC (  ) . Z HC ( N ,V , T  )  V fN

is defined from the literature [2], where V f is the free volume. Thus, Eq. (S4) can be rewritten

as,
T E conf ( N ,V ,T )
( dT )
Z ( N ,V , T )  V fN e  kT 2
(S5)

The free volume V f can be expressed by using the following simple expression [2],

N
V f  V  N  V  (S6)
 max

where  is the excluded volume per fluid molecule and  max is the molecular density of the full-

distributed fluid. Eqs. (S5) and (S6) are substituted into Eq. (S1) to be,
T E conf ( N ,V ,T )
1 ( dT )
Q( N , V , T )   e  Ei ( N ,V ) / kT
 3 N qint
N
(V  N ) N e  kT 2
(S7)
i N!

It is worthwhile to mention that  max is function of the ratio of the pore size (rp ) and fluid

molecule size ( ) , rp /  , whose specific formulation is shown below [3],

6
 max  (1   ) (S8)
 3
c3 ( 0.5rp /  ) c5 ( 0.5rp /  )
  c1  c2 e  c4 e (S9)

where  is the mean porosity of the porous medium initiated by Mueller (2005) [4] and

ci , i  1, 2, 3, 4, 5 are the numerical coefficients obtained from the curve fitting [3]. Eq. (S7)

6
change to be  max  (1  c1 ) when rp /  tends to be infinite (i.e., bulk phase), so,
 3

6(1  c1 ) 6 (1  c1 ) 6(1  c1 )b
 3 3 3 (S10)
 max  N A

where b is the volume parameter of the cubic EOS and N A is the Avogadro constant.

3
From the statistical thermodynamics [1],

 ln Z ( N ,V , T )
P( N ,V , T )  kT ( ) N ,T (S11)
V

1 3 N N
Given that  qint  constant , Eq. (S7) is substituted into Eq. (S11) to be,
N!
T E conf ( N ,V ,T )
( dT )
 ln[(V  N ) e N 2
kT

]
P( N , V , T )  kT ( ) N ,T
V (S12)
NkT Econf ( N , V , T )
 
V  N V

b
Since   and N  n  N A from the literature [1,5], the first term right-hand side of Eq.
NA

nRT
(S12) is converted to be , where n denotes the moles. It is easily found from Eq. (S12)
V  nb

that E conf ( N ,V , T ) is the key part to explicitly illustrate the relationships of P, V, and T. The

configurational energy ( E conf ) is consist of the configurational energy between molecule and

molecule ( Emolecule
conf
molecule ) as well as between the molecule and wall ( Emolecule wall ) , whose
conf

schematic diagrams are shown in Fig. S2,

E conf ( N ,V , T )  Emolecule
conf
molecule  Emolecule wall
conf
(S13)

In this study, the moleculemolecule interactions ( Emolecule

conf
molecule ) are expressed as [6],

kTN 2 U mm (r1 , r2 ...rN )

conf
Emoleculemolecule   g (r;  , T )dV1dV2 dV3 (S14)
2V 2 kT

where g (r;  , T ) is the pair correlation function for molecules interacting through the potential

U (r ) . In the literature, the pair correlation function at low densities was stated clearly [7],

lim g (r;  , T )  e u ( r ) / kT (S15)

 0

4
Here, C  C ( N ,V , T )  r g (r;  , T )dr. Moleculemolecule interactions U mm (r ) are assumed to

be numerically represented through the Lennard-Jones potential, whose numerical equation is,

 LJ  LJ
U mm (r12 )  4 LJ [( )12  ( )6 ] (S16)
r12 r12

where  LJ is the moleculemolecule Lennard‒Jones energy parameter and and  LJ is the

moleculemolecule Lennard‒Jones size parameter.

Thus, Eq. (S14) is re-written as,

kTn2C U (r1 , r2 ...rN )

E conf
moleculemolecule  2  dV1dV2 (S17)
2V ri  kT

Then, the integral part of Eq. (S17) is solved semi-analytically as,

1 U mm (r1 , r2 ...rN ) 4
 dV1dV2  LJ  LJ
3
f ( A) (S18)
V r12  kT kT

c1 c2
f ( A)  c0  
A A

8
where c0 =  , c1 = 3.5622 , and c 2 =  0.6649 . It should be noted that the value of c0 was
9

calculated by solving Eq. (S5) analytically, while the values of c1 and c 2 are obtained from a

non-linear least-square method. The calculated f (A) values from Eq. (S5) and fitting curve by

tuning c1 and c 2 can be found elsewhere [6]. Thus, the moleculemolecule interactions

molecule ) are presented as,

conf
( Emolecule

c1 c
 an 2C  2n 2C LJ LJ
3
(  2)
A A
 molecule  (S19)
conf
Emolecule
V

On the other hand, the moleculewall interactions U mw (r ) are assumed to be well modeled

5
through the square-well potential, which is stated as follows,

 , rij   sw-ij

U mw (rij )    sw-ij , σ sw-ij  rij  ( sw-ij   sw-ij ) (S20)

0, rij  ( sw-ij   sw-ij )

where rij is the distance between the molecule and wall;  sw-ij is the moleculewall square-well

energy parameter;  sw-ij is the moleculewall square-well size parameter; and  sw-ij is the

moleculewall square-well width of interactions, which also denotes as  p to represent the

moleculewall distance in the following sections. Hence, the moleculewall interactions

 wall ) are expressed as [3],

conf
( Emolecule

wall   N  Fp   sw (S21)
conf
Emolecule

where Fp is the fraction of the confined fluid molecules that interact with the pore wall (i.e., in

the square-well region). The local distributions of fluid molecules interacting with the pore wall

are numerically represented by Fp , which is function of the temperature, fluid density, degree of

confinement, and moleculewall interaction potential [3]. An empirical correlation rather than a

complex theoretical model is capable of describing Fp in an accurate and simple way as follows,


Fp  Fpr  (1  Fpr )(1  e  sw / kT )(1  )θ (S22)
 max

(rp   / 2) 2  (rp   / 2   p ) 2
Fpr 
(rp   / 2) 2

rp

p  / 2

where Fpr is the fraction of the random distributed fluid molecules in the square-well region of

6
the pores and  is the geometric term. Thus, Eq. (S22) is substituted into Eq. (S21) to be,

 sw / kT 
 wall  n  N A   sw [ Fpr  (1  Fpr )(1  e )(1  )θ ] (S23)
conf
Emolecule
 max

N N
 ,  max  A
V b

Where  is the volume number density. Hence, the total configurational energy between

molecule and molecule as well as between the molecule and wall is obtained by combing Eqs.

(S13), (S19), and (S23),

c1 c
 an 2C  2n 2C LJ LJ
3
(  2)
E conf ( N , V , T )  A A  n  N   [ F  (1  F )(1  e  sw / kT )(1  nb ) θ ]
A sw pr pr
V V
(S24)

On a basis of Eq. (S12), the generalized expression of the EOS for confined fluids considering

the effects of pore radius, moleculemolecule, and moleculewall interactions is shown as,

nRT n 2C c c
P( N , V , T )   2 [a  2 LJ LJ3
 ( 1  2 )]
V  nb V A A (S25)
nb nb
 (nN A   sw )( 2 )(1  ) 1 (1  Fpr )(1  e  NA sw / RT )
V V

7
References

[1] Abrams DS, Prausnitz JM. Statistical thermodynamics of liquid mixtures: a new

expression for the excess Gibbs energy of partly or completely miscible systems. AIChE J

1975;21:116–28.

[2] Sandler SI. The generalized van der Waals partition function. I. Basic theory. Fluid Phase

Equilib 1985;19:238–57.

[3] Travalloni L, Castier M, Tavares FW, Sandler SI. Thermodynamic modeling of confined

fluids using an extension of the generalized van der Waals theory. Chem Eng Sci

2010;65:3088–99.

[4] Mueller GE. Numerically packing spheres in cylinders. Powder Technol 2005;159:105–10.

[5] Whitson CH, Brule MR. Phase Behavior, Vol. 20. Richardson, Texas Monogr Ser SPE

2000.

[6] Zarragoicoechea GJ, Kuz VA. van der Waals equation of state for a fluid in a nanopore.

Phys Rev E 2002;65:21110.

[7] Islam AW, Patzek TW, Sun AY. Thermodynamics phase changes of nanopore fluids. J

Nat Gas Sci Eng 2015;25:134–9.

[8] Yu Y-X, Gao G-H. Lennard–Jones chain model for self-diffusion of n-alkanes. Int J

Thermophys 2000;21:57–70.

[9] Sharma R V, Sharma KC. The structure factor and the transport properties of dense fluids

having molecules with square well potential, a possible generalization. Phys A Stat Mech

Its Appl 1977;89:213–8.

[10] Zhang P, Hu L, Meegoda JN, Gao S. Micro/nano-pore network analysis of gas flow in

shale matrix. Sci Rep 2015;5:13501.

8
Table S1
Recorded critical properties (i.e., temperature, pressure, and volume), van der Waals equation of state (EOS) constants, Lennard-Jones
and square-well potential parameters of CO2, N2, and C1‒C10 [5,8,9].

Component Tc ( K ) Pc ( Pa ) Vc (m3 / mol) a ( Pa  m 6 / mol 2 ) b (m3 / mol) (/k ) LJ ( K )  LJ (m) (/k ) sw ( K )

CO2 304.2 7.38106 9.4010‒5 1.0210‒48 7.1210‒29 294 2.9510‒10 179
N2 126.2 3.39106 8.9510‒5 3.7810‒49 6.4310‒29 364 3.3210‒10 95
CH4 190.6 4.60106 9.9010‒5 6.3410‒49 7.1410‒29 207 3.5710‒10 105
C2H6 305.4 4.88106 1.4810‒4 1.5410‒48 1.0810‒28 155 3.6110‒10 174
C3H8 369.8 4.25106 2.0310‒4 2.3310‒48 1.5910‒28 120 3.4310‒10 193
i-C4H10 408.1 3.65106 2.6310‒4 3.4210‒48 1.3710‒28 140 3.8510‒10 213
n-C4H10 425.2 3.80106 2.5510‒4 3.8410‒48 1.9410‒28 118 3.9110‒10 211
C5H12 469.6 3.37106 3.0410‒4 5.2710‒48 2.4110‒28 145 3.9610‒10 213
C6H14 507.5 3.29106 3.4410‒4 6.8910‒48 2.9110‒28 199 4.5210‒10 223
C7H16 543.2 3.14106 3.8110‒4 8.5210‒48 3.3810‒28 206 4.7010‒10 231
C8H18 570.5 2.95106 4.2110‒4 1.0410‒47 3.8210‒28 213 4.8910‒10 237
C9H20 598.5 2.73106 4.7110‒4 1.2410‒47 4.4910‒28 220 5.0710‒10 238
C10H22 622.1 2.53106 5.2110‒4 1.4610‒47 5.0710‒28 226 5.2310‒10 245

9
 Porous
medium

Zoom in

r
Nanoscale x
pore

Zoom in

Configuration
energy

Fig. S1 Schematic diagrams of the nano-pore network model [10], nanoscale pore system, and
configuration energy in nanoscale pores in this study.

10
 (a) (b)

Fig. S2 Schematic diagrams of (a) the Lennard‒Jones potential between the molecule and molecule and (b) the square-well potential
between the molecules and wall.

11