CET Volume 86


 
 

 

                                                                    DOI: 10.3303/CET2186211 
 

 
 

 
 
 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Paper Received: 15 September 2020; Revised: 6 February 2021; Accepted: 2 May 2021 
Please cite this article as: De Guido G., Pellegrini L.A., 2021, Phase Equilibria Analysis in the Presence of Solid Carbon Dioxide, Chemical 
Engineering Transactions, 86, 1261-1266  DOI:10.3303/CET2186211 

 CHEMICAL ENGINEERING TRANSACTIONS 
VOL. 86, 2021 

A publication of 

The Italian Association 
of Chemical Engineering 
Online at www.cetjournal.it 

Guest Editors: Sauro Pierucci, Jiří Jaromír Klemeš
Copyright © 2021, AIDIC Servizi S.r.l. 
ISBN 978-88-95608-84-6; ISSN 2283-9216

Phase Equilibria Analysis in the Presence of Solid Carbon 
Dioxide 

Giorgia De Guido*, Laura A. Pellegrini 

GASP - Group on Advanced Separation Processes & GAS Processing, Dipartimento di Chimica, Materiali e Ingegneria 
Chimica “Giulio Natta”, Politecnico di Milano, Piazza Leonardo da Vinci 32, I-20133 Milano, Italy 
 giorgia.deguido@polimi.it 

Carbon dioxide can be found in several gaseous streams, from which it has to be separated in order to meet 
commercial specifications or to comply with environmental regulations. Some examples are natural gas, which 
has to be sweetened to a pipeline-quality gas or to liquefied natural gas, and biogas, which has to be 
upgraded to biomethane or liquefied biomethane. 
In recent years, an intense research has been carried out to develop novel technologies for the separation of 
CO2 from these gaseous streams, considering the new challenges the world has to face. In this respect, the 
need for processing high-CO2 content natural gases for meeting the increased demand for clean energy can 
be mentioned, which requires technologies other than the commercially available ones for a profitable 
exploitation of these low-quality reserves. A growing attention has been devoted to low-temperature/cryogenic 
technologies for this purpose, which requires reliable methods to correctly describe the thermodynamics of 
phase equilibria in the presence of solid CO2 that plays a key role in the design of such processes. This work 
presents a thermodynamic method for the simultaneous stability analysis and multiphase equilibrium 
calculations of CO2 mixtures with hydrocarbon and non-hydrocarbon components. The experimental data 
available in the literature for CO2 frost points and solid-vapour equilibrium conditions have been used to 
validate the proposed method. The calculation results have been also compared with those obtained by using 
the RGibbs tool available in the Aspen Plus® process simulator, obtaining a good agreement between the two 
methods. 

1. Introduction

Single-stage phase equilibrium calculations typically involve specification of a feed composition (global 
composition, zi) and two additional variables, normally selected from temperature (T), pressure (P), vapour 
phase fraction (αV), enthalpy, entropy, or a phase mole fraction. Such specifications, however, have to 
guarantee a unique solution. This is the case for specified T and P, where the solution corresponds to the 
global minimum in the Gibbs energy. A common problem in the calculation of phase equilibria is that the 
number of phases that are present at equilibrium is not known a-priori. Hence, two solution strategies have 
been applied. In the first one, a non-linear programming approach is used to minimize the Gibbs energy 
function with a large number of phases (Castillo and Grossmann, 1981; Lantagne et al., 1988). The second 
approach is to sequentially add a phase in the computations and test the stability of the solution (Michelsen, 
1982). Gupta (Gupta, 1990) developed a new stability criterion for multiphase reacting/non-reacting systems 
that allows the simultaneous calculation of phase stability and flash computations. The proposed solution 
provides a unified set of simultaneous equations that describes phase equilibrium, chemical equilibrium (if this 
applies) and the stability of the system. Gupta (1990) applied the model to compute the phase behaviour of 
systems involving vapour-liquid or vapour-liquid-liquid equilibria. Ballard and Sloan (Ballard and Sloan, 2004) 
have further investigated the model proposed by Gupta (1990) for the simultaneous calculation of stability and 
flash computations to predict equilibrium conditions with gas hydrates. More recently, Tang and co-workers 
(Tang et al., 2019) have extended this model to the solid-liquid-vapour-phase flash calculation of the CH4-CO2 
mixture at a given P, T, and feed gas composition. In their study, the fugacity coefficients of fluid phases (i.e., 
vapour and liquid) are calculated using the GERG-2004 Equation of State (EoS), while the EoS that describes 

