GeolCarp_Vol57_No2_131

GEOLOGICA CARPATHICA, APRIL 2006, 57, 2, 131—142

www.geologicacarpathica.sk

Saturnin, R language script for application of accessory-

mineral saturation models in igneous geochemistry

VOJTĚCH JANOUŠEK

1, 2

Division of Mineralogy and Material Science, Hellbrunnerstraße 34, A-5020 Salzburg, Austria

2Present address:

Czech Geological Survey, Klárov 3/131, 118 21 Prague 1, Czech Republic; janousek@cgu.cz,

Phone: +420 251 085 308; Fax: +420 251 818 748

(Manuscript received January 11, 2005; accepted in revised form June 16, 2005)

Abstract: Accessory minerals play a crucial role in the petrogenesis of granitic rocks. They control the response of
isotopic systems (U—Pb, Lu—Hf, Sm—Nd) and rule the geochemical variation of many important trace and minor
elements during partial melting and fractional crystallization. The experimental effort in the past thirty years has
resulted in the formulation of saturation models with manifold applications in igneous geochemistry. Numerical
simulations show that the zircon saturation thermometer is rather insensitive to analytical errors, the presence of a
minor inherited component, as well as moderate accumulation of feldspars and some ferromagnesian minerals (Ca-
poor pyroxenes, less so olivine). However addition of even minor cumulus Ca-pyroxene or Ca-amphibole would
quickly render the calculated temperatures too low. Apatite saturation thermometry is a poor tool for felsic metaluminous
rocks, being oversensitive to errors in the phosphorus determination as well as the presence of extraneous apatite.
Even stronger element of uncertainty is added for peraluminous lithologies, when the increased apatite solubility is
inadequately accounted for by the current models and where other phosphorus-bearing minerals, most importantly
feldspars and monazite, come into play. The newly developed software Saturnin (http://www.gla.ac.uk/gcdkit/saturnin)
written in the freeware R language performs otherwise tedious calculations of zircon, monazite and apatite saturation
in igneous rocks. While the Windows users would probably find it easier to access the programme as a part of the larger
package GCDkit (http://www.gla.ac.uk/gcdkit), on Macintosh and Linux it can be used as a stand-alone application.

Key words: software, igneous geochemistry, fractional crystallization, partial melting, geothermometry, zircon,
apatite, monazite.

Introduction

Recent studies confirm the crucial role played by accesso-
ry minerals in the petrogenesis of granitoid rocks. These
minerals, for the most part, control the response of isotopic
systems (U—Pb, Lu—Hf, Sm—Nd) and rule the geochemical
variation of many important trace and minor elements (e.g.
Zr, P, REE, U, Th) during partial melting and subsequent
fractionation of the magma (Gromet & Silver 1983; Miller
& Mittlefehldt 1984; Sawka 1988; Watt & Harley 1993;
Bea 1996; Ayres & Harris 1997; Hoskin et al. 2000).

In particular three minerals, vital in geochronology and

isotopic geochemistry, have been the subject of copious
studies of natural and experimental systems: zircon (Wat-
son & Harrison 1983; Harrison & Watson 1983; Watson
1996; Baker et al. 2002; Miller et al. 2003; Hanchar &
Watson 2003; Cherniak & Watson 2003), apatite (Watson
& Capobianco 1981; Harrison & Watson 1984; London
1992; Bea et al. 1992; Pichavant et al. 1992; Wolf & Lon-
don 1994, 1995; Cherniak 2000) and monazite (Rapp &
Watson 1986; Montel 1993; Wolf & London 1995; Cher-
niak et al. 2004).

The experimental effort of these as well as numerous

other authors resulted in formulation of saturation models
with manifold applications in igneous geochemistry. The
current paper provides a critical overview of those for zir-
con, apatite and monazite, a summary of the relevant for-

mulae and their practical implementation in the freeware
R language.

Applications of saturation models in igneous

geochemistry

Fractional crystallization

As the fractionation of accessory minerals has profound

effects on the trace-element budget of the parental magma,
a good understanding of the factors controlling their crys-
tallization is a fundamental issue in igneous petrogenesis
(Evans & Hanson 1993; Bea 1996; Hoskin et al. 2000).
The growth of phases like zircon, apatite or monazite
seems to be governed by geochemical species concentrat-
ed in the accessory but incompatible in the main rock-
forming minerals. These elements, termed essential
structural constituents (ESC) by Hanson & Langmuir
(1978), are Zr in zircon, P in apatite, P and LREE (Ce) in
monazite. At the onset of magma crystallization, such a
trace element would behave incompatibly, gradually in-
creasing its concentration in the magma until a required
saturation threshold is reached (point S in Fig. 1a). For this
segment of the crystallization path (OM—S), the proportion
of the remaining melt in each magma batch has to be re-
duced significantly by crystallization of the major rock-

133

SATURNIN R LANGUAGE FOR APPLICATION OF SATURATION MODELS IN IGNEOUS GEOCHEMISTRY

forming minerals for the magma to become saturated. As a
consequence, the growth of the accessory crystals is possi-
ble only at a temperature significantly below the liquidus.

The crystallization of the accessory phase at tempera-

tures close to liquidus may commence at the point S, when
the element starts acting as an essential structural constitu-
ent. Its concentration is from this point on buffered by the
accessory growth at a saturation level (Hanson & Lang-
muir 1978; Evans & Hanson 1993). However this level is
unlikely to remain constant because of the dropping tem-
perature as well as changes in bulk composition owing to
the fractional crystallization. The changing bulk magma
composition affects some parameters featuring in the satu-
ration equations (e.g. M parameter for zircon, see below).
In any case the variation in the ESC is likely to be much
more limited than for strongly compatible as well as in-
compatible trace elements (Evans & Hanson 1993).

Any extraneous component in the accessory mineral

population (inherited from the source or forming xenocrysts
captured from the early assimilated country rocks) would
be likely to dissolve at the undersaturated segment of the
path  (OM—S) but would tend to be preserved subsequently
(S—DM)  (Harrison & Watson 1983; Watson & Harrison
1984; Miller et al.  2003).

In reality, the growth of an accessory may deviate con-

siderably from the ideal behaviour, the crystallization be-
ing often interrupted by periods of sudden resorption. The
reason can be fluctuations in temperature and bulk chem-
istry of the magma (due to the assimilation/magma mix-
ing) as well as effects of the so-called localized saturation.
As proposed by Bacon (1989), non-equilibrium concen-
tration gradients may form adjacent to rapidly growing
rock-forming minerals, permitting local crystallization of
accessory minerals, even when the bulk chemistry of the
magma would rule this out theoretically. However, if the
temperature can be assumed to have been approximately
constant, and effects of sudden changes in bulk chemistry
and localized saturation can be neglected, the saturation
models could serve several purposes.

Estimating the liquidus magma temperature

In order to approximate the temperature of magma from

which the given accessory mineral crystallized, several as-
sumptions have to be fulfilled:

i. The population of the accessory has no extraneous

component (inheritance or xenocrysts), as otherwise the
obtained temperature would be overestimated. This effect
is the most severe for zircon in rather cold, felsic granite
suites (Chappell et al. 1998; Miller et al. 2003), many of
which have a S-type chemistry (Clemens 2003). The pres-
ence of zircon inheritance can be checked using imaging
techniques, such as cathodoluminescence or BSE micros-
copy (Paterson et al. 1992; Hanchar & Miller 1993; Han-
char & Rudnick 1995; Poller 2000; Nasdala et al. 2003).
Its impact on saturation calculations can be minimized by
corrections employing image analysis or mass balance
based on the degree of U—Pb discordance (Broska et al.
1995; Uher et al. 1998; Hanchar & Watson 2003).

ii. The accessory is distributed homogenously through-

out the igneous body and the rock composition corre-
sponds to that of the magma. For instance, the accumulation
of zircon crystals, alone or enclosed in other major rock-
forming minerals, may dramatically change the Zr concen-
tration of both melt and cumulates. In fact the plutonic
rocks very often represent solidified mixtures of a melt and
accumulated crystals. Moreover the magma composition
may be modified by both assimilation and hybridization
processes. However, as discussed by Miller et al. (2003) and
demonstrated below, the zircon saturation thermometry is,
to a large extent, robust in this respect.

iii. The accessory mineral crystallized close to the liq-

uidus, i.e. the given magma batch was saturated early in its
history. With care this can be tested by structural evi-
dence, such as the presence of inherited cores in zircons
(Fig. 1b—c) or zircon crystals enclosed in early-formed
rock-forming minerals. Moreover, the binary plots of a
fractionation index versus trace element should be charac-
terized by a negative correlation or, better still, a convex
upward trend, indicating that the given accessory started
to crystallize (Evans & Hanson 1993; Hoskin et al. 2000)
(Fig. 1a). In the latter case, the saturation temperatures for
the more basic members of the rock suite (left of the inflec-
tion point) would be far below the true liquidus (Piccoli &
Candela 1994). On the other hand, the more felsic litholo-
gies beyond the saturation point may yield good estimates
of the magma temperatures.

Testing the likelihood of the restitic material survival

In cases when a reliable independent estimate of the mag-

ma temperature exists, it can be compared with the tem-
perature obtained by the saturation thermometry. If the
‘real’ temperature is much lower, some extra component of
the accessory mineral must be present, either inherited
from the source, or accidental, from assimilated country
rocks. In this way it is possible to assess the chance to en-
counter inherited zircon cores, an information vital for the
U—Pb geochronology (e.g. Pidgeon & Aftalion 1978).
Care should be exercised in cases when the temperature
and/or bulk chemistry of the magma changed suddenly,
for instance as a consequence of magma mixing (Elburg
1996; Zeck & Williams 2002; Griffin et al. 2002).

Partial melting

For partial melting, the saturation models enable us to

constrain the concentration of the given ESC in the source
rocks, prior to and after the melt extraction. Watson & Har-
rison (1984) have presented scenarios for two conditions
of partial melting, related to the concentration of the ESC
in the source compared to the saturation level (Fig. 1b—c).

Concentration in the source > saturation level

The melt would be saturated throughout the melting

event and its ESC content buffered at a constant level, in-
dependent of the degree of melting. The amount of ESC

134

JANOUŠEK

retained by the residue would rise with increasing degrees
of the partial melting (Fig. 1b).

Concentration in the source < saturation level

The melt would be saturated in ESC only for a limited

time span, until the source is exhausted. As a conse-
quence, the higher still degrees of melting would result in
melt batches progressively depleted in the ESC. Low con-
tents of the given ESC observed in a little fractionated
magma would have two important corollaries: (1) the pa-
rental rocks were probably ESC-poor, and therefore the
source was exhausted during melting, resulting in an es-
sentially ESC-free residuum, and (2) no inheritance is like-
ly to be encountered in these intrusions (Fig. 1c).

Kinetic factors

In reality, the dissolution of the accessory minerals is

governed by a number of factors, including the tempera-
ture, crystal size, the water activity in the magma, diffusiv-
ity, absolute solubility of the ESC coupled with the
distribution and the degree of undersaturation of the melt.
For instance, experiments and numerical modelling show
that only very small zircon grains are likely to dissolve
during low-temperature anatexis (< 700 °C), while merely
relics of large crystals (> 120

µm) are likely to survive a

rather hot (850 °C) crustal fusion (e.g. Watson 1996).

The above discussion applies only to equilibrium condi-

tions, given that: 1) the dissolution rate of the accessory
was fast relative to the duration of the melting event, 2) the
accessory grains in the source were not physically isolated
from the melt as inclusions within residual minerals.

