next up previous
Next: [1] diagram (other plots) Up: CHNOSZ examples Previous: [3] diagram (1 variable)

[4] diagram (2 variables)

### 2-D diagrams: predominance diagrams ### these use the maximum affinity method ## Fe-S-O at 200 deg C, after Helgeson, 1970 basis(c("Fe", "O2", "S2"))
Fe O S ispecies logact state Fe 1 0 0 2535 0 cr1 O2 0 2 0 67 0 aq S2 0 0 2 3097 0 gas
species(c("iron", "ferrous-oxide", "magnetite", "hematite", "pyrite", "pyrrhotite"))
Fe O2 S2 ispecies logact state name 1 1 0.0 0.0 2535 0 cr1 Fe 2 1 0.5 0.0 1923 0 cr ferrous-oxide 3 3 2.0 0.0 1975 0 cr1 magnetite 4 2 1.5 0.0 1945 0 cr1 hematite 5 1 0.0 1.0 2007 0 cr pyrite 6 1 0.0 0.5 2011 0 cr1 pyrrhotite
# include the high-temperature phase of pyrrhotite species("pyrrhotite", "cr2")
Fe O2 S2 ispecies logact state name 1 1 0.0 0.0 2535 0 cr1 Fe 2 1 0.5 0.0 1923 0 cr ferrous-oxide 3 3 2.0 0.0 1975 0 cr1 magnetite 4 2 1.5 0.0 1945 0 cr1 hematite 5 1 0.0 1.0 2007 0 cr pyrite 6 1 0.0 0.5 2011 0 cr1 pyrrhotite 7 1 0.0 0.5 2012 0 cr2 pyrrhotite
a <- affinity(S2=c(-50, 0), O2=c(-90, -10), T=200)
energy.args: temperature is 200 C energy.args: pressure is Psat energy.args: variable 1 is log_f(S2) at 128 values from -50 to 0 energy.args: variable 2 is log_a(O2) at 128 values from -90 to -10 subcrt: 10 species at 473.15 K and 15.54 bar (wet) subcrt: some points above transition temperature for pyrrhotite cr1 (using NA for G)
diagram(a, fill="heat")
balance: coefficients are moles of Fe in formation reactions diagram: plotting A/2.303RT from affinity(), divided by balancing coefficients
title(main=paste("Fe-S-O, 200 degrees C, 1 bar", "After Helgeson, 1970", sep="\n"))

Image diagram4

 

## pe-pH diagram for hydrated iron sulfides, ## goethite and pyrite, after Majzlan et al., 2006 # add some of these species to the database add.obigt()
add.obigt: using default file: /home/jedick/R/x86_64-slackware-linux-gnu-library/3.3/CHNOSZ/extdata/thermo/OBIGT-2.csv add.obigt: read 305 rows; made 84 replacements, 221 additions, units = cal add.obigt: use data(thermo) to restore default database
basis(c("Fe+2", "SO4-2", "H2O", "H+", "e-"), c(0, log10(3), log10(0.75), 999, 999))
Fe H O S Z ispecies logact state Fe+2 1 0 0 0 2 1109 0.0000000 aq SO4-2 0 0 4 1 -2 24 0.4771213 aq H2O 0 2 1 0 0 1 -0.1249387 liq H+ 0 1 0 0 1 3 999.0000000 aq e- 0 0 0 0 -1 2 999.0000000 aq
species(c("rhomboclase", "ferricopiapite", "hydronium jarosite", "goethite", "melanterite", "pyrite"))
Fe+2 SO4-2 H2O H+ e- ispecies logact state name 1 1.00 2.17 4.40 1.34 -1.00 3587 0 cr rhomboclase 2 4.78 6.00 23.05 -2.34 -4.78 3586 0 cr ferricopiapite 3 3.00 2.00 7.00 -5.00 -3.00 3583 0 cr hydronium jarosite 4 1.00 0.00 2.00 -3.00 -1.00 2539 0 cr goethite 5 1.00 1.00 7.00 0.00 0.00 3584 0 cr melanterite 6 1.00 2.00 -8.00 16.00 14.00 2007 0 cr pyrite
a <- affinity(pH=c(-1, 4, 256), pe=c(-5, 23, 256))
energy.args: temperature is 25 C energy.args: pressure is Psat energy.args: variable 1 is pH at 256 values from -1 to 4 energy.args: variable 2 is pe at 256 values from -5 to 23 subcrt: 11 species at 298.15 K and 1 bar (wet)
d <- diagram(a, main="Fe-S-O-H, after Majzlan et al., 2006")
balance: coefficients are moles of Fe+2 in formation reactions diagram: plotting A/2.303RT from affinity(), divided by balancing coefficients
# the first four species show up along the top of the diagram stopifnot(all.equal(unique(t(d$predominant)[256,]), 1:4)) water.lines(yaxis="pe")
subcrt: 3 species at 298.15 K and 1 bar (wet) info.character: found H2O(liq), also available in gas, cr1, cr2, cr3, cr5, cr6, cr7, cr8, cr9 subcrt: 4 species at 298.15 K and 1 bar (wet) subcrt: 2 species at 298.15 K and 1 bar info.character: found H2O(liq), also available in gas, cr1, cr2, cr3, cr5, cr6, cr7, cr8, cr9 subcrt: 4 species at 298.15 K and 1 bar (wet)
text(3, 22, describe.basis(thermo$basis[2:3,], digits=2, oneline=TRUE)) text(3, 21, describe.property(c("T", "P"), c(25, 1), oneline=TRUE)) # reset the database data(thermo)
thermo$obigt: 1809 aqueous, 3368 total species