1261



the thermodynamic behaviour of solid CO2 is based on the Gibbs free energy method suggested by Jäger and 
Span (Jäger and Span, 2012). However, to our knowledge, this approach has not been applied to the 
simultaneous calculation of phase stability and multiphase equilibria of CO2 mixtures with components other 
than methane. This is of interest for practical applications like low-temperature/cryogenic technologies (De 
Guido and Pellegrini, 2019a) for CO2 removal from low-quality natural gas reserves (Pellegrini et al., 2015) or 
biogas (Qyyum et al., 2020), which are operated at T and P conditions where dry ice can form in 
multicomponent systems. This work contributes to this through the application of the method to the prediction 
of frost points and solid-vapour equilibrium (SVE) conditions of mixtures containing CO2, methane and 
nitrogen, a contaminant that is present in both natural gas (De Guido et al., 2019b) and biogas. 

2. Method

In this section, the theory behind the thermodynamic approach, which has been implemented in a Fortran 
code, is described. 

2.1 Thermodynamic model 

For a system with Nc components and π possible phases, if all π phases are present at equilibrium, the 
following must be true, Eq(1). 

ˆ ˆ
1ˆ ˆ

ir ir ir

ik ikik

f x P
x Pf

φ
φ

= =  for i=1, …, Nc and for k = 1,…, π (1)

In Eq(1), the subscript r refers to a reference phase, îkf to the fugacity of component i in phase k, xik to the

mole fraction of component i in phase k, and îkφ  to the fugacity coefficient of component i in phase k. Eq(1)
can be rewritten making use of the equilibrium constant for component i: 

ˆ
ˆ

ik ir
ik

ir ik

x
K

x
φ
φ

= = for i=1, …, Nc and for k = 1,…, π (2)

As reported by Ballard (Ballard, 2002), it is necessary to seek for an equation that reduces to Eq(2) for phases 
that are present at equilibrium, but not for phases that are not present. Eq(3) can be used for this purpose. 

ˆ 1       
ˆ 1        
ir

ik

if phase k is presentf
if phase k is not presentf

=


<
for i=1, …, Nc and for k = 1,…, π (3)

By multiplying the mole fraction ratio in Eq(2) by the fugacity ratio in Eq(3), Eq(4) is obtained: 

ˆ
exp ln ˆ

ik ik
ik

ir ir

x f
K

x f

 
= − 

  
 for i=1, …, Nc and for k = 1,…, π (4)

Gupta (1990) showed that the natural log of the ratio of fugacities in Eq(4) is equal for all components in a 
given phase k, and referred to this value as the stability of phase k, θk. Hence, rearranging Eq(4): 

kik
ik

ir

x
K e

x
θ= for i=1, …, Nc and for k = 1,…, π (5)

Eq(5) is valid for all phases regardless of that phase’s presence in the system. Gupta (1990) showed that 
defining the mole fraction ratio in this manner is equivalent to minimizing the Gibbs energy of the system 
conditional to: 

0k kk
k k

S
α θ

α θ
= =

+
  for k = 1,…, π (6)

With reference to Eq(6), it can be noticed that: 
• if αk > 0, then phase k is present and θk = 0;
• if αk = 0, then phase k is not present and θk ≠ 0.

1262



By including Eq(5) into the combined overall and components’ material balances and considering the 
difference between the stoichiometric equation for each phase and the one for the reference phase, it is 
possible to derive the objective function in Eq(7). 

