Modelling Supercritical Phase Equilibria : Problems and Pitfalls

This article addresses some recurring difficulties and problems of computing phase equilibria involving supercritical fluid phases. These difficulties prevent a full automatization of thermodynamic calculations and require human interference. Examples are the wrong interpretation of experimental data, phase inversion phenomena, or overlooking phases. While none of these insights are knew, publications dealing with them are scattered and sometimes hard to obtain. This article gives a short overview over some of the most common difficulties and pitfalls.


Introduction
While there exist powerful computer programs for the calculation of the thermodynamic properties of fluids, even for phase equilibria of mixtures, these programs usually cannot regarded as "fool-proof " in the sense that they will always: • arrive at a physically reasonable set of model parameters when fitting a model to experimental data, • and find the correct, stable phases when calculating phase equilibria.
All this means that human guidance for such thermodynamic software is-and will be for many years-indispensible.Of course, this insight is not new.But publications addressing these problems are scattered over several decades of literature.This article will, in the limited space available here, address some of the most common problems and pitfalls of phase equilibrium calculations.
2 Solid-fluid equilibria 2.1 Solid compressibility For a stable two-phase equilibrium it is necessary to have equal temperatures, pressures, and chemical potentials for each component in the coexisting phases, In case of fluid phases, the chemical potentials can be calculated from equations of state.For a (pure) solid phase one can assume that (a) the chemical potential is equal where the integrand is the molar volume of the solid as a function of pressure and temperature.It is advisable not to treat this volume as a constant: Many organic compounds have solid phases whose compressibilities are comparable to those of liquids.Examples are the so-called plastic phases of normal alkanes.A better approach is to assume that the isothermal compressibility of the solid is constant, ( )= ( ) There exist better equations of state for solids, e.g., the Murnaghan equation [1, Section 6.1], but these require more experimental input data.Assuming a non-zero compressibility is a minimal requirement for high-pressure applications (above ca. 10 MPa), and even a crude estimate is better than none.

The (carbon dioxide + caffeine) system
The extraction of drugs, dyes, or flavours from biological material by means of supercritical fluids is an established technique.An important application is the extraction of caffeine from green coffee beans with supercritical carbon dioxide.One of the phase equilibria relevant in this context is the solid-fluid equilibrium of the (carbon dioxide + caffeine) system (The real process also involves water).
Here the thermodynamic modelling is complicated by the fact that anhydrous caffeine has two different solid phases.Fig. 1 shows the phase diagram.
The prediction of the melting pressure curve was attempted with three equations of state, the Peng-Robinson equation (PR) [2,3], which is a cubic equation of state, the "simplified perturbed-hard-chain" equation (SPHCT) [4], a simple non-cubic equation with a more realistic hard-sphere term, and the PC-SAFT equation [5].As explained above, the calculation needs the sublimation pressures and the molar volumes of the solid phases.It turned out that none of the publications of the sublimation pressure curves had spotted the solid-solid transition.While there was a good agreement between several authors on the volume of the high-temperature solid phase s 1 , there is a surprisingly wide range of values for the volume of the low-temperature phase s 2 (see Table 1).
The authors of the X-ray studies state that the s 2 phase of caffeine is polymorphic, i.e., a mixture of different crystal structures having almost the same energy [6].
From a thermodynamic point of view, however, this is rather dubious.Moreover, the crystal volumes of the X-ray studies are not compatible with the slope of the s 1 s 2 transition curve, from which the volume change can be obtained via Clapeyron's equation.
The most likely explanation of the crystallographic results is that the samples were cooled too rapidly: The s 1 → s 2 transition takes some time.Indeed, a careful pyknometric study gave a different result that is compatible with the value estimated by transitiometry [8].The molar volume from Wikipedia, which also is used in several computational studies, is wrong; it is probably the molar volume of the liquid phase at the triple point.For a detailed discussion of the various experimental results the reader is referred to a previous publication [9].
The important insight is that it is not always the most sophisticated experimental technique that gives the best results, but that the outcome of experiments depends very much on the sample preparation.
Fig. 2 shows the complete phase diagram of the (CO 2 + caffeine) system, calculated with the PR equation (Fig. 2(a)) and with the PC-SAFT equation (Fig. 2(b)).The diagram contains the vapor pressure curves of the pure fluids (labeled "(lg) i ", with i = 1 (CO 2 ) or 2 (caffeine)), binary critical curves ("l = g"), and three-phase curves.The calculations were made with the ThermoC software package [11].It turns out that the PR equation yields a phase diagram that is in accordance with the experimental data, whereas the PC-SAFT equation predicts a rather distorted phase diagram with two quadruple points.This, however, is not surprising: The molecular model underlying the PC-SAFT equation is a fully flexible chain of tangential hard spheres-a reasonable model for polymers, but not for the caffeine molecule, which is neither spherical nor flexible.
It is not always the most modern equation of state that gives the best results.Particularly with theory-based equations it is advisable to check if the underlying assumptions are justified.