In general, the principal factor ruling the kinetics of

the zircon, apatite and monazite dissolution seems to be
the water activity in the magma. In normal, rather under-
saturated and hydrous granitic melts spanning from de-
hydratation crustal melting, the dissolution of zircon is
assumed to be fast enough for the equilibrium to be at-
tained. However, considerably slower dissolution in dry
melts (< 2 wt. % H

O) may pose a major problem (Harrison

& Watson 1983; Watson 1996). Monazite and apatite dis-
solutions are mainly controlled by sluggish LREE and
phosphorus diffusion (Harrison & Watson 1984; Rapp &
Watson 1986). In dry melts it is much too slow and thus
most of the apatite and monazite would fail to equilibrate
with the melt. While at elevated water activity the apatite
dissolves reasonably quickly, the dissolution of monazite
would still take up to several million years – too long
compared to the presumed duration of most of the crustal
anatectic events (see Spear & Pyle 2002 and references
therein). The monazite inheritance is indeed a more com-
mon phenomenon than thought previously (Copeland et
al. 1988; Harrison et al. 1995; Cocherie et al. 1998).

In the high-grade metamorphic rocks the great majority

of the larger accessory mineral grains seem to be located at
newly-formed (migrated) grain boundaries of the main
rock-forming minerals and thus probably in contact with
the partial melt during the anatexis (Watson et al. 1989).

