Next: buffer
Up: CHNOSZ examples
Previous: diagram: equilibrium activities of
diagrm> ### predominance diagrams (equal activities of
diagrm> ### species as a function of two variables)
diagrm>
diagrm> ## T-P stability diagram for SiO2
diagrm> ## after Ernst, 1976, Fig. 4.4
diagrm> # the system is SiO2 (one component) but
diagrm> # the basis species require all the elements:
diagrm> # note that basis(c("SiO2","O2")) would identify SiO2(aq)
diagrm> # which is okay but brings in calculations of properties
diagrm> # of water which are relatively slow.
diagrm> basis(c("quartz","O2")) # cr, gas
basis: changed basis to SiO2 O2.
O Si ispecies logact state
SiO2 2 1 2005 0 cr1
O2 2 0 2691 0 gas
diagrm> # browse the species of interest
diagrm> info("SiO2 ")
info: no match for SiO2 .
info: approximately matching species are:
name abbrv formula state
72 SiO2 SiO2 SiO2 aq
1821 amorphous silica Amor-Sl SiO2 cr
1856 chalcedony Cha SiO2 cr
1870 coesite Cos SiO2 cr1
1871 coesite Cos SiO2 cr2
1877 cristobalite Crs SiO2 cr1
1878 cristobalite Crs SiO2 cr2
1879 cristobalite,alpha A-Crs SiO2 cr
1880 cristobalite,beta B-Crs SiO2 cr
2005 quartz Qtz SiO2 cr1
2006 quartz Qtz SiO2 cr2
diagrm> # we'll take the crystalline ones
diagrm> si <- info("SiO2","cr")
info: SiO2 matches these species:
name abbrv formula state ref1
1821 amorphous silica Amor-Sl SiO2 cr HDN+78
1856 chalcedony Cha SiO2 cr HDN+78
1870 coesite Cos SiO2 cr1 HDN+78
1871 coesite Cos SiO2 cr2 HDN+78
1877 cristobalite Crs SiO2 cr1 HDN+78
1878 cristobalite Crs SiO2 cr2 HDN+78
1879 cristobalite,alpha A-Crs SiO2 cr HDN+78
1880 cristobalite,beta B-Crs SiO2 cr HDN+78
2005 quartz Qtz SiO2 cr1 HDN+78
2006 quartz Qtz SiO2 cr2 HDN+78
info: 1821 refers to amorphous silica cr (HDN+78, 5.May.78)
info: 1856 refers to chalcedony cr (HDN+78, 5.May.78)
info: 1870 refers to coesite cr1 (HDN+78, 5.May.78)
info: 1871 refers to coesite cr2 (HDN+78, 5.May.78)
info: 1877 refers to cristobalite cr1 (HDN+78, 5.May.78)
info: 1878 refers to cristobalite cr2 (HDN+78, 5.May.78)
info: 1879 refers to cristobalite,alpha cr (HDN+78, 5.May.78)
info: 1880 refers to cristobalite,beta cr (HDN+78, 5.May.78)
info: 2005 refers to quartz cr1 (HDN+78, 5.May.78)
info: 2006 refers to quartz cr2 (HDN+78, 5.May.78)
diagrm> species(si)
SiO2 O2 ispecies logact state name
1 1 0 1821 0 cr amorphous silica
2 1 0 1856 0 cr chalcedony
3 1 0 1870 0 cr1 coesite
4 1 0 1871 0 cr2 coesite
5 1 0 1877 0 cr1 cristobalite
6 1 0 1878 0 cr2 cristobalite
7 1 0 1879 0 cr cristobalite,alpha
8 1 0 1880 0 cr cristobalite,beta
9 1 0 2005 0 cr1 quartz
10 1 0 2006 0 cr2 quartz
diagrm> # calculate chemical affinities
diagrm> # the do.phases argument is necessary for
diagrm> # dealing with the phase transitions of minerals
diagrm> a <- affinity(T=c(0,2000),P=c(0,30000),do.phases=TRUE)
energy.args: variable 1 is T at 128 increments from 273.15 to 2273.15
energy.args: variable 2 is P at 128 increments from 0 to 30000
subcrt: 12 species at 16384 values of T and P
subcrt: some points below T limits for coesite cr2 (using 999999 for G).
subcrt: some points below T limits for cristobalite cr2 (using 999999 for G).
subcrt: some points below T limits for quartz cr2 (using 999999 for G).
subcrt: some points above T limits for quartz cr1 (using 999999 for G).
subcrt: some points above T limits for amorphous silica cr (using 999999 for G).
subcrt: some points above T limits for chalcedony cr (using 999999 for G).
subcrt: some points above T limits for coesite cr1 (using 999999 for G).
subcrt: some points above T limits for coesite cr2 (using 999999 for G).
subcrt: some points above T limits for cristobalite cr1 (using 999999 for G).
subcrt: some points above T limits for cristobalite cr2 (using 999999 for G).
subcrt: some points above T limits for cristobalite,alpha cr (using 999999 for G).
subcrt: some points above T limits for cristobalite,beta cr (using 999999 for G).
subcrt: some points above T limits for quartz cr1 (using 999999 for G).
subcrt: some points above T limits for quartz cr2 (using 999999 for G).
diagrm> diagram(a)
diagram: immobile component is SiO2
diagram: conservation coefficients are 1 1 1 1 1 1 1 1 1 1
diagrm> title(main="Phases of SiO2\nafter Ernst, 1976")
diagrm> ## Distribution of copper on Eh-pH diagram
diagrm> ## after Fig. 15 of Pourbaix, 1949
diagrm> basis(c("Cu+2","Cl-","H2O","H+","e-"))
basis: changed basis to Cu+2 Cl- H2O H+ e-.
Cl Cu H O Z ispecies logact state
Cu+2 0 1 0 0 2 62 0 aq
Cl- 1 0 0 0 -1 29 0 aq
H2O 0 0 2 1 0 1 0 liq
H+ 0 0 1 0 1 3 0 aq
e- 0 0 0 0 -1 2 0 aq
diagrm> basis("Cl-",-2)
Cl Cu H O Z ispecies logact state
Cu+2 0 1 0 0 2 62 0 aq
Cl- 1 0 0 0 -1 29 -2 aq
H2O 0 0 2 1 0 1 0 liq
H+ 0 0 1 0 1 3 0 aq
e- 0 0 0 0 -1 2 0 aq
diagrm> # two aqueous species, three solid ones
diagrm> # (our CuCl is aq but it was cr in P's Fig. 15)
diagrm> species(c("CuCl","Cu+2","copper","cuprite","tenorite"))
Cu+2 Cl- H2O H+ e- ispecies logact state name
1 1 1 0 0 1 1000 -3 aq CuCl
2 1 0 0 0 0 62 -3 aq Cu+2
3 1 0 0 0 2 1872 0 cr copper
4 2 0 1 -2 2 1885 0 cr cuprite
5 1 0 1 -2 0 2031 0 cr tenorite
diagrm> # (which is equivalent to ...)
diagrm> # species(c("CuCl","Cu+2","Cu","Cu2O","CuO"),c("","","","","","cr"))
diagrm> a <- affinity(pH=c(0,7),Eh=c(-0.1,1))
affinity: temperature is 25 C
energy.args: pressure is Psat
energy.args: variable 1 is pH at 128 increments from 0 to 7
energy.args: variable 2 is Eh at 128 increments from -0.1 to 1
subcrt: 10 species at 298.15 K and 1 bar (wet)
diagrm> # this one corresponds to activity contours of
diagrm> # aqueous species at 10^-3 (the default aq activity in CHNOSZ)
diagrm> diagram(a,color=NULL)
diagram: immobile component is Cu+2
diagram: conservation coefficients are 1 1 1 2 1
diagrm> # here we set activities to unity; aq-cr boundaries change
diagrm> species(c("CuCl","Cu+2"),c(0,0))
Cu+2 Cl- H2O H+ e- ispecies logact state name
1 1 1 0 0 1 1000 0 aq CuCl
2 1 0 0 0 0 62 0 aq Cu+2
3 1 0 0 0 2 1872 0 cr copper
4 2 0 1 -2 2 1885 0 cr cuprite
5 1 0 1 -2 0 2031 0 cr tenorite
diagrm> a <- affinity(pH=c(0,7),Eh=c(-0.1,1))
affinity: temperature is 25 C
energy.args: pressure is Psat
energy.args: variable 1 is pH at 128 increments from 0 to 7
energy.args: variable 2 is Eh at 128 increments from -0.1 to 1
subcrt: 10 species at 298.15 K and 1 bar (wet)
diagrm> diagram(a,add=TRUE,names=NULL,col="blue",color=NULL)
diagram: immobile component is Cu+2
diagram: conservation coefficients are 1 1 1 2 1
diagrm> water.lines()
subcrt: 3 species at 298.15 K and 1 bar (wet)
subcrt: 4 species at 2 values of T and P (wet)
subcrt: 4 species at 2 values of T and P (wet)
diagrm> title(main=paste("H2O-Cl-(Cu); activities of 10^-3 (black)\n",
diagrm+ "or 0 (blue); after Pourbaix, 1949"))
diagrm> ## a pe-pH diagram for hydrated iron sulfides,
diagrm> ## goethite and pyrite, After Majzlan et al., 2006
diagrm> # add some of these species to the database
diagrm> add.obigt()
add.obigt: from /home/jedick/R/x86_64-slackware-linux-gnu-library/2.13/CHNOSZ/extdata/thermo/OBIGT-2.csv added 301 of 301 species ( 82 replacements, 219 new, units = cal )
add.obigt: to restore default database, use data(thermo)
diagrm> 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 1038 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
diagrm> species(c("rhomboclase","ferricopiapite","hydronium jarosite",
diagrm+ "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 3156 0 cr rhomboclase
2 4.78 6.00 23.05 -2.34 -4.78 3155 0 cr ferricopiapite
3 3.00 2.00 7.00 -5.00 -3.00 3152 0 cr hydronium jarosite
4 1.00 0.00 2.00 -3.00 -1.00 3157 0 cr goethite
5 1.00 1.00 7.00 0.00 0.00 3153 0 cr melanterite
6 1.00 2.00 -8.00 16.00 14.00 1917 0 cr pyrite
diagrm> a <- affinity(pH=c(-1,4),pe=c(-5,23))
affinity: temperature is 25 C
energy.args: pressure is Psat
energy.args: variable 1 is pH at 128 increments from -1 to 4
energy.args: variable 2 is pe at 128 increments from -5 to 23
subcrt: 11 species at 298.15 K and 1 bar (wet)
diagrm> diagram(a)
diagram: immobile component is Fe+2
diagram: conservation coefficients are 1 4.78 3 1 1 1
diagrm> water.lines(yaxis="pe")
subcrt: 3 species at 298.15 K and 1 bar (wet)
subcrt: 4 species at 2 values of T and P (wet)
subcrt: 4 species at 2 values of T and P (wet)
diagrm> title(main=paste("Fe-S-O-H After Majzlan et al., 2006",
diagrm+ describe(thermo$basis[2:3,],digits=2),sep="\n"))
diagrm> # reset the database
diagrm> data(thermo)
thermo$obigt has 1800 aqueous, 2925 total species
diagrm> ## cysteine-cysteinate-cystine Eh-pH
diagrm> ## at 25 and 100 deg C
diagrm> basis("CHNOSe")
C H N O S Z ispecies logact state
CO2 1 0 0 2 0 0 69 -3 aq
H2O 0 2 0 1 0 0 1 0 liq
NH3 0 3 1 0 0 0 68 -4 aq
H2S 0 2 0 0 1 0 70 -7 aq
e- 0 0 0 0 0 -1 2 -7 aq
H+ 0 1 0 0 0 1 3 -7 aq
diagrm> species(c("cysteine","cysteinate","cystine"))
CO2 H2O NH3 H2S e- H+ ispecies logact state name
1 3 -4 1 1 10 10 1511 -3 aq cysteine
2 3 -4 1 1 10 9 1512 -3 aq cysteinate
3 6 -8 2 2 18 18 1589 -3 aq cystine
diagrm> a <- affinity(pH=c(5,10),Eh=c(-0.5,0))
affinity: temperature is 25 C
energy.args: pressure is Psat
energy.args: variable 1 is pH at 128 increments from 5 to 10
energy.args: variable 2 is Eh at 128 increments from -0.5 to 0
subcrt: 9 species at 298.15 K and 1 bar (wet)
diagrm> diagram(a,color=NULL)
diagram: immobile component is CO2
diagram: conservation coefficients are 3 3 6
diagrm> water.lines()
subcrt: 3 species at 298.15 K and 1 bar (wet)
subcrt: 4 species at 2 values of T and P (wet)
subcrt: 4 species at 2 values of T and P (wet)
diagrm> a <- affinity(pH=c(5,10),Eh=c(-0.5,0),T=100)
affinity: temperature is 100 C
energy.args: pressure is Psat
energy.args: variable 1 is pH at 128 increments from 5 to 10
energy.args: variable 2 is Eh at 128 increments from -0.5 to 0
subcrt: 9 species at 373.15 K and 1.01 bar (wet)
diagrm> diagram(a,col="red",add=TRUE,names=NULL)
diagram: immobile component is CO2
diagram: conservation coefficients are 3 3 6
diagrm> water.lines(T=convert(100,"K"),col="red")
subcrt: 3 species at 373.15 K and 1.01 bar (wet)
subcrt: 4 species at 2 values of T and P (wet)
subcrt: 4 species at 2 values of T and P (wet)
diagrm> title(main=paste("Cysteine Cysteinate Cystine",
diagrm+ "25 and 100 deg C",sep="\n"))
diagrm> ## Soil Organic Matter CO2-H2O, O2-H2O (after Tardy et al., 1997)
diagrm> # NH3 is conserved, and H2O is on an axis of the diagram
diagrm> # formulas for aqueous species, names for phases ...
diagrm> add.obigt()
add.obigt: from /home/jedick/R/x86_64-slackware-linux-gnu-library/2.13/CHNOSZ/extdata/thermo/OBIGT-2.csv added 301 of 301 species ( 82 replacements, 219 new, units = cal )
add.obigt: to restore default database, use data(thermo)
diagrm> basis(c("NH3","water","CO2","O2"),c(999,999,-2.5,-28))
C H N O ispecies logact state
NH3 0 3 1 0 68 999.0 aq
H2O 0 2 0 1 1 999.0 liq
CO2 1 0 0 2 69 -2.5 aq
O2 0 0 0 2 2609 -28.0 gas
diagrm> # switch to gaseous CO2 (aq is the default)
diagrm> basis("CO2","gas")
C H N O ispecies logact state
NH3 0 3 1 0 68 999.0 aq
H2O 0 2 0 1 1 999.0 liq
CO2 1 0 0 2 2601 -2.5 gas
O2 0 0 0 2 2609 -28.0 gas
diagrm> # load the species of interest
diagrm> species(c("microflore","plantes","acide fulvique",
diagrm+ "acide humique","humine"))
NH3 H2O CO2 O2 ispecies logact state name
1 1 7.0 7 -7.00 3129 0 cr microflore
2 1 25.5 34 -34.25 3130 0 cr plantes
3 1 18.5 30 -27.25 3131 0 cr acide fulvique
4 1 10.5 20 -19.25 3132 0 cr acide humique
5 1 43.5 72 -75.75 3133 0 cr humine
diagrm> # proceed with the diagrams
diagrm> diagram(affinity(H2O=c(-0.6,0.1),CO2=c(-3,-1)))
affinity: temperature is 25 C
energy.args: pressure is Psat
energy.args: variable 1 is H2O at 128 increments from -0.6 to 0.1
energy.args: variable 2 is CO2 at 128 increments from -3 to -1
subcrt: 9 species at 298.15 K and 1 bar (wet)
diagram: immobile component is NH3
diagram: conservation coefficients are 0.999999999999999 1 1 1 1
diagrm> title(main=paste("Relative stabilities of soil organic matter\n",
diagrm+ "after Tardy et al., 1997"))
diagrm> # this is the O2-H2O diagram
diagrm> # diagram(affinity(H2O=c(-1,0.5),O2=c(-28.5,-27.5)))
diagrm> data(thermo)
thermo$obigt has 1800 aqueous, 2925 total species
diagrm> ## Aqueous Aluminum Species F-/OH- (afer Tagirov and Schott, 2001)
diagrm> # in order to reproduce this calculation, we have to
diagrm> # consider the properties of species given by these authors,
diagrm> # which are not the default ones in the database
diagrm> add.obigt()
add.obigt: from /home/jedick/R/x86_64-slackware-linux-gnu-library/2.13/CHNOSZ/extdata/thermo/OBIGT-2.csv added 301 of 301 species ( 82 replacements, 219 new, units = cal )
add.obigt: to restore default database, use data(thermo)
diagrm> # The coefficients on H+ and O2 in all the formation reactions
diagrm> # are zero, so the number of basis species here is three. Al+3
diagrm> # becomes the conservant, and F- and OH- are being plotted ...
diagrm> # so their initial activities don't have to be set.
diagrm> basis(c("Al+3","F-","OH-","H+","O2"),rep(999,5))
Al F H O Z ispecies logact state
Al+3 1 0 0 0 3 2855 999 aq
F- 0 1 0 0 -1 28 999 aq
OH- 0 0 1 1 -1 21 999 aq
H+ 0 0 1 0 1 3 999 aq
O2 0 0 0 2 0 2609 999 gas
diagrm> species(c("Al+3","Al(OH)4-","AlOH+2","Al(OH)2+","Al(OH)3",
diagrm+ "AlF+2","AlF2+","AlF3","AlF4-","Al(OH)2F2-","Al(OH)2F","AlOHF2"))
Al+3 F- OH- H+ O2 ispecies logact state name
1 1 0 0 0 0 2855 -3 aq Al+3
2 1 0 4 0 0 2856 -3 aq Al(OH)4-
3 1 0 1 0 0 2857 -3 aq AlOH+2
4 1 0 2 0 0 2858 -3 aq Al(OH)2+
5 1 0 3 0 0 2859 -3 aq Al(OH)3
6 1 1 0 0 0 2860 -3 aq AlF+2
7 1 2 0 0 0 2861 -3 aq AlF2+
8 1 3 0 0 0 2862 -3 aq AlF3
9 1 4 0 0 0 2863 -3 aq AlF4-
10 1 2 2 0 0 2864 -3 aq Al(OH)2F2-
11 1 1 2 0 0 2865 -3 aq Al(OH)2F
12 1 2 1 0 0 2866 -3 aq AlOHF2
diagrm> # Increase the x- and y- resolution from the default and calculate
diagrm> # and plot predominance limits. Names of charged basis species,
diagrm> # such as "H+", "e-" and the ones shown here, should be quoted
diagrm> # when given as arguments to affinity(). The OH- values shown here
diagrm> # correspond to pH=c(0,14) (at unit activity of water).
diagrm> a <- affinity("OH-"=c(-14,0),"F-"=c(-1,-8),T=200)
affinity: temperature is 200 C
energy.args: pressure is Psat
energy.args: variable 1 is OH- at 128 increments from -14 to 0
energy.args: variable 2 is F- at 128 increments from -1 to -8
subcrt: 17 species at 473.15 K and 15.54 bar (wet)
diagrm> diagram(a)
diagram: immobile component is Al+3
diagram: conservation coefficients are 1 1 1 1 1 1 1 1 1 1 1 1
diagrm> title(main=paste("Aqueous aluminium species, T=200 C, P=Psat\n",
diagrm+ "after Tagirov and Schott, 2001"))
diagrm> # We could do this to overlay boundaries for a different pressure
diagrm> #a.P <- affinity("OH-"=c(-14,0),"F-"=c(-1,-8),T=200,P=5000)
diagrm> #diagram(a.P,names=NULL,color=NULL,col="blue",add=TRUE)
diagrm> # restore thermodynamic database to default
diagrm> data(thermo)
thermo$obigt has 1800 aqueous, 2925 total species
diagrm> ## Fe-S-O at 200 deg C, After Helgeson, 1970
diagrm> basis(c("Fe","O2","S2"))
Fe O S ispecies logact state
Fe 1 0 0 2315 0 cr1
O2 0 2 0 2691 0 gas
S2 0 0 2 2693 0 gas
diagrm> species(c("iron","ferrous-oxide","magnetite",
diagrm+ "hematite","pyrite","pyrrhotite"))
Fe O2 S2 ispecies logact state name
1 1 0.0 0.0 2315 0 cr1 Fe
2 1 0.5 0.0 1914 0 cr ferrous-oxide
3 3 2.0 0.0 1966 0 cr1 magnetite
4 2 1.5 0.0 1936 0 cr1 hematite
5 1 0.0 1.0 1998 0 cr pyrite
6 1 0.0 0.5 2002 0 cr1 pyrrhotite
diagrm> a <- affinity(S2=c(-50,0),O2=c(-90,-10),T=200)
affinity: temperature is 200 C
energy.args: pressure is Psat
energy.args: variable 1 is S2 at 128 increments from -50 to 0
energy.args: variable 2 is O2 at 128 increments from -90 to -10
subcrt: 9 species at 473.15 K and 15.54 bar
subcrt: some points above T limits for pyrrhotite cr1 (ignored).
diagrm> diagram(a,color="heat")
diagram: immobile component is Fe
diagram: conservation coefficients are 1 1 3 2 1 1
diagrm> title(main=paste("Fe-S-O, 200 degrees C, 1 bar",
diagrm+ "After Helgeson, 1970",sep="\n"))
diagrm> ## Temperature-Pressure: kayanite-sillimanite-andalusite
diagrm> # this is a system of a single component (Al2SiO5)
diagrm> # but we still have to use the same number of
diagrm> # basis species as the number of elements
diagrm> basis(c("kyanite","quartz","oxygen"))
basis: changed basis to Al2SiO5 SiO2 O2.
Al O Si ispecies logact state
Al2SiO5 2 5 1 1950 0 cr
SiO2 0 2 1 2005 0 cr1
O2 0 2 0 2691 0 gas
diagrm> species(c("kyanite","sillimanite","andalusite"))
Al2SiO5 SiO2 O2 ispecies logact state name
1 1 0 0 1950 0 cr kyanite
2 1 0 0 2017 0 cr sillimanite
3 1 0 0 1824 0 cr andalusite
diagrm> diagram(affinity(T=c(0,1000,200),P=c(500,5000,200)),color=NULL)
energy.args: variable 1 is T at 200 increments from 273.15 to 1273.15
energy.args: variable 2 is P at 200 increments from 500 to 5000
subcrt: 6 species at 40000 values of T and P
subcrt: some points above T limits for kyanite cr (ignored).
subcrt: some points above T limits for quartz cr1 (ignored).
subcrt: some points above T limits for kyanite cr (ignored).
subcrt: some points above T limits for sillimanite cr (ignored).
subcrt: some points above T limits for andalusite cr (ignored).
diagram: immobile component is Al2SiO5
diagram: conservation coefficients are 1 1 1
diagrm> title(main="Al2SiO5")
diagrm> ## End(No test) # end donttest
diagrm>
diagrm>
diagrm>
Next: buffer
Up: CHNOSZ examples
Previous: diagram: equilibrium activities of