Phase inversion phenomena
Equation ( 1) is symmetric in the sense that there is no preference or indication which of the coexisting phases is the liquid and which is the vapor (or supercritical fluid).If Eq. ( 1) is solved numerically, it is possible that the computer returns the phase compositions in an unexpected order.
In calculations of phase envelopes this problem can be minimized by taking small steps of the control variable; it can be eliminated by integrating the appropriate differential equation of the phase equilibrium [12][13][14].But generally it is necessary to inspect and identify the solutions offered by the computer.

Pycnic inversion
How can we identify the liquid phase?A common answer is that, if there are two coexisting fluid phases, the liquid is the one with the higher molar density.This, however, is not always true.Fig. 3 shows a density-composition diagram of the (hydrogen + hexane) system at 411 K in which the connodes are indicated.The calculations (again based on the PR equation of state) as well as the experimental data [15] indicate that at some point the slopes of the connodes change from negative to positive, which means that the vapor phase attains a higher molar density than the liquid.
The switching of the molar densities of the coexisting fluid phases is called "pycnic inversion".The phenomenon is not rare; in fact, it is rather common for so-called Class III mixtures (see Section 4).For more details the reader is referred to the work of Quiñones-Cisneros [16].

Barotropic inversion
It is tempting to assume that the liquid phase is the one at the bottom of the vessel, and therefore use the mass density (specific gravity) to distinguish vapor and liquid phases.This, however, may be erroneous.Fig. 4 shows an  isothermal phase diagram of the (carbon dioxide + hexadecane) system at 313 K, which means that the carbon dioxide is supercritical.
It turns out that above ca.14 MPa the CO 2 -rich fluid phase attains a higher mass density than the hexadecane-rich liquid phase.A computer program trying to track the phase envelope and using mass densities to identify phases would be "derailed" here.The phenomenon is called "barotropic inversion" or "barotropy".It typically occurs in asymmetric mixtures where the volatile component has small, but heavy molecules, whereas the other component has large, but relatively light molecules.Examples are (CO 2 + alkane) or (SF 6 + alkane) mixtures.Natural oil reservoirs can be inverted if their gas phase is rich in CO 2 and H 2 S (so-called sour gas).A more thorough discussion of barotropy can be found in the literature [16].
Barotropic inversion not only poses a problem to chemical engineers because phases do not flow where they are expected in distillation and extraction columns; it can also hamper the thermodynamic modelling-even though mass usually does not play a role in phase equilibrium calculations.The reason for this is that barotropic inversion can take place in two different ways: Far away from the critical point, where the interfacial tension between the coexisting phases are large, the liquid and the supercritical phase merely change places.Close to the critical point, however, where the interfacial tension is low, the phases may form a fog or an emulsion, and as gravitation is not able to separate the phases, this emulsion can exist for several days.
If an experimental apparatus uses the so-called analytic method, i.e., if samples of the coexisting phases are withdrawn and analyzed, a barotropic two-phase state may be mistaken for a single-phase state.The mismatch and the scatter of the high-pressure experimental data in Fig. 4 suggest a problem of this kind.Fitting interaction parameters of the equation of state to these data would have resulted in grossly distorted phase diagram predictions.
When fitting model parameters to experimental phase equilibrium data, it is important to make sure that the experiments were not affected by barotropic inversion.This may be a problem if the phase equilibrium vessel is not equipped with a window, so that the determination of the phase state is only based on the measurement of the compositions.