Image diagram5

 

## Aqueous Aluminum Species F-/OH-, after Tagirov and Schott, 2001 # some of the species are not in the default databse add.obigt()
add.obigt: using default file: /home/jedick/R/x86_64-slackware-linux-gnu-library/3.3/CHNOSZ/extdata/thermo/OBIGT-2.csv add.obigt: read 305 rows; made 84 replacements, 221 additions, units = cal add.obigt: use data(thermo) to restore default database
# the 999s have no effect on the diagram: # pH and log_a(F-) are plotting variables # O2 is not in the formation reactions # Al+3 is the balanced quantity basis(c("Al+3", "F-", "H+", "O2", "H2O"), c(rep(999, 4), 0))
Al F H O Z ispecies logact state Al+3 1 0 0 0 3 1019 999 aq F- 0 1 0 0 -1 28 999 aq H+ 0 0 1 0 1 3 999 aq O2 0 0 0 2 0 67 999 aq H2O 0 0 2 1 0 1 0 liq
species(c("Al+3", "Al(OH)4-", "AlOH+2", "Al(OH)2+", "Al(OH)3", "AlF+2", "AlF2+", "AlF3", "AlF4-", "Al(OH)2F2-", "Al(OH)2F", "AlOHF2"), "aq")
Al+3 F- H+ O2 H2O ispecies logact state name 1 1 0 0 0 0 1019 -3 aq Al+3 2 1 0 -4 0 4 3376 -3 aq Al(OH)4- 3 1 0 -1 0 1 1179 -3 aq AlOH+2 4 1 0 -2 0 2 3377 -3 aq Al(OH)2+ 5 1 0 -3 0 3 3378 -3 aq Al(OH)3 6 1 1 0 0 0 3379 -3 aq AlF+2 7 1 2 0 0 0 3380 -3 aq AlF2+ 8 1 3 0 0 0 3381 -3 aq AlF3 9 1 4 0 0 0 3382 -3 aq AlF4- 10 1 2 -2 0 2 3383 -3 aq Al(OH)2F2- 11 1 1 -2 0 2 3384 -3 aq Al(OH)2F 12 1 2 -1 0 1 3385 -3 aq AlOHF2
a <- affinity(pH=c(0, 10), "F-"=c(-1, -9), T=200)
energy.args: temperature is 200 C energy.args: pressure is Psat energy.args: variable 1 is pH at 128 values from 0 to 10 energy.args: variable 2 is log_a(F-) at 128 values from -1 to -9 subcrt: 17 species at 473.15 K and 15.54 bar (wet)
diagram(a, fill="heat")
balance: coefficients are moles of Al+3 in formation reactions diagram: plotting A/2.303RT from affinity(), divided by balancing coefficients
title(main=paste("Aqueous aluminium species, T=200 C, P=Psat\n", "after Tagirov and Schott, 2001")) # restore thermodynamic database to default data(thermo)
thermo$obigt: 1809 aqueous, 3368 total species

Image diagram6

 

## Temperature-Pressure: kayanite-sillimanite-andalusite # cf. Fig. 49 of Helgeson et al., 1978 # this is a system of one component (Al2SiO5), but we need the same # number of basis species as elements; and avoid using H2O or other # aqueous species because of T/P limits of the water() calculations; basis(c("corundum", "quartz", "oxygen"))
Al O Si ispecies logact state Al2O3 2 3 0 1884 0 cr SiO2 0 2 1 2014 0 cr1 O2 0 2 0 3095 0 gas
species(c("kyanite", "sillimanite", "andalusite"))
Al2O3 SiO2 O2 ispecies logact state name 1 1 1 0 1959 0 cr kyanite 2 1 1 0 2026 0 cr sillimanite 3 1 1 0 1833 0 cr andalusite
# database has transition temperatures of kyanite and andalusite # at 1 bar only, so we permit calculation at higher temperatures a <- affinity(T=c(200, 900, 99), P=c(0, 9000, 101), exceed.Ttr=TRUE)
energy.args: variable 1 is T at 99 values from 473.15 to 1173.15 K energy.args: variable 2 is P at 101 values from 0 to 9000 bar subcrt: 6 species at 9999 values of T and P subcrt: some points above transition temperature for quartz cr1 (extrapolating G) subcrt: some points above transition temperature for kyanite cr (extrapolating G) subcrt: some points above transition temperature for andalusite cr (extrapolating G)
d <- diagram(a, fill=NULL)
balance: coefficients are moles of Al2O3 in formation reactions diagram: plotting A/2.303RT from affinity(), divided by balancing coefficients
bexpr <- sapply(c("Al2O3", "SiO2", "H2O"), expr.species, simplify=FALSE) btext <- substitute(Al2O3 - SiO2 - H2O, unlist(bexpr)) mtitle(c(as.expression(btext), "after Helgeson et al., 1978")) # find the approximate position of the triple point tp <- find.tp(d$predominant) Ttp <- a$vals[[1]][tp[1, 2]] Ptp <- rev(a$vals[[2]])[tp[1, 1]] points(Ttp, Ptp, pch=10, cex=5)

Image diagram7

 

# some testing of the overall geometry stopifnot(species()$name[d$predominant[1, 1]]=="andalusite") stopifnot(species()$name[d$predominant[1, 101]]=="kyanite") stopifnot(species()$name[d$predominant[99, 101]]=="sillimanite")


next up previous
Next: [1] diagram (other plots) Up: CHNOSZ examples Previous: [3] diagram (1 variable)