Documentation home

Introduction to phase diagrams

File structures, file names and line endings


Calculating phase diagrams


av PT Calculations

activity-composition model documentation

Mineral abbreviations

Input file documentation

Output file documentation

Internally consistent dataset documentation



page last updated: 28.11.2011

Calculating phase diagrams

This section outlines the basic instructions for calculating phase diagrams using thermocalc, It is in part based on the course CD notes and the information in:

Powell, R, Holland, TJB, & Worley, B, 1998. Calculating phase diagrams involving solid solutions via non-linear equations, with examples using THERMOCALC. Journal of Metamorphic Geology 16, 577–588.

Before reading this you should read the introduction to phase diagrams page

In addition to the information here, there are a number of tutorials aimed at different levels of experience in the tutorials page. What this page covers is the principles and logic behind calculating different types of phase diagrams, not the details of using thermocalc to do so. For that you should consult the tutorials page.

Introduction & background

The information below is from the course CD, written by Roger Powell

In the application of equilibrium thermodynamics to the calculation of phase diagrams, there are two approaches that can be followed: one based on the minimisation of Gibbs energy, the other being a derivative equivalent based on the solution of sets of non-linear equations. The main non-linear equations involved are the “equilibrium relationships”: the relationships for balanced chemical reactions between the end-members of phases that are in equilibrium with each other:

0=∆G0 +RT ln K                                                                    (1)     

In this, ∆G0 is the Gibbs energy of the reaction between the pure end-members in the same structure as the phases in which they occur, K is the equilibrium constant, in terms of the activities of the end-members in their phases, T is temperature, and R is the gas constant.

Thermocalc follows this non-linear equation approach. Part of the reason for this was that, at the point when this facility was added to thermocalc in the late ’80s, thermocalc already was using a multi-equilibrium-relationship approach to geothermometry/barometry (the average PT approach). Also it seemed to be the easiest way to draw, for example, PT projections.

The problem of calculating mineral equilibria can be thought of in terms of the number equations (constraints) involved and the number of unknowns being calculated. Average PT and phase diagram calculations may be contrasted, with the thing they have in common that each involves an independent set of reactions between the end-members of the phases involved in the equilibrium.

Phase diagram calculations

The independent set of reactions involves the end-members in the minerals of the equilibrium being calculated. In contrast to the average PT case, the system of equations formed by the equilibrium relationships is either exactly determined (ie all the unknowns can be calculated), or is under-determined, requiring that one or more variables need to be specified, depending on the variance of the system. For an n component model system with p phases, the well-known expression for the variance is np+2. Whereas for pseudosections the non-linear equations are augmented by composition ones, for P–T projections and compatibility diagrams just equilibrium relationships are involved. Consider first the latter, simpler case. If phase k involves ek end-members, then it involves ek − 1 composition variables. For the p phases in the equilibrium, there will be

equation 1

end-members of phases, and

equation 2

composition variables. The number of reactions between the end-members that make up an independent set is the number of end-members minus the number of components, sn (Powell & Holland, 1988). Given that there is a non-linear equation for each reaction in the independent set, these relationships indicate how many unknowns can be solved for, and therefore how many must be set, because the sn equations can only be solved for sn unknowns. The number of things that have to be set in order for an equilibrium to be calculated can be represented in terms of degrees of freedom, equal to the number of unknowns, (s − p)+2, minus the number of equations, s − n, giving n − p +2. So the number of degrees of freedom is just the variance. Setting unknowns may involve setting P and/or T, or setting compositional variables, as would be done, for example, in calculating composition isopleths on a P–T diagram. Therefore, by variance, v:

v = 0. For an invariant equilibrium, a point on a P–T diagram, all the unknowns can be solved for. With no composition variables set, the equilibrium involves n +2 phases; the P–T of the point, and the compositions of all the phases in the equilibrium, can be solved for. For c composition variables set, the equilibrium involves n +2 − c phases. For example, a divariant assemblage involving n phases, with 2 composition variables set, corresponding to the intersection of two isopleths, is (effectively) invariant, and the P–T and remaining composition variables of the equilibrium can be solved for.

v = 1. For a univariant equilibrium, a line on a P–T diagram, if one of the unknowns is set (i.e. one of P, T and the composition variables), then the remaining unknowns can be solved for. With no composition variables set, the equilibrium involves n +1 phases, and, given, say, P, the T and the compositions of all of the phases can be solved for. For c composition variables set, the equilibrium involves n +1 − c phases. For example, a divariant assemblage involving n phases with 1 composition variable set, is a line on a P–T diagram (an isopleth), the equilibrium is (e ectively) univariant and, given, say, P, the T and the remaining compositions of the phases can be solved for.