1

1

( 1)
0

1 ( 1)

k

j

Nc
i ik

k
i

j ij
j
j r

z K e
E

K e

θ

π
θα=

=
≠

−
= =

+ −



 

for k = 1,…, π (7)

Therefore, the proposed method is based on the solution of the following system of (2π-1) equations in the      
(2π-1) unknowns, namely αk (k = 1, …, π) and θk (k = 1, …, π and k ≠ r). 

1

1

1

0   

( 1)
0

    1, , ( 1)

    1, , ( 1)
1 ( 1)

1

k

j

k k
k

k k

Nc
i ik

k
i

j ij
j
j r

r k
k
k r

fS

z K e
E

K e

or k

for k
θ

π
θ

π

π
α θ

α θ

α
π

α α

=

=
≠

=
≠

= … −

= … −




 = =
 +


− = =
 + −




= −








(8) 

The above system (Eq(8)) is solved at a given set of K-values and composition. At the beginning, this requires 
to have initial estimates for molar phase fractions, αk, for stability variables, θk, for K-values, Kik, and for 
composition of phases, xik. The algorithm is started assuming that all phases are present with an equal 
amount of each, and, therefore, the stability variables are all zero (Ballard, 2002), and that K-values are 
composition-independent. Once the above system is solved, using the Newton-Raphson method, it is possible 
to calculate the mole fractions of each component in each phase using Eq(9) and the K-values removing the 
assumption according to which they were assumed composition-independent. Their expressions are reported 
in Table 1, depending on which phase is taken as the reference one. The fugacity coefficients in the vapour 
and liquid phases have been calculated with the Peng-Robinson Equation of State and the solid vapour 
pressure (Psubl) of CO2 is calculated according to the 6-parameter correlation proposed by Jensen et al. 
(2015). 

1
1 ( 1)

k

j

i ik
ik

j ij
j
j r

z K e
x

K e

θ

π
θα

=
≠

=
+ − for i=1, …, Nc (9)

Table 1: Expressions for Kik depending on which phase is taken as the reference (r) phase 

Kik r = V r = L r = S 

KiV 1
iV

iV

x
x

= =   
ˆ ( , , )
ˆ ( , , )

L
iV i L

V
iL i V

x T P
x T P

φ
φ

= =
x
x

 
2

2

for i CO

( )
   for i COˆ ( , , )

V
i
S
i
V subl
i i
S V
i i V

x
x
x P T
x P T Pφ


= → ∞ ≠


= = =
 ⋅ x

KiL 
ˆ ( , , )
ˆ ( , , )

V
iL i V

L
iV i L

x T P
x T P

φ
φ

= =
x
x

 1iL
iL

x
x

= =  
 

2

2

      for i CO

( )
   for i COˆ ( , , )

L
i
S
i
L subl
i i
S L
i i L

x
x
x P T
x P T Pφ


= → ∞ ≠


= = =
 ⋅ x

KiS 

2

2
,

0         for i CO
ˆ ( , , )

 for i CO
( )

V
iS i V

iV subl i

x P T P
x P T

φ

= ≠


⋅
= = =


x
2

2
,

0   for i CO
ˆ ( , , )

=    for i CO
( )

L
iS i L

iL subl i

x P T P
x P T

φ

= ≠


⋅
= =


x  1
S
i
S
i

x
x

= =  

1263



2.2 Application of the thermodynamic model to the prediction of frost points and SVE conditions 

