Next: buffer
Up: CHNOSZ examples
Previous: diagram: 1-D plots
diagrm> ## phosphate predominance f(IS,pH)
diagrm> diagram(affinity(pH=c(6.8,7.2),IS=c(0,0.25),T=T[1]))
affinity: temperature is 25 C
energy.args: pressure is Psat
energy.args: variable 1 is pH at 128 increments from 6.8 to 7.2
energy.args: variable 2 is IS at 128 increments from 0 to 0.25
subcrt: 7 species at 298.15 K and 1 bar (wet)
subcrt: 2 species at 298.15 K and P (wet)
nonideal: 2 species were nonideal.
diagram: immobile component is H3PO4
diagram: conservation coefficients are 1 1
diagrm> title(main=paste("Aqueous phosphate predominance, 25 deg C",
diagrm+ "and 100 deg C (dotted overlay)",sep="\n"))
diagrm> diagram(affinity(pH=c(6.8,7.2),IS=c(0,0.25),T=T[2]),add=TRUE,dotted=2,
diagrm+ names=NULL)
affinity: temperature is 100 C
energy.args: pressure is Psat
energy.args: variable 1 is pH at 128 increments from 6.8 to 7.2
energy.args: variable 2 is IS at 128 increments from 0 to 0.25
subcrt: 7 species at 373.15 K and 1.01 bar (wet)
subcrt: 2 species at 373.15 K and P (wet)
nonideal: 2 species were nonideal.
diagram: immobile component is H3PO4
diagram: conservation coefficients are 1 1
diagrm> ### predominance diagrams (equal activities of
diagrm> ### species as a function of two variables)
diagrm>
diagrm> ## T-P stability diagram for a unary system (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 2202 0 cr1
O2 2 0 2852 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
2018 amorphous silica Amor-Sl SiO2 cr
2053 chalcedony Cha SiO2 cr
2067 coesite Cos SiO2 cr1
2068 coesite Cos SiO2 cr2
2074 cristobalite Crs SiO2 cr1
2075 cristobalite Crs SiO2 cr2
2076 cristobalite,alpha A-Crs SiO2 cr
2077 cristobalite,beta B-Crs SiO2 cr
2202 quartz Qtz SiO2 cr1
2203 quartz Qtz SiO2 cr2
diagrm> # we'll take the crystalline ones
diagrm> t <- info("SiO2","cr")
info: SiO2 matches these species:
name abbrv formula state source1
2018 amorphous silica Amor-Sl SiO2 cr HDN+78
2053 chalcedony Cha SiO2 cr HDN+78
2067 coesite Cos SiO2 cr1 HDN+78
2068 coesite Cos SiO2 cr2 HDN+78
2074 cristobalite Crs SiO2 cr1 HDN+78
2075 cristobalite Crs SiO2 cr2 HDN+78
2076 cristobalite,alpha A-Crs SiO2 cr HDN+78
2077 cristobalite,beta B-Crs SiO2 cr HDN+78
2202 quartz Qtz SiO2 cr1 HDN+78
2203 quartz Qtz SiO2 cr2 HDN+78
info: 2018 refers to amorphous silica cr (HDN+78, 5.May.78).
info: 2053 refers to chalcedony cr (HDN+78, 5.May.78).
info: 2067 refers to coesite cr1 (HDN+78, 5.May.78).
info: 2068 refers to coesite cr2 (HDN+78, 5.May.78).
info: 2074 refers to cristobalite cr1 (HDN+78, 5.May.78).
info: 2075 refers to cristobalite cr2 (HDN+78, 5.May.78).
info: 2076 refers to cristobalite,alpha cr (HDN+78, 5.May.78).
info: 2077 refers to cristobalite,beta cr (HDN+78, 5.May.78).
info: 2202 refers to quartz cr1 (HDN+78, 5.May.78).
info: 2203 refers to quartz cr2 (HDN+78, 5.May.78).
diagrm> species(t)
diagrm> # calculate chemical affinities
diagrm> # the do.phases argument is necessary for
diagrm> # dealing with the phase transitions of minerals
diagrm> t <- 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(t)
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"))
diagrm> # (which is equivalent to ...)
diagrm> # species(c("CuCl","Cu+2","Cu","Cu2O","CuO"),c("","","","","","cr"))
diagrm> t <- 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(t,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))
diagrm> t <- 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(t,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> # we could do it at higher temperatures too
diagrm> # t <- affinity(pH=c(0,7),Eh=c(-0.1,1),T=100)
diagrm> # diagram(t,add=TRUE,dotted=2,col="red")
diagrm>
diagrm> ## Eh-pH diagrams for copper-water-glycine
diagrm> ## After Fig. 2 of Aksu and Doyle, 2001
diagrm> ## We need to add some species and change some Gibbs energies.
diagrm> # copy two rows of the data object (actually one row twice)
diagrm> # and change them as needed
diagrm> i <- info(c("Cu(Gly)+","glycinate","e-","H+"))
info: 796 refers to Cu(Gly)+, Cu(C2H4NO2)+ aq (SK95, AH97b, 3.Sep.06).
info: 563 refers to glycinate, C2H4NO2- aq (SK95, AH97a, 3.Sep.06).
info: 2 refers to e-, Z-1 aq (CHNOSZ, 28.Oct.06).
info: 3 refers to H+ aq (SH88, 06.Nov.97).
diagrm> n <- nrow(thermo$obigt <- rbind(thermo$obigt,thermo$obigt[rep(i[1],2),]))
diagrm> thermo$obigt$name[n-1] <- "Cu(Gly)2-"
diagrm> thermo$obigt$name[n] <- "HCu(Gly)+2"
diagrm> thermo$obigt$formula[n-1] <- makeup(makeup(c(i[1],i[2],i[3]),""),"")
diagrm> thermo$obigt$formula[n] <- makeup(makeup(c(i[1],i[4]),""),"")
diagrm> # In Fig 2b, total log activities of Cu (Cu_T)
diagrm> # and glycine (L_T) are -4 and -1
diagrm> basis(c("Cu+2","H2O","H+","e-","glycine","CO2"),c(999,0,999,999,-1,999))
C Cu H N O Z ispecies logact state
Cu+2 0 1 0 0 0 2 62 999 aq
H2O 0 0 2 0 1 0 1 0 liq
H+ 0 0 1 0 0 1 3 999 aq
e- 0 0 0 0 0 -1 2 999 aq
C2H5NO2 2 0 5 1 2 0 1630 -1 aq
CO2 1 0 0 0 2 0 69 999 aq
diagrm> # solid species
diagrm> species(c("copper","cuprite","tenorite"))
diagrm> # aqueous species
diagrm> species(c("glycinium","glycine","glycinate","Cu+","Cu+2","CuO2-2","HCuO2-",
diagrm+ "Cu(Gly)+","Cu(Gly)2","Cu(Gly)2-","HCu(Gly)+2"),-4)
diagrm> ispecies <- species()$ispecies
diagrm> # update the Gibbs energies using A&D's Table 1 and Table II
diagrm> logK <- c(convert(convert(c(0,-146,-129.7,-384.061,-370.647,-314.833,
diagrm+ 49.98,65.49,-183.6,-258.5,-298.2)*1000,"cal"),"logK"),15.64,10.1,2.92)
diagrm> # do it in order so later species take account of prev. species' values
diagrm> for(i in 1:length(logK)) {
diagrm+ G <- convert(logK[i],"G")
diagrm+ if(i==12) G <- G + thermo$obigt$G[ispecies[8]] +
diagrm+ 2*thermo$obigt$G[ispecies[6]]
diagrm+ if(i==13) G <- G + thermo$obigt$G[ispecies[7]] +
diagrm+ 2*thermo$obigt$G[ispecies[6]]
diagrm+ if(i==14) G <- G + thermo$obigt$G[ispecies[11]]
diagrm+ thermo$obigt$G[ispecies[i]] <- G
diagrm+ } # done with changing Gibbs free energies!
diagrm> # we have to get some leftovers out of there or diagram() gets confused
diagrm> species(c("glycinium","glycine","glycinate"),delete=TRUE)
diagrm> # make a plot to see if it's working
diagrm> ispecies <- ispecies[-(1:6)]
diagrm> afun <- function(cu,gly) {
diagrm+ # from above: our fifth basis species is glycine(-ate,-ium)
diagrm+ basis(rownames(basis())[5],gly)
diagrm+ t <- match(ispecies,species()$ispecies)
diagrm+ species(t,cu)
diagrm+ affinity(pH=c(0,16),Eh=c(-0.6,1.0))
diagrm+ }
diagrm> diagram(afun(-4,-1))
affinity: temperature is 25 C
energy.args: pressure is Psat
energy.args: variable 1 is pH at 128 increments from 0 to 16
energy.args: variable 2 is Eh at 128 increments from -0.6 to 1
subcrt: 17 species at 298.15 K and 1 bar (wet)
diagram: immobile component is Cu+2
diagram: conservation coefficients are 1 2 1 1 1 1 1 1 1 1 1
diagrm> title(main=paste("Aqueous Copper + Glycine, 25 deg C, 1 bar",
diagrm+ "After Aksu and Doyle, 2001 Fig. 2b",sep="\n"))
diagrm> # What's missing? Try glycinate not glycine in reactions at ph > ~9.8
diagrm> basis(c("Cu+2","H2O","H+","e-","glycinate","CO2"),
diagrm+ c(999,0,999,999,-2,999))
C Cu H N O Z ispecies logact state
Cu+2 0 1 0 0 0 2 62 999 aq
H2O 0 0 2 0 1 0 1 0 liq
H+ 0 0 1 0 0 1 3 999 aq
e- 0 0 0 0 0 -1 2 999 aq
C2H4NO2- 2 0 4 1 2 -1 563 -2 aq
CO2 1 0 0 0 2 0 69 999 aq
diagrm> species(c("copper","cuprite","tenorite","Cu+","Cu+2","CuO2-2","HCuO2-",
diagrm+ "Cu(Gly)+","Cu(Gly)2","Cu(Gly)2-","HCu(Gly)+2"))
diagrm> loga_Cu <- -4
diagrm> loga_Gly <- -1
diagrm> diagram(afun(loga_Cu,loga_Gly),color=NULL,col="blue",
diagrm+ names=species()$name,col.names="blue",add=TRUE)
affinity: temperature is 25 C
energy.args: pressure is Psat
energy.args: variable 1 is pH at 128 increments from 0 to 16
energy.args: variable 2 is Eh at 128 increments from -0.6 to 1
subcrt: 17 species at 298.15 K and 1 bar (wet)
diagram: immobile component is Cu+2
diagram: conservation coefficients are 1 2 1 1 1 1 1 1 1 1 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> # the glycine ionization constants could be calculated using
diagrm> # subcrt, here they are taken from A&D Table II
diagrm> abline(v=c(2.35,9.778),lty=3)
diagrm> # now put glycinium (low pH) in the basis
diagrm> basis(c("Cu+2","H2O","H+","e-","glycinium","CO2"),c(999,0,999,999,-2,999))
C Cu H N O Z ispecies logact state
Cu+2 0 1 0 0 0 2 62 999 aq
H2O 0 0 2 0 1 0 1 0 liq
H+ 0 0 1 0 0 1 3 999 aq
e- 0 0 0 0 0 -1 2 999 aq
C2H6NO2+ 2 0 6 1 2 1 1631 -2 aq
CO2 1 0 0 0 2 0 69 999 aq
diagrm> species(c("copper","cuprite","tenorite","Cu+","Cu+2","CuO2-2","HCuO2-",
diagrm+ "Cu(Gly)+","Cu(Gly)2","Cu(Gly)2-","HCu(Gly)+2"))
diagrm> diagram(afun(loga_Cu,loga_Gly),color=NULL,col="green",
diagrm+ names=NULL,col.names="green",add=TRUE)
affinity: temperature is 25 C
energy.args: pressure is Psat
energy.args: variable 1 is pH at 128 increments from 0 to 16
energy.args: variable 2 is Eh at 128 increments from -0.6 to 1
subcrt: 17 species at 298.15 K and 1 bar (wet)
diagram: immobile component is Cu+2
diagram: conservation coefficients are 1 2 1 1 1 1 1 1 1 1 1
diagrm> # let's forget our changes to 'thermo' so the example
diagrm> # below that uses glycine will work as expected
diagrm> data(thermo)
thermo: loaded 1997 aqueous, 3089 total species to thermo$obigt
thermo: loaded 5264 proteins to thermo$ECO
thermo: loaded 6717 proteins to thermo$SGD
thermo: loaded 4155 localizations and 3570 abundances to thermo$yeastgfp
diagrm> ## a pe-pH diagram for hydrated iron sulfides,
diagrm> ## goethite and pyrite, After Majzlan et al., 2006
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 1080 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"))
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)
info: a,b,c of rhomboclase cr are NA; set to 0.
info: a,b,c of ferricopiapite cr are NA; set to 0.
info: V of hydronium jarosite cr are NA; set to 0.
info: V,a,b,c of melanterite cr are NA; set to 0.
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> ## 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"))
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
diagrm> # diagrams, so their activities are nonsensical here.
diagrm> # formulas for aqueous species, names for phases ...
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 2852 -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 2844 -2.5 gas
O2 0 0 0 2 2852 -28.0 gas
diagrm> # load the species of interest
diagrm> species(c("microflore","plantes","acide fulvique",
diagrm+ "acide humique","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)
info: V,a,b,c,d,e,f of microflore cr are NA; set to 0.
info: V,a,b,c,d,e,f of plantes cr are NA; set to 0.
info: V,a,b,c,d,e,f of acide fulvique cr are NA; set to 0.
info: V,a,b,c,d,e,f of acide humique cr are NA; set to 0.
info: V,a,b,c,d,e,f of humine cr are NA; set to 0.
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 would be the O2-H2O diagram
diagrm> # diagram(affinity(H2O=c(-1,0.5),O2=c(-28.5,-27.5)))
diagrm>
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> # (see last example of 'thermo' help page)
diagrm> thermo$opt$level <- 2
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 1588 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 2852 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"))
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 data selection level to default
diagrm> thermo$opt$level <- 1
diagrm> ## Mineral equilibrium activity diagram
diagrm> ## (After Bowers et al., 1984, p. 246)
diagrm> # Consider the system denoted by BJH84 as
diagrm> # HCl-H2O-CaO-CO2-MgO-(SiO2) at 300 deg C and 1000 bar
diagrm> # HCl, CO2, O2 have zero stoichiometric coeffs in the species
diagrm> # Ca+2, Mg+2 are going to be on the diagram
diagrm> # SiO2 will be conserved (the immobile component)
diagrm> # (try changing any of the 999's, the diagram will be unaffected)
diagrm> basis(c("HCl","H2O","Ca+2","CO2","Mg+2","SiO2","O2","H+"),
diagrm+ c(999,0,999,999,999,999,999,-7))
C Ca Cl H Mg O Si Z ispecies logact state
HCl 0 0 1 1 0 0 0 0 854 999 aq
H2O 0 0 0 2 0 1 0 0 1 0 liq
Ca+2 0 1 0 0 0 0 0 2 10 999 aq
CO2 1 0 0 0 0 2 0 0 69 999 aq
Mg+2 0 0 0 0 1 0 0 2 9 999 aq
SiO2 0 0 0 0 0 2 1 0 72 999 aq
O2 0 0 0 0 0 2 0 0 2852 999 gas
H+ 0 0 0 1 0 0 0 1 3 -7 aq
diagrm> # we have to tell CHNOSZ which species to include; there are
diagrm> # others that could be in this system
diagrm> species(c("quartz","talc","forsterite","tremolite","diopside",
diagrm+ "wollastonite","monticellite","merwinite"))
diagrm> # calculate the chemical affinities of formation reactions
diagrm> #t <- affinity("Mg+2"=c(-14,-8),"Ca+2"=c(-12,-4),T=300,P=1000)
diagrm> t <- affinity("Mg+2"=c(-12,-4),"Ca+2"=c(-8,0),T=300,P=1000)
affinity: temperature is 300 C
affinity: pressure is 1000 bar
energy.args: variable 1 is Mg+2 at 128 increments from -12 to -4
energy.args: variable 2 is Ca+2 at 128 increments from -8 to 0
subcrt: 16 species at 573.15 K and 1000 bar (wet)
diagrm> diagram(t)
diagram: immobile component is SiO2
diagram: conservation coefficients are 1 4 1 8 2 1 1 2
diagrm> title(main=paste("HCl-H2O-CaO-CO2-MgO-(SiO2) at 300 deg C,\n",
diagrm+ "1000 bar and pH=7. After Bowers et al., 1984"))
diagrm> # note: BJH84 use a different method for representing
diagrm> # the axes of the diagrams, similar to (a_Ca+2)/(a_H+)^2,
diagrm> # so this is so far an approximate reproduction of their diagram.
diagrm>
diagrm> ## Fe-S-O at 200 deg C
diagrm> ## (After Helgeson, 1970)
diagrm> basis(c("Fe","O2","S2"))
basis: changed basis to Fe O2 S2.
Fe O S ispecies logact state
Fe 1 0 0 2507 0 cr1
O2 0 2 0 2852 0 gas
S2 0 0 2 2854 0 gas
diagrm> species(c("iron","ferrous-oxide","magnetite",
diagrm+ "hematite","pyrite","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> ## Nucleobase - Amino Acid Interaction Eh-H2O
diagrm> ## association of amino acids with
diagrm> ## nucleotides in the genetic code.
diagrm> # for this example we try a unique basis definition
diagrm> basis(c("CO2","H2O","glutamine","e-","H+"),c(-3,0,-3,0,-7))
C H N O Z ispecies logact state
CO2 1 0 0 2 0 69 -3 aq
H2O 0 2 0 1 0 1 0 liq
C5H10N2O3 5 10 2 3 0 1627 -3 aq
e- 0 0 0 0 -1 2 0 aq
H+ 0 1 0 0 1 3 -7 aq
diagrm> species(c("uracil","cytosine","adenine","guanine",
diagrm+ "phenylalanine","proline","lysine","glycine"),"aq")
diagrm> # this loaded four nucleobases and four related amino acids
diagrm> # (coded for by the homocodon triplets)
diagrm> # check out the predominance diagrams
diagrm> a.1 <- affinity(H2O=c(-5,0),Eh=c(-0.5,0))
affinity: temperature is 25 C
energy.args: pressure is Psat
energy.args: variable 1 is H2O at 128 increments from -5 to 0
energy.args: variable 2 is Eh at 128 increments from -0.5 to 0
subcrt: 13 species at 298.15 K and 1 bar (wet)
diagrm> diagram(a.1,color=NULL)
diagram: immobile component is CO2
diagram: conservation coefficients are -1 -3.5 -7.5 -7.5 6.5 2.5 1 -0.5
diagrm> # overlay a different temperature
diagrm> a.2 <- affinity(H2O=c(-5,0),Eh=c(-0.5,0),T=100)
affinity: temperature is 100 C
energy.args: pressure is Psat
energy.args: variable 1 is H2O at 128 increments from -5 to 0
energy.args: variable 2 is Eh at 128 increments from -0.5 to 0
subcrt: 13 species at 373.15 K and 1.01 bar (wet)
diagrm> diagram(a.2,col="red",add=TRUE,names=NULL)
diagram: immobile component is CO2
diagram: conservation coefficients are -1 -3.5 -7.5 -7.5 6.5 2.5 1 -0.5
diagrm> # start make a title for the plot
diagrm> tb <- thermo$basis # includes activities of basis species
diagrm> # exclude those that are on the axes
diagrm> tb <- tb[!((rownames(tb) %in% c("e-","H2O"))),]
diagrm> title(main=paste("Nucleobases and amino acids,",
diagrm+ "T=25 and 100 C at Psat\n",
diagrm+ describe(tb,T=NULL,P=NULL)),cex.main=0.9)
diagrm> ## Temperature-Pressure: kayanite-sillimanite-andalusite
diagrm> # this is a system of a single component (Al2SiO5)
diagrm> # but we are required to provide put more than one
diagrm> # species in the basis definition
diagrm> basis(c("kyanite","quartz","oxygen"))
basis: changed basis to Al2SiO5 SiO2 O2.
Al O Si ispecies logact state
Al2SiO5 2 5 1 2147 0 cr
SiO2 0 2 1 2202 0 cr1
O2 0 2 0 2852 0 gas
diagrm> species(c("kyanite","sillimanite","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")
Next: buffer
Up: CHNOSZ examples
Previous: diagram: 1-D plots