In lower-grade rocks the accessories are mostly included
in mafic silicates (biotite and/or hornblende) that in
course of the dehydratation melting would release them
(Clemens 2003). Another factor favouring the dissolution
of the accessory phase but difficult to quantify would be
the deformation, facilitating a physical fragmentation of
the source/residue (Watson & Harrison 1984).

Overview of the saturation formulae

This section gives a synopsis of the saturation formulae.

Note that all the temperatures are absolute (K) in accord with
most of the original authors (except for Bea et al. 1992).

Zircon

The Zr saturation at a given temperature is defined by

the equation of Watson & Harrison (1983):

where:

Zr/melt

= distribution coefficient for Zr between zircon

and melt; T = absolute temperature (K); M = cationic ratio

in the melt.

Taking into account the theoretical concentration of Zr

in stoichiometric zircon (497,644 ppm) and the definition
of the distribution coefficient, we can obtain saturation
concentrations of Zr (Zr

SAT,

ppm) by:

Computing temperature using the observed Zr concen-

tration:

Fig. 2a—b shows dependence of the Zr saturation con-

centrations upon the M parameter and temperature (°C)
Eq. (1—3). Especially the Fig. 2b confirms the notion of
Miller et al. (2003) that the zircon thermometer is highly
resistant to analytical errors in Zr and major-element deter-
minations, as well as the presence of (limited) zircon in-
heritance. This represents a major advantage of the zircon
saturation thermometry for practical applications.

While the inheritance or cumulus zircon would in-

crease the estimated temperature, the crystal accumula-
tion of main rock-forming minerals would have an
opposite  effect, rendering the temperature estimates too
low. Generally speaking, the M  value of a mixture be-
tween  (1—f)*100 wt. %  liquid (denoted by subscript L) and
f*100 wt. % crystals (C) is given by (cf. Eq. 2):

(1)

. (3)

. (4)

. (5)

(2)

( )

(

)

—

0.85

3.8

497644

12900

⎟⎟⎠

⎞

⎜⎜⎝

⎛

Zr/melt

12900

— 1) +

— 0.85 (

) = — 3.80

ln(

Al.Si

Na + K + 2

Zr/melt

SAT

497644

(

) (

)

fSi

fAl

)

(1 —

)

(1 —

)

+ 2

(

) +