In this work, the calculation of CO2 frost points and SVE conditions has been carried out, respectively, for the 
CO2-CH4 mixture and for the CO2-CH4-N2 mixture. Calculations have been also performed using the RGibbs 
tool available in Aspen Plus® V9.0 (AspenTech, 2016) that, to our knowledge, is the only unit operation able to 
deal with phase equilibria also in the presence of a solid phase if properly set-up for this type of phase 
equilibrium calculations (Pellegrini et al., 2020). As for the calculation of CO2 frost points for the binary 
mixture, the same P and global composition as the experimental ones (Pikaar, 1959; Agrawal and Laverman, 
1995; Le and Trebble, 2007; Zhang et al., 2011) have been specified and the temperature has been varied so 
that the highest value for which the CO2 solidification ratio (defined as the ratio between the solid molar phase 
fraction and the CO2 mole fraction in the feed stream) is less than or equal to 0.001 has been taken as the T 
of frost point. Such calculated temperature (Tcalc) is compared with the one reported in the literature (Texp). As 
for the calculation of SVE conditions for the ternary mixture, the experimental data (Xiong et al., 2015) are 
given in terms of T, P, mole fraction of N2 in the ternary mixture (i.e., 3 mol % or 5 mol %) and composition of 
the vapour phase at equilibrium. In this case, the global composition, which is required as input data in 
addition to P and T, has been assigned considering the given mole fraction of N2 and a mole fraction of CO2 
greater than the one in the vapour phase at equilibrium. The performances of the two methods have been 
assessed by comparing the calculated and the experimental CO2 mole fraction in the vapour phase. 

3. Results and discussion

In this section, results are illustrated and discussed taking into account the average absolute deviation 
(AAD%) calculated according to Eq(10) in order to compare the performances of the two methods. 

.
, ,

1 ,

100
%

.

No points
calc j exp j

j exp j

variable variable
AAD

No points variable=

−
=  (10) 

In Eq(10), variablecalc,j and variableexp,j respectively denote the calculated and experimental values for the 
variable of interest for the j-th point, which is the temperature in the case of frost point calculations and the 
CO2 mole fraction in the vapor phase in the case of SVE calculations. 

3.1 CO2-CH4 binary mixture 

Figure 1 shows the results obtained when using the two methods for the calculation of CO2 frost points of 
different CO2-CH4 mixtures, according to the literature sources. The results are summarized in Table 2 in 
terms of AAD%. The calculated values show a good agreement with the experimental ones and suggest both 
methods are conservative, in most cases (as shown in Figure 1b and in Figure 1d), in predicting the 
temperature at which CO2 freezes out. Higher deviations (Table 2), though still acceptable, have been 
obtained considering the data of Le and Trebble (Le and Trebble, 2007). To better understand the reason for 
this, all the available experimental data have been reported on the same plot (not shown) for assessing the 
consistency of each data set with the other ones. It has been observed that at higher pressures and lower 
amounts of carbon dioxide in the initial mixture (ca. 1 mol %), where the largest disagreement exists, the data 
by Le and Trebble are close to the data by Pikaar (Pikaar, 1959) and are shifted to the right with respect to the 
data by Agrawal and Laverman (Agrawal and Laverman, 1995). On the contrary, at higher amounts of carbon 
dioxide (ca. 3 mol %) the experimental data by the different literature works are very close to each other. 
However, the experimental data provided by Le and Trebble (Le and Trebble, 2007) cover higher pressures (9 
- 25 bar) than those by Pikaar (Pikaar, 1959)      (2 - 18 bar), so it is not possible to safely state whether some 
of the data available in the literature are not reliable. 

Table 2: AAD% (Eq(10)) for the CO2 frost temperature calculated with the two methods considering the 
different literature sources for the experimental data 

Proposed method RGibbs tool 
Pikaar (1959) 0.354 0.633 
Agrawal and Laverman (1995) 1.079 0.873 
Le and Trebble (2007) 1.239 1.364 
Zhang et al. (2011) 0.297 0.375 

1264



a) b)

c) d)

Figure 1: Comparison between the results obtained with the proposed approach (solid lines) and with the 
RGibbs tool (AspenTech, 2016; dashed and dotted lines), and the experimental frost point data for the CO2- 
CH4 system for different CO2 contents (as specified in the labels) as reported by: a) Pikaar (1959); b) Agrawal 
and Laverman (1995); c) Le and Trebble (2007); d) Zhang et al. (2011) 