v = v. For a v-variant equilibrium, if v of the unknowns are set (i.e. v of P, T and the composition variables), then the remaining unknowns can be solved for. If there are no composition variables set, then v ≤ 2, or the equations cannot be solved. With no composition variables set, the equilibrium involves n +2 − v phases; for c composition variables set, the equilibrium involves n +2 − vc phases.


The calculation of P–T and T–x/P–x pseudosections with the non-linear equation approach involves augmenting the non-linear equations formed by the set of equilibrium relationships with a set of equations derived from mass balance constraints. For the specified bulk composition that is being used for the pseudosection, these additional equations relate the mole proportion of each component in the bulk composition with the sum of the calculated mineral compositions multiplied by their modal proportions. In other words, the bulk composition must be able to be made up of an assemblage of the phases of interest. Of course the phases in the assemblage must each have non-negative modal proportions.

For the n-component system considered above, there will be a mass balance equation for each component, giving n equations additional to the equilibrium relationships involved. There are also p additional variables, the modes of the phases. Then the total number of equations is sn + n, i.e. s equations, in sp + p + 2, i.e. s + 2 unknowns. So, if two things are specified, the number of equations equals the number of unknowns regardless of the number of phases involved. This means that the compositions of the minerals and their modes can be calculated for an equilibrium of any variance at given PT, once a bulk composition is specified (see also Spear, 1986). Alternatively if two modes are specified, the compositions of the minerals, the remaining modes and the PT can be calculated. Further, if the modes are specified to be zero, then the PT corresponds to a point where 4 boundary lines come together on a pseudosection (for the variance of the highest variance field at the field being greater than 1).

P-T projections


P-T projections, also known as P-T grids or petrogenetic grids show only invariant and univariant equilibria. Because P-T projections don’t require the input of bulk rock composition (ie are independent of rock composition) their topology is only dependent on the thermodynamics of the phases used. Thus, barring changes in the internally consistent dataset and the a-x models used, P-T projections only rarely need to be calculated and those appropriate for many chemical systems already exist in various publications. Furthermore, the number of chemical systems for which a P-T projection  or a useful P-T projection can be drawn is typically limited to the simpler systems. In larger systems there may not exist univariant or invariant equilibria, or there may be only a scattering of short univariant equilibria stable.  In other large systems, there may be so many univariant and invariant equilibria stable as to render the P-T projection so complex as to be unreadable. However, this can sometimes be alleviated by considering a large number of phases to be in excess but with the corollary that this process reduces the composition space for which the diagram is appropriate.

Calculating P-T projections is usually fairly straightforward but there are places where one can run into trouble (see also the introduction to phase diagrams page). For example, some reactions contain singularities characterised by compositional coplanarities with a phase swapping sides on the reaction either side of the singularity. If one such singularity occurs close to an invariant point it can make thing complex, and chaotic if you miss it. Singularities can also result in two univariants that emanate from an invariant point recrossing each other in an indifferent manner (in this case, one or more of the common phases to both reactions will have different compositions at the P-T conditions of the indifferent crossing).

Typically, the only foolproof way to build a P-T projection is from all the subsystems up. Obviously this is a pretty laborious process if you want to ultimately draw a grid in a large system. Even in the very small systems like MAS, you still have to get the subsystem grid topology right, recognising there are always two possible topologies. Fortunately, there are experiments in these smaller sytems that show the stable assemblages. In slightly larger systems such as KFMASH there are also experiments and information from natural assemblages that can help you decide which equilibria are likely to be stable. For example two alternative topologies around an invariant point are shown below


These two P-T projections are two equally valid (in terms of Schreinemakers rules) alternative topologies and are inversions of each other. Now, the one on the left is the correct grid based on several lines of logic.

Firstly the reactions and possible divariant reactions either side of each reaction are consistent with what we observe in rocks, especially metamorphic belts where we see a metamorphic field gradient (e.g. g-chl assemblages occur at lower grades than st-bi assemblages)

Secondly we could look at the full reactions when we calculate them (ie including the in-excess phases) and see that H2O is on the high-T side of each reaction in the P-T projection on the left and is on the low-T side in the P-T projection on the right. As metamorphic reactions are typically dehydration reactions (higher T rocks typically have less water) we can infer that the P-T projection on the left is the correct one. The same logic can be applied to which side melt occurs in high-T, P-T projections.

