berman: Thermodynamic Properties of Minerals

Description Usage Arguments Details Value References Examples

View source: R/berman.R

Description

Calculate thermodynamic properties of minerals using the equations of Berman (1988).

Usage

1
2
  berman(name, T = 298.15, P = 1, thisinfo = NULL, check.G = FALSE,
         calc.transition = TRUE, calc.disorder = TRUE, units = "cal")

Arguments

name

character, name of mineral

T

numeric, temperature(s) at which to calculate properties (K)

P

numeric, pressure(s) at which to calculate properties (bar)

thisinfo

dataframe, row for mineral from thermo()$obigt

check.G

logical, check consistency of G, H, and S?

calc.transition

logical, include calculation of polymorphic transition properties?

calc.disorder

logical, include calculation of disordering properties?

units

character, energy units, cal or J

Details

This function calculates the thermodynamic properties of minerals at high \P and \T using equations given by Berman (1988). These minerals should be listed in thermo()$obigt with the state cr and chemical formula, and optionally an abbreviation and references, but all other properties set to NA.

The standard state thermodynamic properties and parameters for the calculations are stored in data files under extdata/Berman, or can be read from a user-created file (if available) named berman.csv in the working directory.

The equation used for heat capacity is \Cp = k0 + k1*\T\S-0.5 + k2*\T\S-2 + k3*\T\S-3 + k4*\T\S-1 + k5*\T + k6*\T\S2. This is an extended form Eq. 4 of Berman (1988) as used in the winTWQ program (Berman, 2007). The equation used for volume is \IV(\P, \T) / \IV(1 bar, 298.15 K) = 1 + v1 * (\T - 298.15) + v2 * (\T - 298.15)\S2 + v3 * (\P - 1) + v4 * (\P - 1)\S2 (Berman, 1988, Eq. 5, with terms reordered to follow winTWQ format). The equations used for lambda transitions follow Eqs. 8-14 of Berman (1988). The equation used for the disorder contribution between Tmin and Tmax is \Cp[dis] = d0 + d1*\T\S-0.5 + d2*\T\S-2 + d3*\T + d4*\T\S2 (Berman, 1988, Eq. 15). The parameters correspond to Tables 2 (GfPrTr, HfPrTr, SPrTr, VPrTr), 3a (k0 to k3), 4 (v1 to v4), 3b (transition parameters: Tlambda to dTH), and 5 (disorder parameters: Tmax, Tmin, d1 to d4 and Vad) of Berman (1988). Following the winTWQ data format, multipliers are applied to the volume parameters only (see below). Note that VPrTr is tabulated in J bar\S-1 mol\S-1, which is equal to 10 cm\S3 mol\S-1.

A value for GfPrTr is not required and is only used for optional checks (see below). Numeric values (possibly 0) should be assigned for all of HfPrTr, SPrTr, VPrTr, k0 to k6 and v1 to v4. Missing (or NA) values are permitted for the transition and disorder parameters, for minerals where they are not used. The data files have the following 30 columns:

name mineral name (must match an entry with a formula but NA properties in thermo()$obigt)
GfPrTr standard Gibbs energy at 298.15 K and 1 bar (J mol\S-1) (Benson-Helgeson convention)
HfPrTr standard enthalpy at 298.15 K and 1 bar (J mol\S-1)
SPrTr standard entropy at 298.15 K and 1 bar (J mol\S-1 K\S-1)
VPrTr standard volume at 298.15 K and 1 bar (J mol\S-1)
k0 ... k6 k0 (J mol\S-1 K\S-1) to k6
v1 v1 (K\S-1) * 10\S5
v2 v2 (K\S-2) * 10\S5
v3 v3 (bar\S-1) * 10\S5
v4 v4 (bar\S-2) * 10\S8
Tlambda Tλ (K)
Tref Tref (K)
dTdP d\T / d\P (K bar\S-1)
l1 l1 ((J/mol)\S0.5 K\S-1)
l2 l2 ((J/mol)\S0.5 K\S-2)
DtH \DeltaTt\IH (J mol\S-1)
Tmax temperature at which phase is fully disordered (TD in Berman, 1988) (K)
Tmin reference temperature for onset of disordering (t in Berman, 1988) (K)
d0 ... d4 d0 (J mol\S-1 K\S-1) to d4
Vad constant that scales the disordering enthalpy to volume of disorder (d5 in Berman, 1988)