3.2 CO2-CH4-N2 ternary mixture 

By comparing SVE data at two nitrogen contents (3 and 5 mol %), Xiong et al. (2015) investigated the effect of 
nitrogen on SVE conditions in the CO2-CH4-N2 mixture. Figure 2 shows the results obtained with the two 
methods. Both of them exhibit a good agreement with the experimental data, to a higher extent for the 
proposed method (AAD% = 11.44 %) rather than for the RGibbs tool (AAD% = 15.17 %). In particular, the 
curves confirm the observations made on the basis of the experimental data (Xiong et al., 2015), namely that 
the nitrogen addition, at least up to 5 mol %, has little effect on CO2 freezing conditions, even if with its 
increasing content the maximum pressure increases enabling the mixture to keep in the solid-vapour region at 
higher pressures. 

a) b)

Figure 2: Comparison between the results obtained with the proposed approach (solid lines), with the RGibbs 
tool (AspenTech, 2016; dashed and dotted lines) and the experimental data by Xiong et al. (2015) at 
increasing temperatures (as indicated by the arrow) of 153.15 K, 168.15 K, 178.15 K, 188.15 K, 193.15 K for 
the:     a) CO2-CH4-3 mol % N2 mixture; b) CO2-CH4-5 mol % N2 mixture

1% CO2

5% CO2

3% CO2

20% CO2

0

5

10

15

20

25

30

35

40

45

150 155 160 165 170 175 180 185 190 195 200 205 210 215

P
re

s
s

u
re

 [
b

a
r]

Temperature [K]

0.12% CO2

0.97% CO2

1.8% CO2 3.07% CO2

10.67% CO2

0

5

10

15

20

25

30

130 150 170 190 210

P
re

s
s

u
re

 [
b

a
r]

Temperature [K]

1% CO2

2.93 % CO2

1.9% CO2

0

5

10

15

20

25

30

35

165 170 175 180 185 190

P
re

s
s

u
re

 [
b

a
r]

Temperature [K]

33.4% CO2

17.8% CO2

10.8% CO2

54.2% CO2

42.4% CO2

0

5

10

15

20

25

30

35

40

45

50

190 195 200 205 210 215

P
re

s
s

u
re

 [
b

a
r]

Temperature [K]

0

5

10

15

20

25

0 5 10 15 20 25 30 35 40

P
re

s
s

u
re

 [
b

a
r]

mol % CO2

T

0

5

10

15

20

25

0 5 10 15 20 25 30 35 40

P
re

s
s

u
re

 [
b

a
r]

mol % CO2

T

1265



4. Conclusions

This work deals with a novel thermodynamic method for the simultaneous stability analysis and multiphase 
equilibrium calculations of CO2 mixtures with hydrocarbons and non-hydrocarbon components, taking into 
account the solid phase in addition to fluid phases. The proposed method allows performing isothermal-
isobaric flash computations in multiphase systems at given temperature, pressure and their global composition 
without knowing a-priori the number and the type of phases present at equilibrium. It has been shown that it is 
able to reliably predict frost point temperatures and the vapour phase composition at solid-vapour equilibrium 
conditions, respectively, for the CO2-CH4 and CO2-CH4-N2 mixtures, with average absolute deviations lower 
than 1.3 % and 11.5 % in the two cases. The application of the proposed method is useful for the study and 
correct design of the recently developed CO2 low-temperature/cryogenic removal processes. 

References 

Agrawal G.M., Laverman R.J., 1995, Phase behavior of the methane-carbon dioxide system in the solid-vapor 
region, Advances in Cryogenic Engineering, 19, 327-338. 

AspenTech, 2016, Aspen Plus®, AspenTech, Burlington (MA), United States. 
Ballard A.L., 2002, A non-ideal hydrate solid solution model for a multi-phase equilibria program, PhD 