Entropic inversion
In contrast to a common belief, one cannot safely assume that the liquid is always the phase with the lower entropy.Fig. 6 shows the (total) molar entropy of the coexisting phases of the (methane + propane) system at 250 K, computed with the PR equation of state.Evidently, the molar entropy of the liquid exceeds that of the vapor for pressures above 15 MPa.
The "entropic inversion" phenomenon can easily be understood by realizing that the molar entropy consists of two contributions-a configurational one, which depends on the molar volume, and an intramolecular one, which depends on the rotational and vibrational degrees of freedom of the molecules.While (in this example) the liquid has the smaller molar volume and hence the lower configurational entropy, it contains more propane-and propane has a higher ideal-gas heat capacity and hence a higher intramolecular entropy.With increasing pressure,  the molar volume of the gas phase drops rapidly, and so its "advantage" of configurational entropy dwindles, until the liquid phase "wins" because of its higher propane content.
This renders the molar entropy unreliable as a criterion for identifying liquid phases.If a distinction between liquid and vapor phases is necessary, the presently best method is probably that of Venkatarathnam [19], who uses the sign of the temperature derivative of the isothermal compressibility, ∂ ∂ ( ) , to distinguish between liquid-like and vapor-like supercritical phases.

Phase diagram classes 4.1 Nomenclature
A computer program designed to calculate two-phase vapor-liquid equilibria will not necessarily discover liquid-liquid or gas-gas equilibria, let alone three-phase states.The question how one can ensure that all relevant equilibria have been found is linked to the classification of fluid phase diagrams.A first attempt was made in the seminal work of Scott and van Konynenburg [20], who studied mixtures described with the van der Waals equation of state and identified five classes, "I" to "V".Their schematic pT diagrams are shown in the top two rows of Fig. 6.Van Konynenburg also added a Class VI to their list, because it was known to exist, although it cannot be obtained with the van der Waals equation.This was somewhat unfortunate, for Class VI is in fact a set of classes.Most of them have been observed experimentally, as has been Class VIII.
Fig. 6 contains many more classes, but is far from being a complete list of the possible phase diagrams.Unfortunately, it is not possible to present the complete history of the phase diagram classification in this article; the reader is referred to the review of van Konynenburg et al. [21] or to textbooks [1,22] and the literature cited therein.
It is evident, however, that the nomenclature of phase diagram classes by van Konynenburg and Scott is not adequate for handling all possible phase diagram classes: It does not provide names for some subclasses of VI and VII, and class names like III-A** do not really describe the pT diagram and are hard to remember.In contrast to this, the rational nomenclature of Bolz et al. [23] is easy to apply.It starts from the observation that, in a pT phase diagram, it is possible to identify paths ("sequences") consisting of alternating critical curves and three-phase curves llg.The simplest sequence is of course a single critical curve.Then a phase diagram class name is constructed as follows: 1.One starts at the pure-fluid critical point with the higher temperature, determines the number of critical curve segments in the sequence originating at that critical point, and indicates the target of the sequence: P for the critical point of the other pure component, Z for low temperatures and pressures ("zero Kelvin"), C for a (hypothetical) compact state, or Q for a quadruple state lllg.The rational name of Class I thus is 1 P (a single critical curve connecting the pure-fluid critical points), and that of Class V is 2 P (a sequence consisting of 2 critical curves, connected by a threephase curve).2. If the target of the first sequence is not P, a symbol for the sequence originating at the other pure-fluid critical point is constructed in the same way.Thus Class III becomes 1 C 1 Z (one critical curve running to high pressures, a sequence consisting of a critical curve and a three-phase curve running to the origin).3.For critical curves which are not part of these sequences, a letter is added that graphically represents the shape of the curve: "l" for a critical curve running from high pressures to a critical endpoint, "n" for a curve running between two endpoints at low pressures (passing through a pressure maximum), and "u" for a critical curve coming from and running to high pressures.Class II thus becomes 1 P l. 4. If there is critical azeotropy, an "A" is appended to the symbol of the affected critical curve.A"Q" can be used to indicate the existence of a quadruple point, if this has not already been mentioned as a sequence target.This way even the variants or subclasses of the Classes VI and VII can be clearly distinguished.