The function outputs apparent Gibbs energies according to the Benson-Helgeson convention (\Delta\IG = \Delta\IH - \T\Delta\IS) using the entropies of the elements in the chemical formula of the mineral to calculate \Delta\IS (cf. Anderson, 2005). If check.G is TRUE, the tabulated value of GfTrPr (Benson-Helgeson) is compared with that calculated from HfPrTr - 298.15*DSPrTr (DSPrTr is the difference between the entropies of the elements in the formula and SPrTr in the table). A warning is produced if the absolute value of the difference between tabulated and calculated GfTrPr is greater than 1000 J/mol.

Providing thisinfo avoids searching for the mineral in thermo()$obigt, potentially saving some running time. If the function is called with missing name, the parameters for all available minerals are returned.

Value

A data frame with T (K), P (bar), G, H, S, and Cp expressed in the given units (cal or J), and V (cm\S3 mol\S-1).

References

Anderson, G. M. (2005) Thermodynamics of Natural Systems, 2nd ed., Cambridge University Press, 648 p. http://www.worldcat.org/oclc/474880901

Berman, R. G. (1988) Internally-consistent thermodynamic data for minerals in the system Na\s2O-K\s2O-CaO-MgO-FeO-Fe\s2O\s3-Al\s2O\s3-SiO\s2-TiO\s2-H\s2O-CO\s2. J. Petrol. 29, 445-522. https://doi.org/10.1093/petrology/29.2.445

Berman, R. G. and Aranovich, L. Ya. (1996) Optimized standard state and solution properties of minerals. I. Model calibration for olivine, orthopyroxene, cordierite, garnet, and ilmenite in the system FeO-MgO-CaO-Al\s2O\s3-TiO\s2-SiO\s2. Contrib. Mineral. Petrol. 126, 1-24. https://doi.org/10.1007/s004100050233

Berman, R. G. (2007) winTWQ (version 2.3): A software package for performing internally-consistent thermobarometric calculations. Open File 5462, Geological Survey of Canada, 41 p. https://doi.org/10.4095/223425

Helgeson, H. C., Delany, J. M., Nesbitt, H. W. and Bird, D. K. (1978) Summary and critique of the thermodynamic properties of rock-forming minerals. Am. J. Sci. 278-A, 1–229. http://www.worldcat.org/oclc/13594862

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
# other than the formula, the parameters aren't stored in
# thermo()$obigt, so this shows NAs
info(info("quartz", "cr"))
# properties of alpha-quartz (aQz) at 298.15 K and 1 bar
berman("quartz")
# Gibbs energies of aQz and coesite at higher T and P
T <- seq(200, 1300, 100)
P <- seq(22870, 31900, length.out=length(T))
G_aQz <- berman("quartz", T=T, P=P)$G
G_Cs <- berman("coesite", T=T, P=P)$G
# that is close to the univariant curve (Ber88 Fig. 4),
# so the difference in G is close to 0
DGrxn <- G_Cs - G_aQz
stopifnot(all(abs(DGrxn) < 100))

### compare mineral stabilities in the Berman and Helgeson datasets
### on a T - log(K+/H+) diagram, after Sverjensky et al., 1991
### (doi:10.1016/0016-7037(91)90157-Z)
## set up the system: basis species
basis(c("K+", "Al+3", "quartz", "H2O", "O2", "H+"))
# use pH = 0 so that aK+ = aK+/aH+
basis("pH", 0)
# load the species
species(c("K-feldspar", "muscovite", "kaolinite",
          "pyrophyllite", "andalusite"), "cr")