+ 2

)(

(1 —

135

SATURNIN R LANGUAGE FOR APPLICATION OF SATURATION MODELS IN IGNEOUS GEOCHEMISTRY

These effects were examined numerically on a small

dataset from the mainly granodioritic Kozárovice intrusion,
Central Bohemian Pluton, Czech Republic, and associated
monzonitic rocks (Lučkovice and Zalužany: Janoušek et al.
2000a,b and unpublished data). The granodiorites seem to
have been saturated in zircon throughout their history
(Fig. 2c), forming two groups of zircon saturation tempera-
tures, ~ 750 °C and ~ 780 °C (Fig. 2d).

In modelling, initial melt composition corresponding to

the sample Koz-9 (Janoušek et al. 2000b) and accumulation
of various ideal minerals were assumed (Le Maitre 1982)
(Fig. 2d). As shown by simple calculations, the presence of
feldspars-rich cumulates would pull the projection points
nearly vertically down (typical felsic granitoids: M ~ 1.3, to-
nalites: M ~ 1.9: Miller et al. 2003; Koz-9: M = 2.25; K-feld-
spar: M = 1.67; plagioclase: M = 1.67—2.50, increasing with
% An). The Ca-rich ferromagnesian minerals with negligi-
ble contents of Al

(Ca-pyroxenes, Ca-amphiboles) have

theoretically very high M (close to infinity as the denomi-
nator in Eq. (2) approaches zero) and thus would shift the M
of the mixtures strongly to the right (see Di in Fig. 2d). The
biotite would be characterized by moderate M (2.57—2.67,
only slightly increasing with the Mg contents) and its effect
will be similar to that of anorthite.

The Eq. (5) can be used to assess the effects of accumu-

lation of ferromagnesian minerals free of alkalis and calci-
um (Na

+ K

+ 2Ca

= 0). If also either Si

or Al

is zero,

the equation changes to:

(6)

with limits

Fig. 2. a—b – Zirconium saturation level (ppm) as a function of cationic ratio M = (Na + K + 2Ca) / (Al.Si) and temperature (

°C) (Watson &

Harrison 1983). c – Binary plot SiO

(wt. %)—Zr (ppm) for mainly granodioritic Kozárovice intrusion, Central Bohemian Pluton, and as-

sociated monzonitic rocks (Janoušek et al. 2000a,b and unpublished data); shown are also tentative fractionation paths. d – Binary plot M
vs. Zr (ppm) for the same dataset. Modelled are effects of up to 90 % accumulation of pure mineral phases K-feldspar (Or), albite (Ab),
Fe- and Mg-biotite (Bt, Ph), olivine with 50 % of forsterite component (Ol Fo

5 0

) and diopside (Di). Assumed initial composition corre-

sponds to the sample Koz-9; theoretical mineral compositions are from Le Maitre (1982).

)

lim(

→

)

lim(

→

and

(

)

fAl

)

(1 —

(

)

fSi

)

(1 —

136

JANOUŠEK

respectively, leading to:

if Si

= 0 and

if Al

= 0.

Thus for Ca-poor pyroxenes (En, Fs: Si = 0.5, Al = 0) and

Koz-9 as a melt composition the theoretical intersection
with the x axis (Eq. 8) would be 2.44, i.e. the points would
be shifted nearly vertically. For minerals with lower Si
(e.g. olivine, Si = 0.33, Al = 0) the intersection would be
further right (3.66).

Even though the results of such a numerical modelling

should be viewed with caution, as there is unfortunately a
dearth of experimental data for compositions other than
rhyolite glasses (Hanchar & Watson 2003), moderate accu-
mulation of albite and K-feldspar would have little effect on
the calculated saturation temperatures for normal granitic
compositions. The projection points simply move more or
less parallel to the isotherms. On the other hand, accumula-
tion of Ca-rich ferromagnesian minerals (diopsidic pyrox-
ene or magnesiohornblende) may have pronounced effects
even if as little as ~ 10 % of the cumulate is present.

Monazite

Monazite saturation temperature as a function of the

LREE concentration and water activity is defined by
(Montel 1993):

where:

D = cationic ratio

for La, Ce, Pr, Nd, Sm and Gd;

Xmz = mole fraction of the REE-phosphates in monazite;
H

O = assumed water content in the magma.

Apatite

Saturation P

concentration

For metaluminous rocks (A/CNK

≤1):

The saturation level of P

according to Harrison &

Watson (1984) (P

HW) is calculated using a combina-

tion of the expressions:

where:
T = absolute temperature (K); D

= distribution coefficient

for phosphorus between apatite and melt; SiO

= weight

fraction of silica in the melt (wt. % SiO

/100).

For peraluminous rocks (A/CNK >1):

The P

concentrations obtained according to the

model of Harrison & Watson (1984) (P

HW, see above)

can be corrected for higher phosphorus solubility in the
peraluminous melts by the following equations:

(Bea et al. 1992)

(Pichavant et al. 1992).

Apatite saturation temperatures

For metaluminous rocks (A/CNK

≤1):

(Harrison & Watson 1984).

For peraluminous rocks (A/CNK >1):

From equations (11—13) we obtain an expression that

needs to be solved for T (in K) by iterations:

(Bea et al. 1992)

as is the formula spanning from Eqs. (11), (12) and (14):

(Pichavant et al. 1992).

In Fig. 3a—b is shown dependence of the P

saturation

levels (Harrison & Watson 1984) on temperature and SiO

At the first glimpse it is obvious that the apatite saturation
thermometry for metaluminous, more siliceous composi-
tions is extremely sensible to apatite accumulation or
even small errors in the P

determinations.

However, the melts with A/CNK > 1 were shown to be

capable of dissolving much more apatite than predicted
by the formula of Harrison & Watson (1984). The reason
therefore is origin of aluminium phosphate complexes
in peraluminous aluminosilicate melts, in which P

(8)

(7)

(9)

; (10)

(11)

(13)

(14)

(15)

(16)

(17)

(12)

and

+ 9.31

— 0.5)

2.4(SiO

— 3.1 — 1

— 0.5)

00(SiO

8400 + 264

— 3.22 Si

—5900

— 1)

+ (

A/CNK

)

(

lim

→

)

(

lim

→

( )

)

— ln(

+ 0.3879

9.5 + 2.34

13318

REE

+ 2

100

Xmz

at.weight

REE

∑

— 0.5)

2.4(SiO

— 3.1 — 1

— 0.5)

00(SiO

8400 + 264

)

(

273

—

— 1)

6429(

ACNK

HWe

SiO

—

5900

—

— 1)

(

A/CNK

)

—

2.4(SiO

+ 3.1 + 1

— 0.5)

00(SiO

8400 + 264

⎟

⎟
⎠

⎞

⎜

⎜
⎝

⎛

T(K)

= 42

—0.5)

SiO

+3.1+12.4(

—0.5)

(SiO

8400+26400

—

— 273.15

6429(ACNK—

138

JANOUŠEK

bonded with Al

3 +

to form AlPO

species (Mysen et al.

1999 and references therein). These effects are shown in
Fig. 3c—d for the two most elaborated models, published
by Bea et al. (1992) and Pichavant et al. (1992). The latter
seems to give solubilities significantly lower than the
former. Wolf & London (1994) found the phosphorus sol-
ubility to be linearly dependent on the A/CNK values:

However, their experiments were conducted only at a

single temperature (T = 750 °C).

The apatite thermometry in felsic, strongly peraluminous

melts is furthermore flawed by the fact that feldspars can in-
corporate large amounts of phosphorus via a berlinite sub-
stitution (e.g. London et al. 1990). In Ca-poor magmas
plagioclase would crystallize early, further increasing the
P/Ca ratio in the melt. In extreme cases of low-T, highly
fractionated, Ca-poor and strongly peraluminous magmas,
the formation of aluminium phosphate complexes may pro-
hibit the apatite saturation completely, with all phosphorus
being allocated to other phosphates and/or feldspars (Picco-
li & Candela 2002). The calculated apatite saturation tem-
peratures in these cases would be of course meaningless.

Saturnin: a R language script for saturation

calculations

System overview

To the author’s knowledge, a flexible and comprehensi-

ble tool that would perform all the necessary saturation cal-
culations has so far been lacking. As shown by Janoušek
(2000) and Grunsky (2002), the freeware R language (R De-
velopment Core Team, 2003) [1] is an ultimate environ-
ment for writing geochemical recalculation and plotting
routines that are compact and platform independent. Within

this rich environment, a set of scripts, called Saturnin, has
been developed. They have been included, in a slightly
modified form, as a plugin for the larger system called GCD-
kit for Windows (Janoušek et al. 2003) [2]. The scripts have
been tested in R for Windows, version 2.1.1. They should
run not only under Windows 9x/2000/NT/XP (higher ver-
sions are strongly recommended) but also, theoretically
nearly without modifications, on Mac OS X or Linux.

Obtaining, installing and running Saturnin

If not present already, the R environment needs to be in-

stalled first from the nearest mirror [1]. For Windows, a
complete distribution is packed in a single file
‘R-X.X.X-win32.exe’, where X.X.X is the version number
(e.g. R-2.2.1-win32.exe). Run the executable file and select
the required items as well as the target directory.

The second step would be to download a single text file

with the code of Saturnin from [3], the Geologica Car-
pathica Web page [4] or request it by email from the au-
thor. In Windows, the code is started by invoking the
File|Source R code command from the RGUI menu, in
other systems using the standard R command source. It
might be wise to change the working directory, in Win-
dows invoking the menu item File|Change dir…, other-
wise using the R command setwd:

setwd(”c:/data”)

The programme requires plain text (ASCII) files in

which each row represents one sample and individual
items (major and minor elements, Zr, and, or LREE con-
centrations) are in columns delimited by tabulators
(Fig. 4). The first line contains labels for the data columns
(except that for the first column that is automatically as-
sumed to contain the sample names). The first row there-
fore should have one item less than all the following ones.
The columns can be given in an arbitrary order, missing
values (replaced in R by ‘NA’, standing for ‘not available’)
for elements/oxides not directly involved in the calcula-

Fig. 4. Example of a tab-delimited text file formatted for Saturnin.

(18)

A/CNK

—

139

SATURNIN R LANGUAGE FOR APPLICATION OF SATURATION MODELS IN IGNEOUS GEOCHEMISTRY

tions are allowed. However, column names must have sim-
ple form to be treated properly by automatic routines (i.e.
‘SiO2’ instead of ‘SiO2 [wt. %]’ etc.) and must be unique,
as must be the sample names. When a data set is loaded
into the system, all numerical data should be allocated in
a data frame ‘WR’ as follows:

load.file()
Generally speaking, the functions in R can be called

simply by typing their names with required parameters in
parentheses. If some of the parameters are omitted, defaults
from the function’s definitions are taken:

zr.saturation(T=850)

zr.saturation()

x<–c(3,3,3,2.5,2,2,7,3)

# assuming 8 samples

mz.saturation(H2O=x)
The results of calculations are printed on the screen and

stored in a matrix/vector ‘results’.

In Windows, this variable can be subsequently copied

to a clipboard using the function r2clip:

r2clip()
and pasted into any other programme, such as MS Ex-

cel, for further processing. Note that all temperatures in the
output have been converted to the centigrade (°C).

Conclusions

1 – Zircon saturation thermometer is rather robust to

analytical errors, presence of minor inheritance, as well as
moderate accumulation of feldspars and some ferromagne-
sian minerals (Ca-poor pyroxenes, to a lesser extent oliv-
ine). On the other hand, addition of Ca-rich pyroxene or
Ca-rich amphibole cumulate renders quickly the calculat-
ed temperatures too low.

2 – Apatite saturation thermometry is a poor tool for

felsic metaluminous rocks, being exceedingly sensitive to
small analytical errors in the phosphorus determination as
well as the presence of extraneous apatite. An even stron-
ger element of uncertainty is added for peraluminous
lithologies, when the increased apatite solubility is poorly
constrained by the current models and where other P-bear-
ing minerals, most importantly feldspars and monazite,
come to the play.

3 – Saturnin is a newly developed software that performs

otherwise rather tedious zircon, monazite and apatite satura-
tion calculations in igneous geochemistry. It is written in a
freeware R language and designed to be platform indepen-
dent. While the Windows users would probably find it easier
to access the programme as a part of the larger package GCD-
kit (http://www.gla.ac.uk/gcdkit), on Macintosh and Linux it
can be invoked as a stand-alone application.

4 – The programme can be downloaded from the

WWW or can be obtained upon request from the author,
who hopes that it will promote the power and beauty of
the R language in a wider geochemical community.

Acknowledgments:

This manuscript originated during the

research stay of VJ at the Institute of Mineralogy, Univer-

sity of Salzburg, in the framework of the FWF Project
15133—GEO (to Fritz Finger). Generous support from the
Czech Grant Agency (GAČR 205/03/0040) is gratefully
acknowledged. The author is indebted to the R Core Team
and other people from the R community for their efforts in
continuous developing R. The paper benefited from help-
ful reviews by Igor Petrík and two anonymous colleagues.
The author also thanks Igor Broska (Academy of Sciences,
Bratislava) for organizing a workshop explaining the prin-
ciples of saturation calculations, Fritz Finger for com-
ments on an early version of the manuscript, Colin Farrow
(University of Glasgow) and Vojtěch Erban (Czech Geo-
logical Survey, Prague) for interesting discussions about
the R language. Vojtěch Erban also helped considerably
by testing the code. Moreover, Angela Storkey (La Trobe
University, Bundoora, Australia) assisted with painstaking
debugging of the complex apatite saturation script.

References

Ayres M. & Harris N. 1997: REE fractionation and Nd-isotope dis-

equilibrium during crustal anatexis: constraints from Himalay-
an leucogranites. Chem. Geol. 139, 249—269.

Bacon C.R. 1989: Crystallization of accessory phases in magmas by

local saturation adjacent to phenocrysts. Geochim. Cosmochim.
Acta 53, 1055—1066.

Baker D.R., Conte A.M., Freda C. & Ottolini L. 2002: The effect of

halogens on Zr diffusion and zircon dissolution in hydrous
metaluminous granitic melts. Contr. Mineral. Petrology 142,
666—678.

Bea F. 1996: Residence of REE, Y, Th and U in granites and crust-

al protoliths; implications for the chemistry of crustal melts.
J. Petrology 37, 521—552.

Bea F., Fershtater G.B. & Corretgé L.G. 1992: The geochemistry of

phosphorus in granite rocks and the effects of aluminium.
Lithos 29, 43—56.

Broska I., Gerdes A., Haunschmid B., Schindlmayr A. & Finger

F. 1995: Magma temperatures in the southern Bohemian
batholith estimated on the basis of zircon solubility. Terra
Abstracts EUG 8, Strasbourg, 143.

Chappell B.W., Bryant C.J., Wyborn D., White A.J.R. & Will-

iams I.S. 1998: High- and low-temperature I-type granites.
Resource Geol. 48, 225—235.

Cherniak D.J. 2000: Rare earth element diffusion in apatite.

Geochim. Cosmochim. Acta 64, 3871—3885.

Cherniak D.J. & Watson E.B. 2003: Diffusion in zircon. In: Han-

char J.M. & Hoskin P.W.O. (Eds.): Zircon. Mineralogical So-
ciety of America and Geochemical Society Reviews in
Mineralogy and Geochemistry 53, Washington, 113—143.

Cherniak D.J., Watson E.B., Grove M. & Harrison T.M. 2004: Pb

diffusion in monazite: a combined RBS/SIMS study. Geochim.
Cosmochim. Acta 68, 829—840.

Clemens J.D. 2003: S-type granitic magmas – petrogenetic issues,

models and evidence. Earth Sci. Rev. 61, 1—18.

Cocherie A., Legendre O., Peucat J.J. & Kouamelan A.N. 1998:

Geochronology of polygenic monazites constrained by in situ
microprobe Th—U—total lead determination: implications for
lead behaviour in monazite. Geochim. Cosmochim. Acta 62,
2475—2497.

Copeland P., Parrish R.R. & Harrison T.M. 1988: Identification of

inherited radiogenic Pb in monazite and its implications for
U—Pb systematics. Nature 333, 760—763.

Elburg M.A. 1996: U—Pb ages and morphologies of zircon and mi-

crogranitoid enclaves and peraluminous host granite: evidence
for magma mingling. Contr. Mineral. Petrology 123, 177—189.

140

JANOUŠEK

Evans O.C. & Hanson G.N. 1993: Accessory-mineral fractionation

of rare-earth element (REE) abundances in granitoid rocks.
Chem. Geol. 110, 69—93.

Griffin W.L., Wang X., Jackson S.E., Pearson N.J., O’Reilly S.Y.,

Xu X. & Zhou X. 2002: Zircon chemistry and magma mixing,
SE China: In-situ analysis of Hf isotopes, Tonglu and Pingtan
igneous complexes. Lithos 61, 237—269.

Gromet L.P. & Silver L.T. 1983: Rare earth element distribution

among minerals in a granodiorite and their petrogenetic impli-
cations. Geochim. Cosmochim. Acta 47, 925—939.

Grunsky E.C. 2002: R: a data analysis and statistical programming

environment – an emerging tool for the geosciences. Com-
put. and Geosci. 28, 1219—1222.

Hanchar J.M. & Miller C.F. 1993: Zircon zonation patterns as re-

vealed by cathodoluminescence and backscattered electron
images: implications for interpretation of complex crustal his-
tories. Chem. Geol. 110, 1—13.

Hanchar J.M. & Rudnick R.L. 1995: Revealing hidden structures:

the application of cathodoluminescence and back-scattered
electron imaging to dating zircons from lower crustal xeno-
liths. Lithos 36, 289—303.

Hanchar J.M. & Watson E.B. 2003: Zircon saturation thermometry.

In: Hanchar J.M. & Hoskin P.W.O. (Eds.): Zircon. Mineralog-
ical Society of America and Geochemical Society Reviews in
Mineralogy and Geochemistry 53, Washington, 89—112.

Hanson G.N. & Langmuir C.H. 1978: Modelling of major elements

in mantle-melt systems using trace element approaches.
Geochim. Cosmochim. Acta 42, 725—741.

Harrison T.M. & Watson E.B. 1983: Kinetics of zircon dissolution

and zirconium diffusion in granitic melts of variable water
content. Contr. Mineral. Petrology 84, 66—72.

Harrison T.M. & Watson E.B. 1984: The behavior of apatite during

crustal anatexis: equilibrium and kinetic considerations.
Geochim. Cosmochim. Acta 48, 1467—1477.

Harrison T.M., McKeegan K.D. & Le Fort P. 1995: Detection of

inherited monazite in the Manaslu leucogranite by

208

Pb/

232

ion microprobe dating; crystallization age and tectonic impli-
cations. Earth Planet. Sci. Lett. 133, 271—282.

Hoskin P.W.O., Kinny P.D., Wyborn D. & Chappell B.W. 2000:

Identifying accessory mineral saturation during differentiation
in granitoid magmas: an integral approach. J. Petrology 41,
1365—1396.

Janoušek V. 2000: R – an alternative to spreadsheets and special

software for geochemical calculations and plotting. Geolines
10, 34—35.

Janoušek V., Bowes D.R., Braithwaite C.J.R. & Rogers G. 2000a: Mi-

crostructural and mineralogical evidence for limited involve-
ment of magma mixing in the petrogenesis of a Hercynian
high-K calc-alkaline intrusion: the Kozárovice granodiorite,
Central Bohemian Pluton, Czech Republic. Trans. Roy. Soc. Ed-
inb., Earth Sci. 91, 15—26.

Janoušek V., Bowes D.R., Rogers G., Farrow C.M. & Jelínek E.

2000b: Modelling diverse processes in the petrogenesis of a
composite batholith: the Central Bohemian Pluton, Central Eu-
ropean Hercynides. J. Petrology 41, 511—543.

Janoušek V., Farrow C.M. & Erban V. 2003: GCDkit: new PC soft-

ware for interpretation of whole-rock geochemical data from
igneous rocks. Geochim. Cosmochim. Acta A67, 186.

Le Maitre R.W. 1982: Numerical petrology. Elsevier, Amsterdam,

1—281.

London D. 1992: Phosphorus in S-type magmas: the P

content

of feldspars from peraluminous granites, pegmatites and rhyo-
lites. Amer. Mineralogist 77, 126—145.

London D., Černý P., Loomis J.L. & Pan J.J. 1990: Phosphorus in

alkali feldspars of rare-element granitic pegmatites. Canad.
Mineralogist 28, 771—786.

Miller C.F. & Mittlefehldt D.W. 1984: Extreme fractionation in fel-

sic magma chambers; a product of liquid-state diffusion or
fractional crystallization? Earth Planet. Sci. Lett. 68, 151—158.

Miller C.F., McDowell S.M. & Mapes R.W. 2003: Hot and cold

granites? Implications of zircon saturation temperatures and
preservation of inheritance. Geology 31, 529—532.

Montel J.M. 1993: A model for monazite/melt equilibrium and ap-

plication to the generation of granitic magmas. Chem. Geol.
110, 127—146.

Mysen B.O., Holtz F., Pichavant M., Beny J.M., Montel J.M. & Holtz

F. 1999: The effect of temperature and bulk composition on the
solution mechanism of phosphorus in peraluminous haplogra-
nitic magma. Amer. Mineralogist 84, 1336—1345.

Nasdala L., Zhang M., Kempe U., Panczer G., Gaft M., Andrut M.

& Plötze M. 2003: Spectroscopic methods applied to zircon.
In: Hanchar J.M. & Hoskin P.W.O. (Eds.): Zircon. Mineralog-
ical Society of America and Geochemical Society Reviews in
Mineralogy and Geochemistry 53, Washington, 427—467.

Paterson B.A., Stephens W.E., Rogers G., Williams I.S., Hinton

R.W. & Herd D.A. 1992: The nature of zircon inheritance in
two granite plutons. Trans. Roy. Soc. Edinb., Earth Sci. 83,
459—471.

Piccoli P. & Candela P. 1994: Apatite in felsic rocks: a model for

the estimation of initial halogen concentrations in the Bishop
Tuff (Long Valley) and Tuolumne Intrusive Suite (Sierra Ne-
vada batholith) magmas. Amer. J. Sci. 294, 92—135.

Piccoli P. & Candela P. 2002: Apatite in igneous systems. In: Kohn

M.J., Rakovan J. & Hughes J.M. (Eds.): Phosphates:
Geochemical, geobiological, and materials importance. Miner-
alogical Society of America Reviews in Mineralogy and
Geochemistry 48, Washington, 255—292.

Pichavant M., Montel J.M. & Richard L.R. 1992: Apatite solubility

in peraluminous liquids: experimental data and extension of
the Harrison—Watson model. Geochim. Cosmochim. Acta 56,
3855—3861.

Pidgeon R.T. & Aftalion M. 1978: Cogenetic and inherited zircon

U—Pb systems in granites: Palaeozoic granites of Scotland and
England. In: Bowes D.R. & Leake B.E. (Eds.): Crustal evolution
in northwestern Britain and adjacent regions. Geol. J. 10, 183—220.

Poller U. 2000: A combination of single zircon dating by TIMS

and cathodoluminescence investigations on the same grain: the
CLC method – U-Pb geochronology for metamorphic rocks.
In: Pagel M., Barbin V., Blanc P. & Ohnenstetter D. (Eds.):
Cathodoluminescence in geosciences. Springer, Berlin, 1—514.

Rapp R.P. & Watson E.B. 1986: Monazite solubility and dissolution

kinetics; implications for the thorium and light rare earth
chemistry of felsic magmas. Contr. Mineral. Petrology 94,
304—316.

Sawka W.N. 1988: REE and trace element variations in accessory

minerals and hornblende from the strongly zoned McMurry
Meadows Pluton, California. Trans. Roy. Soc. Edinb., Earth
Sci. 79, 157—168.

Spear F.S. & Pyle J.M. 2002: Apatite, monazite and xenotime in

metamorphic rocks. In: Kohn M.J., Rakovan J. & Hughes J.M.
(Eds.): Phosphates: Geochemical, geobiological, and materials
importance. Mineralogical Society of America Reviews in Min-
eralogy and Geochemistry 48, Washington 293—336.

Uher P., Breiter K., Klečka M. & Pivec E. 1998: Zircon in highly

evolved Hercynian Homolka granite, Moldanubian zone,
Czech Republic: Indicator of magma source and petrogenesis.
Geol. Carpathica 49, 151—160.

Watson E.B. 1996: Dissolution, growth and survival of zircons dur-

ing crustal fusion; kinetic principles, geological models and
implications for isotopic inheritance. Trans. Roy. Soc. Edinb.,
Earth Sci. 87, 43—56.

Watson E.B. & Capobianco C.J. 1981: Phosphorus and the rare

earth elements in felsic magmas: an assessment of the role of
apatite. Geochim. Cosmochim. Acta 45, 2349—2358.

Watson E.B. & Harrison T.M. 1983: Zircon saturation revisited:

temperature and composition effects in a variety of crustal
magma types. Earth Planet. Sci. Lett. 64, 295—304.

Watson E.B. & Harrison T.M. 1984: Accessory minerals and the

geochemical evolution of crustal magmatic systems: a summa-
ry and prospectus of experimental approaches. Phys. Earth

141

SATURNIN R LANGUAGE FOR APPLICATION OF SATURATION MODELS IN IGNEOUS GEOCHEMISTRY

Planet. Inter. 35, 19—30.

Watson E.B., Vicenzi E.P. & Rapp R.P. 1989: Inclusion/host rela-

tions involving accessory minerals in high-grade metamorphic
and anatectic rocks. Contr. Mineral. Petrology 101, 220—231.

Watt G.R. & Harley S.L. 1993: Accessory phase controls on the

geochemistry of crustal melts and restites produced during wa-
ter-undersaturated partial melting. Contr. Mineral. Petrology
114, 550—566.

Wolf M.B. & London D. 1994: Apatite dissolution into peralumi-

nous haplogranitic melts; an experimental study of solubilities
and mechanisms. Geochim. Cosmochim. Acta 58, 4127—4146.

Wolf M.B. & London D. 1995: Incongruent dissolution of REE- and

Sr-rich apatite in peraluminous granitic liquids; differential ap-
atite, monazite, and xenotime solubilities during anatexis.
Amer. Mineralogist 80, 765—775.

Zeck H.P. & Williams I.S. 2002: Inherited and magmatic zircon

from Neogene Hoyazo cordierite dacite, SE Spain – anatectic
source rock proven

ance and magmatic evolution. J. Petrology

43, 1089—1104.

Internet references

[1] The Comprehensive R Archive Network (CRAN).
          URL  http://cran.r-project.org.
[2]  Geochemical Data Toolkit (GCDkit) for Windows.
          URL  http://www.gla.ac.uk/gcdkit.
[3] Saturnin.
          URL  http://www.gla.ac.uk/gcdkit/saturnin.
[4]   Geologica Carpathica.
   URL http://www.geologicacarpathica.sk.

Auxiliary functions

require(”MASS”)

major<-c(”SiO2”,”TiO2”,”Al2O3”,”Fe2O3”,”FeO”,”MnO”,”MgO”,

”CaO”,”Na2O”,”K2O”,”P2O5”,”Li2O”)

LREE<-c(”La”,”Ce”,”Pr”,”Nd”,”Sm”,”Gd”)

WR<-data.matrix(0)

# Loading data file
load.file

<-function(){

if(.Platform$OS.type==”windows”){

filename<-file.choose()

}else{

filename<-readline(”Enter the filename: ”)

}

x<-read.table(filename,sep=”\t”,dec=”.”,

check.names=F, fill=T) # For decimal commas: dec=”,”

# Ensures that all the necessary variables are there,

# even if empty

add.on

<-function(param,what){

if(any(colnames(WR)==param))return(WR)

options(show.error.messages = F)

try(WR<-cbind(WR,what))

try(colnames(WR)[ncol(WR)]<-param)

options(show.error.messages = T)

return(WR)

}

col.names<-c(major,LREE,”Zr”)

matrix(nrow=nrow(x),ncol=length(col.names),

dimnames=list(rownames(x),col.names))
WR<-cbind(x,y[,setdiff(colnames(y),colnames(x))])

WR<-add.on(”Li2O”,WR[,”Li”]/0.46452/1e4)

# Recalculates Li (ppm) to LiO2 (wt%) if necessary

WR<-add.on(”FeO”,WR[,”FeOt”])

WR<-add.on(”FeO”,WR[,”FeO*”])

WR<-add.on(”A/CNK”,acnk(millicat(WR)))

WR[WR<=0]<-NA # Missing values

print(WR)

WR<<-WR

return(WR)

}

# Calculates millications
millicat

<-function(what=WR){

MW<c(60.0848,79.8988,101.96128,159.6922,71.88464,70.9374,

40.3044,56.0794,61.97894,94.1954,141.94452,29.8814)

# Molecular weights for ’major’

names(MW)<-major

fact<-c(1,1,2,2,1,1,1,1,2,2,2,2)

# Number of cations per molecule

names(fact)<-major

z<-t(apply(what[,major],1,function(x){x/MW[major]*

fact[major]*1000}))

return(z)

}

# Calculates A/CNK
acnk

<-function(what){

z<-what[,”Al2O3”]/(what[,”Na2O”]+what[,”K2O”]+

2*what[,”CaO”])

return(z)

}

# Normalizes a matrix to a given sum
normalize2total

<-function(what,total=100){

z<-t(apply(what,1,function(x,y){x/sum(x,na.rm=T)*y},

y=total))

return(z)

}

# Filters out from matrix ‘where’ rows in which exist all columns
specified in ‘what’

filter.out

<-function(where,what){

i<-apply(where[,what]),1,function(x){all(!is.na(x))})

z.names<-rownames(where)[i]

z<-as.matrix(where[i,what])

rownames(z)<-z.names

colnames(z)<-what; mode(z)<-”numeric”

return(z)

}

r2clip

<-function(what=results){

if(is.null(what)) stop(”No data available!”)

filename<-file(”clipboard”,open=”w”)

what<-as.matrix(what)

if(ncol(what)==1)colnames(what)<-””

write.matrix(cbind(rownames(what),what),filename,sep=”\t”)

close(filename)

}

Appendix: saturnin.r — programme code

[See http://www.gla.ac.uk/gcdkit/saturnin/saturnin.pdf for detailed documentation]

142

JANOUŠEK

Zircon saturation

Parameters:
cats   numeric matrix; whole-rock data recast to millications
T   assumed magma temperature (default = 750 °C)
Zr   vector with Zr concentrations (ppm)
Returns a matrix ‘results’ with the following columns:
M     cationic ratio (Eq. 2)
Zr   observed Zr concentrations (ppm) (Eqs. 1—3)
TZr.sat.C zircon saturation temperatures (°C) (Eq. 4)

zr.saturation

<-function(cats=millicat(WR),T=750,

Zr=subset(WR, Zr>0, ”Zr”)){

T<-T+273.15

cats<-cats[rownames(Zr),]

x<-normalize2total(cats)

M<-(x[,”Na2O”]+x[,”K2O”]+2*x[,”CaO”])/(x[,”Al2O3”]*

x[,”SiO2”])*100

DZr1<-exp(-3.8-0.85*(M-1)+12900/T)

Zr.sat<-497644/DZr1

DZr<-497644/Zr

TZr.sat.C<-12900/(log(DZr)+3.8+0.85*(M-1))-273.15

y<-cbind(M,Zr,round(Zr.sat,1),round(TZr.sat.C,1))

colnames(y)<-c(”M”,”Zr”,”Zr.sat”,”TZr.sat.C”)

results<<-y

return(y)

}

Monazite saturation

Parameters:

cats   numeric matrix; whole-rock data recast to millications
H2O         assumed water contents of the magma (default = 3 wt. %)
Xmz         mole fraction of the REE-phosphates in monazite
  (default = 0.83)

Returns a matrix ‘results’ with the following columns:

Dmz cationic ratio (Eq. 10)
Tmz.sat.C monazite saturation temperature (°C) (Eq. 9)

mz.saturation

<-function(cats=millicat(WR),H2O=3,

Xmz=0.83){

MW.REE<-c(138.9055,140.12,140.9077,144.24,151.4,

154.25)# Atomic weights for LREE

names(MW.REE)<-LREE

REE<-filter.out(WR,LREE)

cats<-cats[rownames(REE),]

# Get only samples for which all LREE are available

ree<-t(t(REE)/MW.REE)

reex<-apply(ree,1,sum,na.rm=T)/Xmz

x<-normalize2total(cats,100)

D<-(x[,”Na2O”]+x[,”K2O”]+2*x[,”CaO”])/x[,”Al2O3”]*

1/(x[,”Al2O3”]+x[,”SiO2”])*100

T.calc<-13318/(9.5+2.34*D+0.3879*sqrt(H2O)-log(reex))

-273.15

y<-cbind(D,round(T.calc,1))

colnames(y)<-c(”Dmz”,”Tmz.sat.C”)

results<<-y

return(y)

}

Apatite saturation

Parameters:

Si SiO

contents in the melt (wt. %)

ACNK vector with A/CNK (mol %) values
P2O5 vector with P

concentrations (wt. %)

T assumed magma temperature (default = 750 °C)

Returns a matrix ‘results’ with the following columns:

A/CNK A/CNK values
Tap.sat.C.H+W saturation T of Harrison & Watson (1984) (°C)
(Eq. 15)

For peraluminous rocks (A/CNK > 1) also:

Tap.sat.C.Bea saturation T of Bea et al. (1992) (°C) (Eq. 16)
Tap.sat.C.Pich saturation T of Pichavant et al. (1992) (°C) (Eq. 17)

ap.saturation

<-function(Si=WR[,”SiO2”],

ACNK=WR[,”A/CNK”],P2O5=data.matrix(WR)[,”P2O5”],T=750){

Si<-Si/100

T<-T+273.15

# Harrison and Watson (1984)

A<-8400+(Si-0.5)*26400

B<-3.1+12.4*(Si-0.5)

D.HW<-exp(A/T-B)

P2O5.HW<-42/D.HW

T.HW<-A/(log(42/P2O5)+B)-273.15

# A general routine that solves non-linear equation

# for T (deg C) by bisection method

solve.T

<-function(fun,tmin=0,tmax=NULL){

T.calc<-NULL

for(i in 1:length(Si)){

if(ACNK[i]>1){

ttold<-0;tt<-1

if(is.null(tmax))tt.max<-T.HW[i]+273 else

tt.max<-tmax

# H+W temperature is a feasible max. estimate

tt.min<-tmin

while(abs(ttold-tt)>0){

ttold<-tt

tt<-(tt.max-tt.min)/2+tt.min

expr<-gsub(”Si”,Si[i],fun)

expr<-gsub(”ACNK”,ACNK[i],expr)

expr<-gsub(”T”,tt,expr)

pp<-eval(parse(text=as.expression(expr)))

if(pp>P2O5[i])tt.max<-tt else tt.min<-tt

}

T.calc[i]<-tt

}else{

T.calc[i]<-NA

}

return(T.calc-273.15)

}

# Pichavant et al. (1992)

T.PV<-solve.T(”42/exp((8400+(Si-0.5)*26400)/T-3.1-

12.4*(Si-0.5))+(ACNK-1)*exp(-5900/T-3.22*Si+9.31)”)

# Bea et al. (1992)

T.Bea<-solve.T(”42*exp(((ACNK-1)*6429)/(T-273.15)-

(8400+(Si-0.5)*26400)/T+3.1+12.4*(Si-0.5))”)

y<-cbind(T.HW,T.Bea,T.PV);y<-round(y,1)

y<-cbind(round(ACNK,2),y)

colnames(y)<-c(”A/CNK”,”Tap.sat.C.H+W”,

”Tap.sat.C.Bea”,”Tap.sat.C.Pich”)

rownames(y)<-names(P2O5)

results<<-y

return(y)

}