The best way to find out more about running THERMOCALC is to work through the annotated THERMOCALC logfiles, and actually running them yourself to make sure you can see what is going on. As it says in the getting started section, running THERMOCALC is all about
once you have the phase definitions set up in the datafile and the preferences file is configured to your liking.
For a discussion of further aspects of generating phase diagrams using THERMOCALC, look at run documentation here, in combination with the annotated THERMOCALC logfiles.
Running THERMOCALC to calculate phase diagrams involves a number of steps:
Choose a model system in which to do the calculations. A model system is just the chemical system, usually specified in terms of oxides, in which the equilibria to be calculated can be represented. This then specifies which phases, and the substitutions within them, that will be able to be considered. Obviously these phases can only involve those end-members that occur in the thermodynamic dataset (or linear combinations of those end-members). For example, the model system normally used to consider metapelites is K2O-FeO-MgO-Al2O3-SiO2-H2O (or KFMASH) (Thompson, 1957), allowing most of the critical minerals, and the FeMg(-1) and (Fe,Mg)SiAl(-1)Al(-1) (Tschermak's) substitutions in them, to be considered.
In the following sections, there is information on calculating the main types of phase diagrams with THERMOCALC.
Projections are phase diagrams which show all of the invariant and univariant equilibria for all of the bulk compositions in a model system. They are the familar petrogenetic grids of the literature:
THERMOCALC output takes the form, in this case for the g+chl=st+bi univariant, the staurolite isograd reaction
P(kbar) T(°C) x(st) x(g) x(chl) y(chl) x(bi) y(bi) x(mu) y(mu)
4.0 540 0.977 0.965 0.792 0.608 0.819 0.448 0.915 0.879
21g + 31chl + 71mu = 10st + 71bi + 117q + 104H2O
6.0 572 0.937 0.911 0.604 0.593 0.638 0.417 0.812 0.860
19g + 32chl + 71mu = 10st + 72bi + 115q + 112H2O
8.0 599 0.884 0.846 0.466 0.579 0.496 0.386 0.713 0.841
17g + 34chl + 72mu = 10st + 73bi + 114q + 118H2O
10.0 622 0.818 0.770 0.361 0.565 0.383 0.355 0.619 0.821
15g + 35chl + 73mu = 10st + 73bi + 113q + 123H2O
giving not only the PT line for the reaction, but the way the phase compositions, and the balanced chemical reaction, vary along the line.
These are the simplest diagrams to generate in that it is very easy to calculate all the invariants and univariants for a system. However, THERMOCALC does not determine the stability of equilibria. When all the univariants, for example,in a system are calculated, many or most will be metastable. The stability of the equilibria in a system can be determined by applying Schreinemakers Rule,that the metastable extension of the reaction not involving phase i (i-out or [i]) lies between i-producing reactions (eg Zen, 1966), to the univariants around each invariant, then putting the invariants together. This process might be started in subsystems, as stable full system univariants emanate from stable subsystem invariants, eg the KFMASH univariant ctd + st = ky + g emanates from the KFASH invariant [ bi,cd,chl] at high P in the above Figure. Thus, it is useful to determine the grids for KMASH and KFASH before determining the KFMASH grid.
The approach of determining the stable univariants by Schreinemaker's Rule can sometimes be short-circuited by the method illustrated in the following Table
no T(°C) inc reaction 1 538 * ctd + ky = st + chl 2 555 * ctd = st + g + chl 3 557 1 ctd = st + bi + ky 4 559 1 ctd = st + g + bi 5 563 1 ctd + chl = st + bi 6 565 1 ctd = st + cd + bi 7 572 * g + chl = st + bi 8 572 1 ctd = g + chl + ky 9 575 1 ctd = g + bi + ky 10 581 1 ctd + ky = st + cd 11 581 1 ctd + chl = bi + ky 12 583 2 g + chl = bi + ky 13 585 * st + chl = bi + ky 14 588 1 ctd = st + cd + g 15 591 1 ctd = cd + g + ky 16 603 3 chl = st + cd + bi 17 604 1 ctd + chl = cd + g 18 604 1 ctd + bi = cd + g 19 604 1 chl = ctd + cd + bi 20 604 2 chl = cd + g + bi 21 605 1 g + chl = ctd + bi 22 619 1 st + ctd = g + ky 23 626 3 chl + ky = st + cd 24 627 2 chl + ky = cd + g 25 630 3 st + chl = cd + g 26 644 1 chl + ky = ctd + cd 27 649 2 st + chl = g + ky 28 660 5 st + cd = g + ky 29 670 * st + bi = g + ky 30 734 4 st + bi = cd + g
Table 1. All the reactions at 6 kbars in KFMASH, ordered in T (except for 4 reactions that were "not in range"). The bulleted reactions are the stable ones (see text).
Here all of the univariants at a particular P are ordered in T, with each reaction written low T \rightarrow high T (as given by THERMOCALC). Noting that the reaction i = j + k + l, means that i cannot be stable at higher T, so that any reaction involving i must necessarily be metastable, then reaction 2 in Table 1 excludes reactions with 1 in the "inc" column. Noting that i + j = k + l means that i + j cannot be stable at higher T, so that any reaction involving i + j (other than with i + j alone being on the high T side of the reaction) must be metastable, then reaction 7 excludes reactions with 2 in the "inc" column, 13 excludes reactions with 3 in the "inc" column, and 29 excludes reactions with 4 in the "inc" column. The same lines of logic apply with decreasing T, so 29 also excludes reactions with 5 in the "inc" column. Thus the stable reactions in the above Table are the ones marked with a bullet (*). This method is developed here for an effective ternary system but can be generalised to apply to any size system. Note that, if certain equilibria are not found by THERMOCALC, for example for starting guess reasons, then the resulting stability relationships will be incomplete. In pathological situations (for example systems with singularities, Worley et al., 1997) disconnected i + jtrivariants may be stable, making the method difficult to apply.
Having determined the P-T projection for a system, with the establishment of the stability relationships that this entails, makes the phase diagrams of the following sections easier to calculate and draw.
Compatibility diagrams show the mineral assemblages, and ranges of mineral solid solutions, at a specified P-T, for all of the bulk compositions in a model system. To be useful compatibility diagrams need to involve a maximum of 3 apices (or be the equivalent orthogonal diagram), ie be an effective ternary, so that it can be represented in 2D. For this to be true usually involves a reduction in the effective size of the system by stipulating that certain phases are present in all of the equilibria being considered. The reduction of KFMASH to an effective ternary (at conditions below the upper amphibolite facies, where muscovite + quartz breaks down) is normally done by considering that muscovite, quartz and H2O are "in excess" (Thompson, 1957), and the phase relationships are projected from these phases onto Al2O3--FeO--MgO (AFM). When these projections are made, the phases must be projected from the calculated compositions of the phases involved in each equilibria (eg Thompson, 1979).
One important practical aspect of compatibility diagrams, beyond the (geological) choice of which phases to project from, is the choice of projection plane, the plane onto which the phase relationships are to be projected. It is quite reasonable for phases to plot at infinity, as does K-feldspar in AFM, but if this can be prevented it is helpful, particularly if more than one phase plots there, or if that phase is a solid solution. More seriously, the projection plane should be chosen to be between the projecting phases and the "nearest" solid solutions (Powell and Sandiford, 1987; Worley and Powell, 1997). See Worley and Powell (1997) for a discussion of this problem in relation to choosing a projection plane for NCKFMASH. In fact, the classic AFM projection for KFMASH fails at lower temperatures as a consequence of the projection plane being further away from muscovite than is biotite, for example in KAlO2--(FeO,MgO)--Al2O3. At lower temperature, on the AFM compatibility diagram, biotite in equilibrium with alkali feldspar plots at more Al-rich compositions than biotite in equilibrium with chlorite, for Fe-rich and Mg-rich compositions.
THERMOCALC is set up to calculate compatibility diagram information in a straightforward way. Having decided which phases to project from (if it is necessary), and what the compositions are of the apices of the projection plane, then the information is added to the datafile in the form of scripts. Either these scripts can be entered directly, or THERMOCALC can be run with the script "project yes" in the datafile. In this case the projection information is prompted for, and the required script information is placed in the logfile for copying into the datafile. The scripts to instruct THERMOCALC to provide AFM coordinates for the phases in each equilibrium it calculates is
fluidexcess yes
setexcess mu q
project yes
% [script to enter the projection plane info]
% -------------------- Al2O3 MgO FeO K2O
setproject A 1 0 0 0
setproject F 0 0 1 0
setproject yes M 0 1 0 0
% -------------------------------------------
The necessary calculations usually involve all of the divariants being calculated for the full system, giving the information for the interior of the compatibility diagram, as well as all of the divariants for the sub-systems, giving the information for the edges of the diagram. Establishing which divariants are stable can be done most easily via the relationships in the P-T projection.
THERMOCALC output takes the form, in this case for the g+chl+bi divariant
<==================================================>
phases : g, chl, bi, (mu, q, fluid) or [st,ky]
--------------------------------------------------------------------
P(kbar) T(°C) x(g) x(chl) y(chl) x(bi) y(bi) x(mu) y(mu)
6.0 560.0 0.939 0.690 0.571 0.720 0.381 0.864 0.837
proj A F M mu phase
g 0.250 0.704 0.046 -0.250
chl 0.190 0.559 0.251 -0.167
bi -0.228 0.871 0.356 0.500 -0.500
<==================================================>
with the AFM coordinates of the phases given in projection from the calculated composition of the mu.
In the AFM examples, the boundaries of the one phase fields are close to being straight lines between the apices of the tie triangles, so these boundaries are drawn as such above. In more complex systems, these boundaries can be quite curved. The easiest way to trace such trivariant-quadrivariant boundaries in compatibility diagrams is via the trivariant equilibrium and set isopleths, looking at an appropriate range of isopleths for one of the composition variables for the phase of the one phase field. For example, to trace the boundary of the chlorite one phase field on the low Al side of the above compatibility diagram, calculations on the trivariant, chl-bi(-mu-q-H2O), with x(chl) set to say 0 through to 0.75 at 0.05 intervals, will give the coordinates of the ends of the tie lines between biotite and chlorite. As well as defining the low Al side of the chlorite one phase field, this will give the high Al side of the biotite one phase field.
In complex systems, for example the NCKFMASH system considered by Worley et al. (1997), the reduction of the system so that it is effective ternary requires considering five phases to be in excess (eg plagioclase, biotite, muscovite, quartz and H2O). There comes a point when such reductions cannot be done in a useful way. Further, such reductions may conceal more than they reveal if projected bulk compositions migrate too much with changing conditions. This can happen if most of the mode of the bulk composition is being projected from, particularly if projection from several solid solutions is involved. In this case, it is likely that an x-x pseudosection, showing the phase relationships for a bulk composition plane through the system, would be more useful. Such diagrams can be calculated with THERMOCALC, but not in a straightforward way.
P-T pseudosections
A P-T pseudosection is a P-T diagram on which just the phase relationships that apply to (or are "seen" by) a specified bulk composition are shown. In the same way that a compatibility diagram can have one phase (quadrivariant), two phase (trivariant), and three phase (divariant) fields, so can P-T pseudosections. In addition, P-T pseudosections can involve those parts of a P-T projection which are seen by the bulk composition: univariant lines and invariant points.
To access pseudosection facilities in THERMOCALC, the "pseudosection" script in the datafile must be set. Then, running THERMOCALC, the user is prompted for the bulk composition that is to be used. A bulk composition could, for example, be derived from a whole rock analysis, or be an estimated equilibration volume composition constructed from the analysed compositions of phases in an assemblage. This bulk composition appears in the logfile in an appropriate form, and can be copied across to the datafile. The script to instruct THERMOCALC to provide pseudosection information with the bulk composition used in the examples is
pseudosection yes
% [script to enter the bulk composition info]
% ----------- Al2O3 MgO FeO K2O ---------
setbulk yes 41.89 18.19 27.29 12.63
% ------------------------------------------------
setmodeiso yes
zeromodeiso yes
The latter two scripts access a facility that is central to calculating pseudosections. Various points and lines on pseudosections can be understood in terms of particular modes going to zero (so "zeromodeiso"). For example, in the above Figure, the lower P extent of the univariant is where the modes of g and st go to zero. The notation used for a such a point is chl bi (st g), indicating that the equilibrium involves all four phases, but with the modes of the phases in brackets set to zero. It can be calculated, using the above script, by specifying the univariant equilibrium and specifying that the modes of g and st are set to zero at the appropriate prompt in THERMOCALC. All the points and lines in the above Figure are calculated in this way (apart from the univariant lines and invariant points). Note that on such pseudosections, variance changes by 2 between the highest and lowest variance fields at points.
THERMOCALC output takes the form, in this case for the bi chl (st g) point that terminates the seen part of the g+chl=st+bi univariant as it goes into the bi+chl trivariant
Al2O3 MgO FeO K2O
41.89 18.19 27.29 12.63
<==================================================>
phases : st, g, chl, bi, (mu, q, fluid)
P(kbar) T(°C) x(st) x(g) x(chl) y(chl) x(bi) y(bi) x(mu) y(mu)
6.31 576.8 0.9292 0.9014 0.5797 0.5906 0.6134 0.4120 0.7962 0.8573
mode st g chl bi mu
0 0 0.4439 0.1019 0.4543
<==================================================>
with the PT coordinate of the point given, as well as the modes of the phases that apply at the point.
The usual procedure followed when constructing a P-T pseudosection is to start by plotting the chosen bulk composition on compatibility diagrams for a series of P-T conditions corresponding to those of the pseudosection. This indicates which mineral assemblages, ie which fields, should be expected. Then the stable univariants are examined; where a univariant is seen, the modes of the divariants either side of the reaction are given. Along the reaction, changes in the divariant seen indicate the crossing of a trivariant termination on the univariant. Such terminations can then be calculated directly. Then the resulting divariant-trivariant boundariescan be followed by calculating the divariant with the mode of the appropriate phase set to zero. Looking at chl bi (g), the mode of bi progressively decreases up P, eventually going to zero: this is chl (g bi), the corner where the chl quadrivariant terminates against the divariant. The corner is calculated by specifying the divariant and setting the modes of g and bi to zero. So P-T pseudosections are built up, each point and line at a time, once the parts of the stable univariants that are seen have been determined.
Note that divariants, in particular, can be very narrow. As a consequence care is needed, meaning that a very fine P or T step may be needed in looking at univariants, to ensure that narrow divariants are not missed. Once identified, the width of a divariant does not affect the calculation of its boundaries, as they are specified separately in terms of the modes that go to zero there, in going into the adjacent trivariants.
Although the stability relationships can usually be established directly from those in the PT projection---specifically in providing a starting point from which a pseudosection can grow---sometimes, particularly if no part of a stable univariant is seen by the bulk composition, an alternative is needed. This is provided by a Gibbs energy minimisation facility in which all the possible assemblages seen by the bulk composition are calculated. The stable one is the one of lowest Gibbs energy. Such output is:
P(kbar) T(°C) st g chl bi mu ky q H2O var
9.0 565 - - * - * - * * 4
9.0 570 - - * - * - * * 4
9.0 575 - - * - * - * * 4
9.0 580 - * * - * - * * 3
9.0 585 - * * - * - * * 3
9.0 590 - * * - * - * * 3
9.0 595 - * * * * - * * 2
9.0 600 - * * * * - * * 2
9.0 605 - * * * * - * * 2
9.0 610 - * * * * - * * 2
9.0 615 * * - * * - * * 2
9.0 620 * * - * * - * * 2
Table 2. Gibbs energy minimisation by looking at the Gibbs energies of all divariant, trivariant and quadrivariant equilibria at 9 kbar and 565 through to 620^\circC for the bulk composition: mol% Al2O3 = 41.89, mol% MgO =18.19, mol% FeO = 27.29, and mol% K2O = 12.63 calculated by THERMOCALC.
A f(P,T)-x pseudosection is a diagram on which the phase relationships that apply to (or are "seen" by) a specified bulk composition line along a P-T line (ie f(P,T)) through the model system are shown. Normally such a pseudosection is drawn at constant P, making it a T-x pseudosection, or at constant T, making it a P-x pseudosection. In the following, for simplicity, they will be described in terms of a T-x pseudosection as the logic and approach are equivalent for both. For example:
In the same way that a P-T pseudosection can have one phase (quadrivariant), two phase (trivariant), and three phase (divariant) fields, so can T-x pseudosections. On T-x pseudosections univariant (and invariant) equilibria occur as horizontal lines whose length is controlled by the range of bulk compositions along the composition line that "see" the equilibria.
To access these pseudosection facilities in THERMOCALC, the "pseudosection" script in the datafile must be set. Then, running THERMOCALC, the user is prompted for the ends of a bulk composition line that is to be used. Alternatively, a composition line can be generated from a single bulk composition used, for example, in a P-T pseudosection, by for example varying Fe--Mg through that composition. For example, using the bulk composition from the previous section as a basis for exploration of Fe-Mg variation, the scripts to instruct THERMOCALC to provide pseudosection information are
pseudosection yes
% [script to enter the bulk composition info]
% ----------- Al2O3 MgO FeO K2O ----------
setbulk yes 41.89 45.48 0 12.63
setbulk yes 41.89 0 45.48 12.63
% ------------------------------------------------
setmodeiso yes
zeromodeiso yes
In the previous bulk composition, MgO = 18.19 and FeO = 27.29; adding these gives 45.48, so the bulk composition line extends from MgO = 45.48 and FeO = 0, to MgO = 0 and FeO = 45.48.
The same approach is followed in constructing a T-x pseudosection as in constructing a P-T pseudosection, with the main difference being that any calculation is repeated for each bulk composition. As before, the stable univariants are examined first, looking to see the composition range, if any, over which a univariant is seen. As in the P-T case, the modes of the divariants either side of the reaction are given for each bulk composition seen. The T-x pseudosection is built up in the same way as for the P-T pseudosection by following divariant-trivariant lines, recognising the existing of qudrivariant fields, following the trivariant-quadrivariant boundaries, and so on. One difference arises in the calculation of points in that THERMOCALC calculates the trace of such a point with varying P, and the user must interpolate to get the T and x values for the fixed P of the pseudosection.
Again, note that fields can be very narrow. As a consequence, small sections of the bulk composition line may need to be examined in looking at univariants, to ensure that narrow divariants are not missed.
or go back to