## start with the data from Helgeson et al., 1978
add.obigt("SUPCRT92")
# calculate affinities in aK+ - temperature space
# exceed.Tr: enable calculations above stated temperature limit of pyrophyllite
res <- 400
a <- affinity(`K+` = c(0, 5, res), T = c(200, 650, res), P = 1000, exceed.Ttr = TRUE)
# make base plot with colors and no lines
diagram(a, xlab = ratlab("K+", molality = TRUE), lty = 0, fill = "terrain")
# add the lines, extending into the low-density region (exceed.rhomin = TRUE)
a <- affinity(`K+` = c(0, 5, res), T = c(200, 650, res), P = 1000, 
              exceed.Ttr = TRUE, exceed.rhomin = TRUE)
diagram(a, add = TRUE, names = FALSE, col = "red", lwd = 1.5)
# the list of references:
ref1 <- thermo.refs(species()$ispecies)$key
## now use the (default) data from Berman, 1988
# this resets the thermodynamic database
# without affecting the basis and species settings
obigt()
# we can check that we have Berman's quartz
# and not coesite or some other phase of SiO2
iSiO2 <- rownames(basis()) == "SiO2"
stopifnot(info(basis()$ispecies[iSiO2])$name == "quartz")
# Berman's dataset doesn't have the upper temperature limits,
# so we don't need exceed.Ttr here
a <- affinity(`K+` = c(0, 5, res), T = c(200, 650, res), P = 1000, exceed.rhomin = TRUE)
diagram(a, add = TRUE, names = FALSE, col = "blue", lwd = 1.5)
# the list of references:
ref2 <- thermo.refs(species()$ispecies)$key
ref2 <- paste(ref2, collapse = ", ")
# add legend and title
legend("top", "low-density region", text.font = 3, bty = "n")
legend("topleft", describe.property(c("P", "IS"), c(1000, 1)), bty = "n")
legend("left", c(ref1, ref2),
       lty = c(1, 1), lwd = 1.5, col = c(2, 4), bty = "n")
title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O", "HCl")), line = 1.8)
title(main = "Helgeson and Berman minerals, after Sverjensky et al., 1991",
      line = 0.3, font.main = 1)
# cleanup for next example
reset()

# make a P-T diagram for SiO2 minerals (Ber88 Fig. 4)
basis(c("SiO2", "O2"), c("cr", "gas"))
species(c("quartz", "quartz,beta", "coesite"), "cr")
a <- affinity(T=c(200, 1700, 200), P=c(0, 50000, 200))
diagram(a)

## Getting data from a user-supplied file
## Ol-Opx exchange equilibrium, after Berman and Aranovich, 1996
E.units("J")
species <- c("fayalite", "enstatite", "ferrosilite", "forsterite")
coeffs <- c(-1, -2, 2, 1)
T <- seq(600, 1500, 50)
Gex_Ber88 <- subcrt(species, coeffs, T=T, P=1)$out$G
# add data from BA96
datadir <- system.file("extdata/Berman/testing", package="CHNOSZ")
add.obigt(file.path(datadir, "BA96_obigt.csv"))
thermo("opt$Berman" = file.path(datadir, "BA96_berman.csv"))
Gex_BA96 <- subcrt(species, coeffs, T=seq(600, 1500, 50), P=1)$out$G
# Ber88 is lower than BA96 at low T
stopifnot((Gex_BA96 - Gex_Ber88)[1] > 0)
# the curves cross at about 725 deg C (BA96 Fig. 8)
# (actually, in our calculation they cross closer to 800 deg C)
stopifnot(T[which.min(abs(Gex_BA96 - Gex_Ber88))] == 800)
# reset the database (thermo()$opt$E.units, thermo()$obigt, and thermo()$opt$Berman)
reset()

CHNOSZ documentation built on April 6, 2020, 3:01 a.m.