Dissertation, Colorado School of Mines, USA. 
Ballard A.L., Sloan Jr E.D., 2004, The next generation of hydrate prediction: Part III. Gibbs energy 

minimization formalism, Fluid Phase Equilibria, 218, 15-31. 
Castillo J., Grossmann I.E., 1981, Computation of phase and chemical equilibria, Computers & Chemical 

Engineering, 5, 99-108. 
De Guido G., Pellegrini L.A., 2019a, Application of Conventional and Novel Low-temperature CO2 Removal 

Processes to LNG Production at Different CO2 Concentrations in Natural Gas, Chemical Engineering 
Transactions, 74, 853-858. 

De Guido G., Messinetti F., Spatolisano E., 2019b, Cryogenic Nitrogen Rejection Schemes: Analysis of Their 
Tolerance to CO2, Industrial & Engineering Chemistry Research, 58(37), 17475-17488. 

Gupta A.K., 1990, Steady state simulation of chemical processes, PhD Dissertation, University of Calgary, 
Calgary, Canada. 

Jäger A., Span R., 2012, Equation of state for solid carbon dioxide based on the Gibbs free energy, Journal of 
Chemical & Engineering Data, 57, 590−597. 

Jensen M.J., Russell C.S., Bergeson D., Hoeger C.D., Frankman D.J., Bence C.S., Baxter L.L., 2015, 
Prediction and validation of external cooling loop cryogenic carbon capture (CCC-ECL) for full-scale coal-
fired power plant retrofit, International Journal of Greenhouse Gas Control, 42, 200-212. 

Lantagne G., Marcos R., Cayrol B., 1988, Computation of complex equilibria by nonlinear optimization, 
Computers & Chemical Engineering, 12, 589-599. 

Le T.T., Trebble M.A., 2007, Measurement of Carbon Dioxide Freezing in Mixtures of Methane, Ethane, and 
Nitrogen in the Solid-Vapor Equilibrium Region, Journal of Chemical and Engineering Data, 52, 683-686. 

Michelsen M.L, 1982, The isothermal flash problem. Part I. Stability, Fluid Phase Equilibria, 9, 1-19. 
Pellegrini L.A., Langè S., Baccanelli M., De Guido G., 2015, Techno-economic analysis of LNG production 

using cryogenic vs conventional techniques for natural gas purification, Offshore Mediterranean 
Conference and Exhibition, 25-27 March 2015, Ravenna, Italy. 

Pellegrini L.A., De Guido G., Ingrosso S., 2020, Thermodynamic Framework for Cryogenic Carbon Capture, 
Computer Aided Chemical Engineering, 48, 475-480. 

Pikaar M.J., 1959, A study of phase equilibria in hydrocarbon-CO2 systems, PhD Dissertation, Imperial 
College of Science and Technology, London. 

Qyyum M.A., Haider J., Qadeer K., Valentina V., Khan A., Yasin M., Aslam M., De Guido G., Pellegrini L.A., 
Lee M., 2020, Biogas to liquefied biomethane: Assessment of 3P's–Production, processing, and prospects, 
Renewable and Sustainable Energy Reviews, 119, 109561. 

Tang L., Li C., Lim S., 2019, Solid−Liquid−Vapor Equilibrium Model Applied for a CH4−CO2 Binary Mixture, 
Industrial and Engineering Chemistry Research, 58, 18355-18366. 

Xiong X., Lin W., Jia R., Song Y., Gu A., 2015, Measurement and Calculation of CO2 Frost Points in CH4 + 
CO2/CH4 + CO2 + N2/CH4 + CO2 + C2H6 Mixtures at Low Temperatures, Journal of Chemical and 
Engineering Data, 60, 3077−3086. 

Zhang L., Burgass R., Chapoy A., Tohidi B., Solbraa E., 2011, Measurement and Modeling of CO2 Frost 
Points in the CO2-Methane Systems, Journal of Chemical and Engineering Data, 52, 2971–2975. 

1266