Calculating P-T projections

This basically involves calculating univariant and invariant equilibria and following Schreinemakers rules. If you are building the P-T projection up from subsystems, the invariant points in the subsystems with one less component than the full system (ie KFASH or KMASH to KFMASH) will all be full system univariants that will form full system invariants where they cross.

For calculating these invariants and the univariant lines that emanate from them, first calculate the invariant point of interest. To calculate the univariant lines, a good trick is to type in the invariant assemblage and then set the variance to 1 at the prompt. Here thermocalc will calculate each line sequentially by automatically omitting one phase. It is best to work in a pretty small P-T window (e.g. 0.5 kbar by 30°C or even smaller) At this stage you can plot the lines up in drawpd not worrying about the invariant point and the relative stability of lines across it.

Before trying to sort out the stable and metastable parts of the reaction you should look carefully at the output checking for singularities, the reaction you want is that at the invariant point.  Doing Schreinemakers analysis from the two crossing reactions is normally fairly straightforward as the stable extensions are those that lead back to the subsystem invariants. If you are starting without the subsystem information you will need to decide that one reaction’s extension is stable and which way around the reaction is.

What I normally do is sketch the bundle with a pencil and rule, label it with the reactions and the ‘out’ phases. I then partially erase the metastable extensions so they are dashed then draw it properly in drawpd by calculating each line over a larger P-T range.

Doing all this, you will eventually come across some reactions that sit very close to each other and are hard to visually separate. Then you need to look carefully at the output. Thermocalc may also present degenerate reactions in a way that is not immideiately obvious. For example in KFMASH at high-T at say a [bi] invariant point there are only two other K2O-bearing phases (ksp & liq). What you will find is that the [ksp] and [liq] reactions are the same reaction and actually in FMASH. This reaction will pass through the invariant point unhindered.

As you build up the P-T projection, you will need to check each intersection between lines and calculate each invariant and univariant bundle the results. Keep in mind that any intersection that results in apparent negative variance is an indifferent crossing where the compositional space of the two reactions involved do not intersect.

A tutorial for constructing P-T projections that outlines in more detail the use of thermocalc to do this will eventually (I promise…really) appear on the tutorials page.

P-T, T-x & P-x pseudosections


Pseudosections are multivariant phase diagrams that show the mineral equilibria that are ‘seen’ by a given bulk rock composition or compositional vector. For a more detailed definition and explantion see the introduction to phase diagrams page. Pseudosections are typically composed of a series of multivariant fields, but may also commonly show segments of univariant equilibria ‘seen’ by a particular rock composition. In constructing these, the main calculations are the boundaries where a particular assemblage field ends via the gain or loss of a phase and points where such boundaries intersect. In addition pseudosections can be contoured for mineral mole proportions (modes) and mineral composition (composition isopleths).

Below is a section on calculating pseudosections from the course CD written by Roger Powell. Again this deals with the logic involved in constructing such diagrams and not the mechanics of running thermocalc in detail. There is a basic level tutorial for pseudosection calculations in the tutorials page 


Calculating P-T, T-x & P-x pseudosections

Pseudosections involve the invariant and univariant equilibria they inherit from PT projections, as well as additional boundary lines and points. In a PT pseudosection these inherited equilibria are just those parts of the PT projection “seen” by the bulk composition being considered. In a Tx pseudosection these inherited equilibria occur as horizontal lines, spanning the range of bulk compositions that “sees” an equilibrium. The boundary lines and points separate fields of different variance; across lines the variance changes by one, and through points by two. The boundary points and lines can be most easily discussed with a new notation. Labelling is done in terms of the lower variance assemblage involved, with the names of phases with zero modes given in brackets. So for example in Fig. 1 (iii), going from the g–chl–bi divariant field, through the point at 8 kbar and 580°C, into the chl quadrivariant field, the modes of bi and g both go to zero at the point. The point is therefore labelled “chl (bi g)”. Similarly, in going down T at 10 kbar from the g–chl–bi divariant field into the g–chl trivariant field into the chl quadrivariant field, the two lines crossed are labelled “chl g (bi)” and “chl (g)”.

The labelling of boundaries (and points) in terms of the phases whose mode(s) go to zero is shown in Figs 1 and 2. In thermocalc terminology such a boundary line is called “effectively” univariant (regardless of the variance of the fields) because, in terms of variables, the boundary is a line. The effective univariance comes from the the combining of the bulk composition with the equlibrium relationships.