Global phase diagrams
An overview over the phase diagram classes that can be obtained with a given equation of state (plus its set of mixing rules) can be obtained with a so-called global phase diagram.In such a diagram the parameters of an equation of state are treated as freely variable parameters, and each point corresponds to a (normal) phase diagram.An example is shown in Fig. 7, namely the global phase diagram of the Carnahan-Starling-Redlich-Kwong (CSRK) equation of state for equal-sized molecules [24].
The global phase diagram is shown in terms of dimensionless coordinates, The curves in the global phase diagram represent loci of boundary states, e.g., hypothetical phase diagrams marking the transition between two classes.As an example, Fig. 8 illustrates the boundary states via which Class 2 P (= IV) can evolve to Classes 1 P l (= II) via a tricritical point, to 1 C 1 Z (= III) via a double critical endpoint, and to 2 P (= V) via a zero-Kelvin endpoint.There is not enough space here to state the mathematical criteria of these boundary states.But such criteria exist, and from them global phase diagrams like Fig. 7 can be calculated, which give an overview over the possible phase diagram classes.The interested reader is referred to the literature [1,21,24].
Mixtures of similar compounds usually obey Berthelot's combining rule, T T T ij ii jj * * * ≈ ( ) . In the global phase diagram, this rule corresponds to a semicircular curve, which intersects, or comes close to, the domains of the phase diagram classes 1 P , 1 P l, 2 P , 2 P l, 1 C 1 Z , and their azeotropic subclasses: These classes are most likely to be observed experimentally.
For mixtures of molecules of unequal sizes the global λζ phase diagram loses its symmetry.
The triangular shaded area in Fig. 7, named shield region because of its shape, contains domains of rather complicated phase diagram classes containing quadruple points lllg.It has been argued that the shield region is fictitious and cannot be reached experimentally [25], but an experimental study of (water + n-alkane) mixtures [26] shows a  For the modelling of fluid phase equilibria it is advisable to look up the global phase diagrams of the model used-they are available for several popular equations of state.If the parameters of the model lie safely inside one of the major domains of the global phase diagram, one can be reasonable sure of the phase diagram class.But if the parameter set lies close to a domain boundary or-even worse-in one of the sensitive regions, caution is advised.
Independent of the availability of a global phase diagram it is good technique to compute the critical curves of a mixture (using an algorithm that finds all of them!), because this will hint the phase diagram class and indicate the existence of perhaps unexpected demixing phenomena.

Conclusion
In this article we have attempted to address some problems that can make the modelling of fluid phase equilibria with equations of state difficult: proper assessment of the quality and reliability of experimental data, difficulties stemming from wrong identification of fluid phases, and finally ensuring that all relevant phase equilibria have been found.The latter is a difficult problem indeed, and even today the total number of phase diagram classes for binary fluid mixtures is not really known, let alone that of ternary or higher mixtures.
On the other hand, even for ternary mixtures the situation is not hopeless: In many practical applications there is one volatile component together with a few chemically similar heavy components, so that the mixtures can be treated as pseudobinary at least as far as the phase diagram classes are concerned; an example is the removal of fatty acids from olive oil by extraction with supercritical carbon dioxide.Or there is a mixed extracting fluid and a single product.In either case, the study of binary phase diagrams is very useful.
In our examples we have mostly used rather simple equations of state.But as even such equations capture essential features of phase equilibria, we suppose that the advice given in this work is useful for work with more advanced thermodynamic models, too.
to the molar Gibbs energy, and (b) that the Gibbs energy of solid and vapor are equal at the sublimation pressure.Consequently, the chemical potential of a pure solid is

Fig. 1
Fig. 1 High-pressure phase diagram of pure anhydrous caffeine.○: experimental data (transitiometry); : triple points; curves: Clapeyron equation (C) or equations of state.The sublimation curves are not visible at this scale.

Fig. 5
Fig. 5 Barotropic inversion.Left side: far away from the critical point, yielding two separate phases with inverted positions; right side: close to a critical point, yielding a fog-like two-phase state.

Fig. 6
Fig.6 Entropic inversion in the (methane + propane) system at 250 K, calculated with the PR equation of state.
the d ij can be regarded as characteristic interaction energy densities of the mixture components.The T ij * and b ij are the energy and the covolume parameters of the CSRK equation,

Fig. 9
Fig.9 Class IV and its neighbour classes, together with the boundary states (middle column).For an explanation of the symbols and curves see Fig.6; in addition: •: tricritical point, : double critical endpoint, ■: zero-Kelvin endpoint.

Table 1
Molar volumes of the solid phases of anhydrous caffeine.