In terms of a coherent view of pseudosection construction features, different rules appear to apply across univariants, with this being the one situation in which variance does not change by one (or, alternatively, the number of phase in the assemblage changing by one) across a “boundary”. In the (proper) univariant case, divariant fields occur across the univariant, with the swapping of (just) one phase in the assemblage. However, the univariant case can be made to comply with all others by representing the univariant as an (infinitely) thin field involving the full univariant assemblage. Then the standard way of thinking about boundaries in pseudosections applies to univariants as well.




Points in pseudosections are where two modes are zero with respect to the lowest variance assemblage at the point. At a point there is one n+1, two n, and one n− 1 variance fields. The two n− 1 to n variance field boundaries involve one each of the phases whose modes are zero at the point. In terms of calculation, these two boundaries terminate at the point, whereas the other two boundaries (n to n+1 variance field boundaries) continue through the point. Considering a compatibility diagram, like the portion of AFM in (a) in Fig 3, and considering that the phase relations move with respect to the bulk composition (star) with changing PT, the variance relationships around points in PT pseudosections should be clear(er), as well as accounting for which lines will have metastable extensions.

Getting started in drawing a pseudosection is particularly straightforward if a univariant that is known to be stable (from the PT projection) is “seen” by the bulk composition (see below). Otherwise it is a matter of inspired guesswork,  building a T-x or P-x pseudosection from an existing diagram to your new composition or Gibbs energy minimisation (via dogmin) can be used as a last resort.


Figure 3.

Calculating along a line, 3 things can happen:

  1. The line is cut off by a lower variance field, with a new phase joining the equilibrium. If this line is a boundary between an n and an n +1 variant fields, then at the point where the line is cut off, an n − 1, two n and an n +1 variant fields meet. Case (b) in Fig 3.

  2. The line dies, with a mode going to zero. If this line is a boundary between an n and an n +1 variant fields, then at the point where the line is cut off, an n, two n +1 and an n +2 variant fields meet. Case (c) in the Fig 3.

  3. thermocalc fails to calculate the line beyond a certain point, but no mode is going to zero. This signals a failure of the non-linear equation solver in thermocalc. The usual resolution is to vary the starting guesses in the datafile, using ones that correspond to the values along the univariant that thermocalc could calculate.

Calculating a point, 2 things can happen:

  1. the point is calculated successfully

  2. thermocalc fails to calculate the point. This can arise for two reasons: 1. The qualitative assessment of the pseudosection geometry is flawed, and the point is not there to be found. This can easily happen locating where trivariants touch univariants, when a narrow field has been missed, for example or where two phases were approaching zero along a line and you have picked the wrong one.  2. A failure of the non-linear equation solver in thermocalc is signalled. As before, the usual resolution is to vary the starting guesses

Compatibility diagrams


Compatibility diagrams are a very useful way of visualising the mineral equilibria that occur for different rock compositions at a given P-T. While such diagrams typically use pure oxides as the apices (e.g. the classic AFM compatibility diagram), this need not be the case. The key point in producing such diagrams is effectively reducing the composition space of the system to a ternary.

Calculating Compatibility diagrams

The constructional features of compatibility diagrams are illustrated in figure below. The tie triangles (divariant equilibria) in the full system (within the triangle) as well as the divariant equilibria in the subsytems (edges of the triangle) specify much of the geometry, including the apices of the one-phase fields (quadrivariant equilibria). The edges of the one phase fields can be determined via the adjacent trivariant equilibria. In cases such as this, in which Fe–Mg is the dominant substitution, the edges are close to being straight; in other systems, such edges can be strongly curved.

figure 4

Although the most commonly used compatibility diagrams use pure oxides as the apices, it is often useful to use combinations of oxides, such as endmember compositions. This is particularly useful for reducing the composition space you wish to show. For example, to highlight the phase relationships for Fe-rich metapelites you could use the axes K2O-FeO-Fe0.75Mg0.25O as shown below.

figure 5

The use of a “endmember” as an axis is shown below for metabasic rocks. Here the axes are CaAl2O4–Ca2/7Fe5/7O–Ca2/7Mg5/7O, with axes at the base being the endmember compositions of ferroactinolite and tremolite (minus the in-excess phase components). The diagram below is also a good example of where phases can plot at infinity. Here the clinopyroxenes (omphacite & diopside) plot at infinity below the diagram and the tie lines involving these phases are simply drawn as arrows.

figure 6