Diagrams with multiple metals

options(width = 80)
## use pngquant to optimize PNG images
knit_hooks$set(pngquant = hook_pngquant)
pngquant <- "--speed=1 --quality=0-25"
if (!nzchar(Sys.which("pngquant"))) pngquant <- NULL
## logK with a thin space 20200627
logK <- "log&thinsp;<i>K</i>"

## Resolution settings
# Change this to TRUE to make high-resolution plots
# (default is FALSE to save time in CRAN checks)
hires <- FALSE
res1.lo <- 150
res1.hi <- 256
res1 <- ifelse(hires, res1.hi, res1.lo)
res2.lo <- 200
res2.hi <- 400
res2 <- ifelse(hires, res2.hi, res2.lo)

This vignette was compiled on r Sys.Date() with CHNOSZ version r sessionInfo()$otherPkgs$CHNOSZ$Version.

This vignette is associated with a paper that has been published in Applied Computing and Geosciences (Dick, 2021). Please consider citing that paper if you use the functions or diagrams described here.

The plots in this vignette were made using the following resolution settings, which can be changed if desired (low resolutions are used to save time in CRAN checks):

cat(paste0(ifelse(hires, "# ", ""), "res1 <- ", res1.lo))
cat(paste0(ifelse(hires, "", "# "), "res1 <- ", res1.hi))
cat(paste0(ifelse(hires, "# ", ""), "res2 <- ", res2.lo))
cat(paste0(ifelse(hires, "", "# "), "res2 <- ", res2.hi))

Basic diagrams in CHNOSZ are made for reactions that are balanced on an element (see Equilibrium in CHNOSZ) and therefore represent minerals or aqueous species that all have one element, often a metal, in common. The package documentation has many examples of diagrams for a single metal appearing in different minerals or complexed with different ligands, but a common request is to make diagrams for multiple metals. This vignette describes some methods for constructing diagrams for multi-metal minerals and other multi-element systems. The methods are mashing, mixing, mosaic stacking, and secondary balancing.


Mashing or simple overlay refers to independent calculations for two different systems that are displayed on the same diagram.

This example starts with a logf~O2~--pH base diagram for the C-O-H system then overlays a diagram for S-O-H. The second call to affinity() uses the argument recall feature, where the arguments after the first are taken from the previous command. This allows calculations to be run at the same conditions for a different system. This feature is also used in other examples in this vignette.

par(mfrow = c(1, 2))
species(c("CH4", "CO2", "HCO3-", "CO3-2"))
aC <- affinity(pH = c(0, 14), O2 = c(-75, -60))
dC <- diagram(aC, dx = c(0, 1, 0, 0), dy = c(0, 1, 0, 0))
species(c("H2S", "HS-", "HSO4-", "SO4-2"))
aS <- affinity(aC)  # argument recall
dS <- diagram(aS, add = TRUE, col = 4, col.names = 4, dx = c(0, -0.5, 0, 0))
aCS <- mash(dC, dS)
srt <- c(0, 0, 90, 0, 0, 0, 90, 0, 0, 0)
cex.names <- c(1, 1, 0.8, 1, 1, 1, 1, 1, 1, 1)
dy <- c(0, 0, 0, -0.2, 0, 0, 0, 0, 0, 0)
diagram(aCS, srt = srt, cex.names = cex.names, dy = dy)
legend("topright", legend = lTP(25, 1), bty = "n")

The second diagram is just like the first, except the function mash() is used to label the fields with names of species from both systems, and a legend is added to indicate the temperature and pressure.

Note that these are predominance diagrams, so they show only the species with highest activity; there is in fact a distribution of activities of aqueous species that is not visible here.

Tip: the names of the fields in the second diagram come from aCS$species$name, which are expressions made by combining aC$names and aS$names. If you prefer plain text names without formatting, add format.names = FALSE to all of the diagram() calls.

Mixing 1

As shown above, mashing two diagrams is essentially a simple combination of the two systems. Although it is easy to make such a diagram, there is no interaction between the systems. If there is a possibility of forming bimetallic species, then additional considerations are needed to account for the stoichiometry of the mixture. The stoichiometry can be given as a fixed composition of both metals; then, all combinations of (mono- and/or bimetallic) species that satisfy this compositional constraint are used as the candidate "species" in the system. This is the same type of calculation as that described for binary Pourbaix diagrams in the Materials Project.

This example makes a Pourbaix diagram for the Fe-V-O-H system that is similar to Figure 1 of @SZS_17. Before getting started, it may be helpful to clarify some terminology. In the materials science community, materials are characterized by several energies (among others): 1) the formation energy from the elements in their reference state, 2) the energy above the convex hull, which is zero for stable materials, and greater than zero for metastable materials, and 3) the Pourbaix energy difference (ΔG~pbx~), which refers to the energy of a given material with respect to the stable solids and aqueous species as a function of pH and Eh. The parallel terminology used in CHNOSZ is that aqueous species or minerals have a 1) standard Gibbs energy of formation from the elements, ΔG° = f(T, P), which is available from the OBIGT database, 2) standard Gibbs energy of reaction from the stable species, and 3) affinity of formation from the basis species, A = -ΔG = f(T, P, and activities of all species). As used in CHNOSZ, "formation reaction" refers to formation from the basis species, not from the elements. The basis species are not in general the stable species, so we begin by identifying the stable species in the system; the difference between their affinities and the affinity of any other species corresponds to -ΔG~pbx~.

First we need to assemble the standard Gibbs energies of the solids and aqueous species. For solids, values of formation energy from the elements (in eV/atom) computed using density functional theory (DFT) were retrieved from the Materials API [@OCJ_15] and are converted to units of J/mol. The Materials Project (MP) website also provides these values, but with fewer decimal places, which would lead to a small rounding error in the comparison of energy above the hull at the end of this example. For aqueous species, values of standard Gibbs energy of formation from the elements at 25 °C (in J/mol) are taken mostly from @WEP_82 augmented with data for FeO~2~^-^ from @SSWS97 and FeO~4~^-2^ from @Mis73. Adapting the method described by @PWLC12, a correction for each metal is calculated from the difference between the DFT-based formation energy and the standard Gibbs energy of a representative material; here we use values for Fe~3~O~4~ (magnetite) and V~3~O~5~ from @WEP_82. This correction is then applied to all of the aqueous species that have that metal. Finally, mod.OBIGT() is used to add the obtained energies to the OBIGT database in CHNOSZ.

This code produces species indices in the OBIGT database for Fe- and V-bearing aqueous species (iFe.aq, iV.aq), solids (iFe.cr, iV.cr), and bimetallic solids (iFeV.cr), which are used in the following diagrams.

Now we set up the plot area and assign activities of aqueous species to 10^-5^, which is the default value for diagrams on the MP website (from the page for a material: "Generate Phase Diagram" -- "Aqueous Stability (Pourbaix)"). The following commands compute Eh-pH diagrams for the single-metal systems Fe-O-H and V-O-H. The pH and Eh ranges are made relatively small in order to show just a part of the diagram. The diagrams are not plotted, but the output of diagram() is saved in dFe and dV for later use.

par(mfrow = c(1, 3))
loga.Fe <- -5
loga.V <- -5
# Fe-O-H diagram
basis(c("VO+2", "Fe+2", "H2O", "e-", "H+"))
species(c(iFe.aq, iFe.cr))
species(1:length(iFe.aq), loga.Fe)
aFe <- affinity(pH = c(4, 10, res1), Eh = c(-1.5, 0, res1))
dFe <- diagram(aFe, plot.it = FALSE)
# V-O-H diagram
species(c(iV.aq, iV.cr))
species(1:length(iV.aq), loga.V)
aV <- affinity(aFe)  # argument recall
dV <- diagram(aV, plot.it = FALSE)

# Calculate affinities for bimetallic species
aFeV <- affinity(aFe)  # argument recall
dFeV <- diagram(aFeV, plot.it = FALSE, bold = TRUE)
# 1:1 mixture (Fe:V)
a11 <- mix(dFe, dV, dFeV, c(1, 1))
# Adjust labels 20210219
iV2O3 <- info("V2O3")
iFeO <- info("FeO", "cr")
iFe3V <- info("Fe3V")
srt <- rep(0, nrow(a11$species))
srt[a11$species$ispecies == paste(iFeO, iV2O3, sep = ",")] <- 90
srt[a11$species$ispecies == paste(iV2O3, iFe3V, sep = ",")] <- -13
diagram(a11, min.area = 0.01, srt = srt)
title("Fe:V = 1:1")
label.figure(lTP(25, 1), xfrac = 0.12)
# 1:3 mixture
a13 <- mix(dFe, dV, dFeV, c(1, 3))
srt <- rep(0, nrow(a13$species))
srt[a13$species$ispecies == paste(iFeO, iV2O3, sep = ",")] <- 90
srt[a13$species$ispecies == paste(iV2O3, iFe3V, sep = ",")] <- -13
diagram(a13, min.area = 0.01, srt = srt)
title("Fe:V = 1:3")
# 1:5 mixture
a15 <- mix(dFe, dV, dFeV, c(1, 5))
iFeV3 <- info("FeV3")
srt <- rep(0, nrow(a15$species))
srt[a15$species$ispecies == paste(iFeO, iV2O3, sep = ",")] <- 90
srt[a15$species$ispecies == paste(iV2O3, iFe3V, sep = ",")] <- -13
srt[a15$species$ispecies == paste(iV2O3, iFeV3, sep = ",")] <- -13
diagram(a15, min.area = 0.01, srt = srt)
title("Fe:V = 1:5")

Then we calculate the affinities for the bimetallic species and save the output of diagram() in dFeV, again without making a plot, but formatting the names in bold. Note that diagram() uses different colors for regions with two solids, one solid, and no solids, including some transparency to show the underlying water stability region that is plotted first.

Now we have all the ingredients needed to combine the Fe-bearing, V-bearing, and bimetallic species to generate a given composition. The mix() function is used to calculate the affinities of formation from basis species for all combinations of aqueous species and minerals that satisfy each of three different compositions. Finally, the diagram()s are plotted; the min.area argument is used to remove labels for very small fields. Regarding the legend, it should be noted that although the DFT calculations for solids are made for zero temperature and zero pressure [@SZS_17], the standard Gibbs energies of aqueous species [e.g. @WEP_82] are modified by a correction term so that they can be combined with DFT energies to reproduce the experimental energy for dissolution of a representative material for each metal at 25 °C and 1 bar [@PWLC12].

In these diagrams, changing the Fe:V ratio affects the fully reduced metallic species. In the 1:1 mixture, the FeV~3~ + Fe~3~V assemblage is predicted to be stable instead of FeV. This result is unlike Figure 1 of @SZS_17 but is consistent with the MP page for FeV where it is shown to decompose to this assemblage. On the other hand, FeV~3~ is stable in the 1:3 mixture. For an even higher proportion of V, the V + FeV~3~ assemblage is stable, which can be seen for instance in the Pourbaix diagram linked from the MP page for FeV~5~O~12~.

Let's make another diagram for the 1:1 Fe:V composition over a broader range of Eh and pH. The diagram shows a stable assemblage of Fe~2~O~3~ with an oxidized bimetallic material, Fe~2~V~4~O~13~.

layout(t(matrix(1:3)), widths = c(1, 1, 0.2))
par(cex = 1)
# Fe-bearing species
basis(c("VO+2", "Fe+2", "H2O", "e-", "H+"))
species(c(iFe.aq, iFe.cr))$name
species(1:length(iFe.aq), loga.Fe)
aFe <- affinity(pH = c(0, 14, res2), Eh = c(-1.5, 2, res2))
dFe <- diagram(aFe, plot.it = FALSE)
# V-bearing species
species(c(iV.aq, iV.cr))$name
species(1:length(iV.aq), loga.V)
aV <- affinity(aFe)  # argument recall
dV <- diagram(aV, plot.it = FALSE)
# Bimetallic species
aFeV <- affinity(aFe)  # argument recall
dFeV <- diagram(aFeV, plot.it = FALSE, bold = TRUE)
# 1:1 mixture (Fe:V)
a11 <- mix(dFe, dV, dFeV, c(1, 1))
# Adjust labels 20210219
iV2O3 <- info("V2O3")
iFe3V <- info("Fe3V")
iVO4m3 <- info("VO4-3")
iFe2O3 <- info("Fe2O3")
srt <- rep(0, nrow(a11$species))
srt[a11$species$ispecies == paste(iV2O3, iFe3V, sep = ",")] <- -13
srt[a11$species$ispecies == paste(iFe2O3, iVO4m3, sep = ",")] <- 90
d11 <- diagram(a11, min.area = 0.01, srt = srt)
water.lines(d11, col = "orangered")

# Calculate affinity of FeVO4
aFeVO4 <- affinity(aFe)  # argument recall
# Calculate difference from stable species
aFeVO4_vs_stable <- aFeVO4$values[[1]] - d11$predominant.values
# Overlay lines from diagram on color map
diagram(a11, fill = NA, names = FALSE, limit.water = FALSE)
opar <- par(usr = c(0, 1, 0, 1))
col <- rev(topo.colors(128)) # No hcl.colors() in R < 3.6.0
if(getRversion() >= "3.6.0") col <- rev(hcl.colors(128, palette = "YlGnBu", alpha = 0.8))
image(aFeVO4_vs_stable, col = col, add = TRUE)
diagram(a11, fill = NA, add = TRUE, names = FALSE)
water.lines(d11, col = "orangered")

imax <- arrayInd(which.max(aFeVO4_vs_stable), dim(aFeVO4_vs_stable))
pH <- d11$vals$pH[imax[1]]
Eh <- d11$vals$Eh[imax[2]]
points(pH, Eh, pch = 10, cex = 2, lwd = 2, col = "gold")
stable <- d11$names[d11$predominant[imax]]
text(pH, Eh, stable, adj = c(0.5, -1), cex = 1.2, col = "gold")

# Make color scale 20210228
par(mar = c(3, 0, 2.5, 2.7))
levels <- 1:length(col)
plot.window(xlim = c(0, 1), ylim = range(levels), xaxs = "i", yaxs = "i")
rect(0, levels[-length(levels)], 1, levels[-1L], col = rev(col), border = NA)
# To get the limits, convert range of affinities to eV/atom
arange <- rev(range(aFeVO4_vs_stable))
# This gets us to J/mol
Jrange <- convert(arange, "G")
# And to eV/atom
eVrange <- Jrange / 1.602176634e-19 / 6.02214076e23 / 6
ylim <- formatC(eVrange, digits = 3, format = "f")
axis(4, at = range(levels), labels = ylim)
mtext(quote(Delta*italic(G)[pbx]*", eV/atom"), side = 4, las = 0, line = 1)

We then compute the affinity for formation of a metastable material, in this case triclinic FeVO~4~, from the same basis species used to make the previous diagrams. Given the diagram for the stable Fe-, V- and bimetallic materials mixed with the same stoichiometry as FeVO~4~ (1:1 Fe:V), the difference between their affinities of formation and that of FeVO~4~ corresponds to the Pourbaix energy difference (-ΔG~pbx~). This is plotted as a color map in the second diagram. (See the source of this vignette for the code used to make the scale bar.)

Now we locate the pH and Eh that maximize the affinity (that is, minimize ΔG~pbx~) of FeVO~4~ compared to the stable species. In agreement with @SZS_17, this is in the stability field of Fe~2~O~3~ + Fe~2~V~4~O~13~.

plot(1:10) # so we can run "points" in this chunk
imax <- arrayInd(which.max(aFeVO4_vs_stable), dim(aFeVO4_vs_stable))
pH <- d11$vals$pH[imax[1]]
Eh <- d11$vals$Eh[imax[2]]
points(pH, Eh, pch = 10, cex = 2, lwd = 2, col = "gold")
stable <- d11$names[d11$predominant[imax]]
text(pH, Eh, stable, adj = c(0.3, 2), cex = 1.2, col = "gold")
(Apbx <- range(aFeVO4_vs_stable[d11$predominant == d11$predominant[imax]]))

Although one point is drawn on the diagram, FeVO~4~ has the same Pourbaix energy difference with respect to the entire Fe~2~O~3~ + Fe~2~V~4~O~13~ field, as shown by the range() command (the values are dimensionless values of affinity, A/(RT) = -ΔG~pbx~/(RT)). This can occur only if the decomposition reaction has no free O~2~ or H~2~, and means that in this case ΔG~pbx~ in the Fe~2~O~3~ + Fe~2~V~4~O~13~ field is equal to the energy above the hull.

To calculate the energy above the hull "by hand", let's set up the basis species to be the stable decomposition products we just found. O~2~ is also needed to make a square stoichiometric matrix (i.e. same number of elements and basis species), but it does not appear in the reaction to form FeVO~4~ from the basis species. subcrt() is used to automatically balance the formation reaction for 1 mole of FeVO~4~ and calculate the standard Gibbs energy of the reaction. We first test that r logK of the reaction (calculated with convert(), which divides ΔG° by -RT) is the same as the dimensionless affinity for FeVO~4~ calculated above. Then, the value of ΔG° in J/mol is converted to eV/mol, and finally eV/atom.

b <- basis(c("Fe2O3", "Fe2V4O13", "O2"))
J_mol <- subcrt("FeVO4", 1, T = 25)$out$G
stopifnot(all.equal(rep(convert(J_mol, "logK"), 2), Apbx))
eV_mol <- J_mol / 1.602176634e-19
eV_atom <- eV_mol / 6.02214076e23 / 6
round(eV_atom, 3)
stopifnot(round(eV_atom, 3) == 0.415)

This is equal to the value for the energy above the hull / atom for triclinic FeVO~4~ on the MP website (0.415 eV, accessed on 2020-11-09 and 2021-02-19). This shows that we successfully made a round trip starting with the input formation energies (eV/atom) from the Materials API, to standard Gibbs energy (J/mol) in the OBIGT database, and back out to energy above the hull (eV/atom).

The concept of using the stable minerals and aqueous species to calculate reaction energetics is formalized in the mosaic() function, which is described next. Because this example modified the thermodynamic data for some minerals that are used below, we should restore the default OBIGT database before proceeding to the next section.


Mosaic Stacking 1

For a function that implements the workflow described below, see stack_mosaic() (added in CHNOSZ 2.0.0).

A mosaic diagram shows the effects of changing basis species on the stabilities of minerals. The Fe-S-O-H system is a common example: the speciation of aqueous sulfur species affects the stabilities of iron oxides and sulfides. Examples of mosaic diagrams with Fe or other single metals are given elsewhere.

A mosaic stack is when predominance fields for minerals calculated in one mosaic diagram are used as input to a second mosaic diagram, where the minerals are now themselves basis species. The example here shows the construction of a Cu-Fe-S-O-H diagram.

First we define the conditions and basis species. It is important to put Cu^+^ first so that it will be used as the balance for the reactions with Cu-bearing minerals (which also have Fe). Pyrite is chosen as the starting Fe-bearing basis species, which will be changed as indicated in Fe.cr.

logaH2S <- -2
T <- 200
pH <- c(0, 14, res2)
O2 <- c(-48, -33, res2)
basis(c("Cu+", "pyrite", "H2S", "oxygen", "H2O", "H+"))
basis("H2S", logaH2S)
S.aq <- c("H2S", "HS-", "HSO4-", "SO4-2")
Fe.cr <- c("pyrite", "pyrrhotite", "magnetite", "hematite")
Fe.abbrv <- c("Py", "Po", "Mag", "Hem")

Now we calculate affinities for minerals in the Fe-S-O-H system that take account of the changing aqueous sulfur species in S.aq. The result is used to make different layers of the diagram (1 and 2 are both made by the first call to diagram()):

  1. Water stability region (gray shading)
  2. Predominance fields for the aqueous S species (blue text and dashed lines)
  3. Stability areas for the Fe-bearing minerals (black text and lines)
mFe <- mosaic(S.aq, pH = pH, O2 = O2, T = T)
diagram(mFe$A.bases, lty = 2, col = 4, col.names = 4, italic = TRUE, dx = c(0, 1, 0, 0), dy = c(-1.5, 0, 1, 0))
dFe <- diagram(mFe$A.species, add = TRUE, lwd = 2, names = Fe.abbrv, dx = c(0, 0.5, 0, 0), dy = c(-1, 0, 0.5, 0))
FeCu.cr <- c("chalcopyrite", "bornite")
Cu.cr <- c("copper", "cuprite", "tenorite", "chalcocite", "covellite")
FeCu.abbrv <- c("Ccp", "Bn", "Cu", "Cpr", "Tnr", "Cct", "Cv")
species(c(FeCu.cr, Cu.cr))
mFeCu <- mosaic(list(S.aq, Fe.cr), pH = pH, O2 = O2,
              T = T, stable = list(NULL, dFe$predominant))
col <- c("#FF8C00", rep(2, 6))
lwd <- c(2, 1, 1, 1, 1, 1, 1)
dy = c(0, 0, 0, 0, 0, 1, 0)
diagram(mFeCu$A.species, add = TRUE, col = col, lwd = lwd, col.names = col, bold = TRUE, names = FeCu.abbrv, dy = dy)
TPS <- c(describe.property(c("T", "P"), c(T, "Psat")), expression(sum(S) == 0.01*m))
legend("topright", TPS, bty = "n")
title("Cu-Fe-S-O-H (minerals only)", font.main = 1)

Next we load the Cu-bearing minerals and calculate their affinities while changing both the aqueous sulfur species and the Fe-bearing minerals whose stability fields were just calculated. The latter step is the key to the mosaic stack and is activated by supplying the calculated stabilities of the Fe-bearing minerals in the stable argument. This is a list whose elements correspond to each group of changing basis species given in the first argument. The NULL means that the abundances of S-bearing aqueous species are calculated according to the default in mosaic(), which uses equilibrate() to compute the continuous transition between them ("blending"). Because the Fe-bearing minerals are the second group of changing basis species (Fe.cr), their stabilities are given in the second position of the stable list. The result is used to plot the last layer of the diagram:

  1. Stability areas for Cu-bearing minerals (red text and lines; orange for chalcopyrite)

After that we add the legend and title.

This diagram has a distinctive chalcopyrite "hook" surrounded by a thin bornite field. Only the chalcopyrite-bornite reaction in the pyrite field is shown in some published diagrams [e.g. @And75;@Gio02], but diagrams with a similar chalcopyrite wedge or hook geometry can be seen in @BBR77 and @Bri80.

Mosaic Stacking 2

The previous diagram shows the relative stabilities of minerals only. The next diagram adds aqueous species to the system. The position of the boundaries between the stability fields for minerals and aqueous species are calculated for a given activity for the latter, in this case 10^-6^.

After running the code to make this diagram, we can list the reference keys for the minerals and aqueous species.

minerals <- list(Fe.cr = Fe.cr, Cu.cr = Cu.cr, FeCu.cr = FeCu.cr)
aqueous <- list(S.aq = S.aq, Fe.aq = Fe.aq, Cu.aq = Cu.aq)
allspecies <- c(minerals, aqueous)
iall <- lapply(allspecies, info)
allkeys <- lapply(iall, function(x) thermo.refs(x)$key)

The next code chunk prepends @ to the reference keys and uses the chunk option results = 'asis' to insert the citations into the R Markdown document in chronological order.

allyears <- lapply(iall, function(x) thermo.refs(x)$year)
o <- order(unlist(allyears))
cat(paste(paste0("@", unique(unlist(allkeys)[o])), collapse = "; "))

Mixing 2

The previous diagram shows a stability boundary between chalcopyrite and bornite but does not identify the stable assemblages that contain these minerals. This is where mix() can help. Following the workflow described in Mixing 1, we first calculate individual diagrams for Fe-S-O-H and Cu-S-O-H, which are overlaid on the first plot and saved in dFe and dCu. We then calculate the affinities for the bimetallic Cu and Fe minerals and run them through diagram() without actually making a plot, but save the result in dFeCu. Then, we combine the results using mix() to define different proportions of Fe and Cu.

These diagrams show that changing the amounts of the metals affects the stability of minerals involved in reactions with chalcopyrite. At a 1:1 ratio of Fe:Cu, chalcopyrite is a stable single-mineral assemblage. At a 2:1 ratio, pyrite, pyrrhotite, or magnetite can coexist in a two-phase assemblage with chalcopyrite. At a 1:2 ratio, an assemblage consisting of the two bimetallic minerals (chalcopyrite and bornite) is stable.

Mosaic Stacking 3

The results of a mosaic stack can also be processed with mash() to label each region with the minerals from both systems. For this example, the speciation of aqueous sulfur is not considered; instead, the fugacity of S~2~ is a plotting variable. The stable Fe-bearing minerals (Fe-S-O-H) are used as the changing basis species to make the diagram for Cu-bearing minerals (Cu-Fe-S-O-H). Then, the two diagrams are mashed to show all minerals in a single diagram. Greener colors are used to indicate minerals with less S~2~ and more O~2~ in their formation reactions.

The resulting diagram is similar to Figure 2 of @Sve87; that diagram also shows calculations of the solubility of Cu and concentration of SO~4~^-2^ in model Cu ore-forming fluids. The solubility() function can be used to calculate the total concentration of Cu in different complexes in solution (listed in the iaq argument). The bases argument triggers a mosaic() calculation, so that the solubility corresponds that that of stable minerals at each point on the diagram. The pH for these calculations is set to 6, and the molality of free Cl^-^, which affects the formation of the Cu chloride complexes, is estimated based on the composition of fluids from Table 2 of @Sve87 (ca. 80000 mg Cl / kg H~2~O) and the NaCl() function in CHNOSZ. This also gives an estimated ionic strength, which is used in the following mosaic() and affinity() calls to calculate activity coefficients.

After running the code above, we can inspect the value of calc to show the estimated ionic strength and activity of Cl^-^; the latter is very close to unity.

# Ionic strength
# Logarithm of activity of Cl-
log10(calc$m_Cl * calc$gam_Cl)

The thick magenta lines indicate the 35 ppm contour for Cu and SO~4~^-2^. The first plot shows a lower Cu solubility in this region compared to Figure 2 of @Sve87. The difference is likely due to lower stabilities of Cu(I) chloride complexes in the default OBIGT database, compared to those available at the time [@Hel69]. For the second plot, the standard Gibbs energies of CuCl~2~^-^ and CuCl~3~^-2^ are adjusted so that the logK for their dissociation reactions at 125 °C matches values interpolated from Table 5 of @Hel69. Here are the logK values after the adjustment, followed a reset() call to compare the values with the default database, which is also used for later examples in this vignette. (T was set to 125 above.)

# logK values interpolated from Table 5 of Helgeson (1969)
subcrt(c("CuCl2-", "Cu+", "Cl-"), c(-1, 1, 2), T = T)$out$logK
subcrt(c("CuCl3-2", "Cu+", "Cl-"), c(-1, 1, 3), T = T)$out$logK
# Default OBIGT database
subcrt(c("CuCl2-", "Cu+", "Cl-"), c(-1, 1, 2), T = T)$out$logK
subcrt(c("CuCl3-2", "Cu+", "Cl-"), c(-1, 1, 3), T = T)$out$logK

The higher stability of these complexes from @Hel69 causes the 35 ppm contour to move closer to the position shown by @Sve87.

Interestingly, the calculation here also predict substantial increases of Cu concentration of Cu at high f~S2~ and low f~O2~, due to the formation of bisulfide complexes with Cu. The aqueous species considered in the calculation can be seen like this:


CuHS and Cu(HS)~2~^-^ can be excluded by removing S from the retrieve() call above (i.e. only c("O", "H", "Cl") as the elements in possible ligands); doing so precludes a high concentration of aqueous Cu in the highly reduced, sulfidic region.

The third plot for the concentation of SO~4~^-2^ is simply made by using affinity() to calculate the affinity of its formation reaction as a function of f~S2~ and f~O2~ at pH 6 and 125 °C, then using solubility() to calculate the solubility of S~2~(gas), expressed in terms of moles of SO~4~^-2^ in order to calculate parts per million (ppm) by weight.

Secondary Balancing

Predominance diagrams in CHNOSZ are made using the maximum affinity method, where the affinities of formation reactions of species are divided by the balancing coefficients [@Dic19]. Usually, these balancing coefficients are taken from the formation reactions themselves; for example, if they are the coefficients on the Fe-bearing basis species, then the reactions are said to be "balanced on Fe".

Some diagrams in the literature are made with secondary balancing constraints in addition to the primary ones. For example, reactions of Fe-bearing minerals are balanced on Fe, and reactions of Cu-bearing minerals are balanced on Cu; these are both primary balancing coefficients. Then, reactions between all minerals are balanced on H^+^ as the secondary balancing coefficients. The core concept is to apply the secondary balance while also maintaining the primary balance; a method to do this has been implemented in the rebalance() function.

Different parts of the script to make the diagrams are described below; press the button to show the entire script.

We first define basis species to contain both Cu- and Fe-bearing species. The \emph{x} axis is the ratio of activities of Fe^+2^ and Cu^+^; the label is made with ratlab().

We then calculate the diagrams for the primary balancing coefficients, for the groups of only Fe-, only Cu-, and only Fe+Cu-bearing minerals. It is obvious that the first two systems are balanced on Fe and Cu, respectively, but the third has a somewhat unusual balance: H^+^. See Reaction 4 of @MH85 for an example.

Now comes the secondary balancing, where all reactions, not only that between bornite and chalcopyrite, are balanced on H^+^. We first rebalance the diagrams for the Fe- or Cu-bearing minerals to make diagram D. Note that after secondary balancing with rebalance(), the argument balance = 1 should be used in diagram() to prevent further balancing. This is because rebalance() preserves the primary balancing for Fe- and Cu-bearing minerals (internally the "plotvals" components of dFe and dCu).

Then we rebalance diagrams D and C to make the final diagram in E. The fields in this diagram are labeled with mineral abbreviations from the OBIGT database.

The final diagram is like one shown in Figure 5 of @Bri80 and Figure 5 of @MH85.

Challenge: Although the diagram here is drawn only for H~2~S in the basis species, take it a step further and make a mosaic diagram to account for the stability of HSO~4~^-^ at high oxygen fugacity.

Other Possibilities

Conceptually, the methods described above treat different metal-bearing elements as parts of distinct chemical systems that are then joined together. Other methods may be more suitable for considering multiple metals (or other elements) in one system.

Balancing on a Non-Metal

As shown in the secondary balancing example, there is no requirement that the balancing coefficients come from a metal-bearing species. It is possible to make diagrams for minerals with different metallic elements simply by using a non-metallic element as the primary balance. Here is an example for the Cu-Fe-S-O-H system. The reactions are balanced on O~2~, which means that no O~2~ appears in the reaction between any two minerals, but Fe^+2^ and/or Cu^+^ can be present, depending on the chemical composition. Saturation limits are shown for species that have no O~2~ in their formation reactions.

basis(c("Fe+2", "Cu+", "hydrogen sulfide", "oxygen", "H2O", "H+"))
basis("H2S", 2)
species(c("pyrite", "magnetite", "hematite", "covellite", "tenorite",
          "chalcopyrite", "bornite"))
a <- affinity("Cu+" = c(-8, 2, 500), "Fe+2" = c(-4, 12, 500), T = 400, P = 2000)
names <- info(a$species$ispecies)$abbrv
d <- diagram(a, xlab = ratlab("Cu+"), ylab = ratlab("Fe+2"), balance = "O2", names = names)
title(bquote("Cu-Fe-S-O-H; 1° balance:" ~ .(expr.species(d$balance))))
# Add saturation lines
species(c("pyrrhotite", "ferrous-oxide", "chalcocite", "cuprite"))
asat <- affinity(a)  # argument recall
names <- asat$species$name
names[2] <- "ferrous oxide"
diagram(asat, type = "saturation", add = TRUE, lty = 2, col = 4, names = names)
legend("topleft", legend = lTP(400, 2000), bty = "n")

This example was prompted by Figure 3 of @MH85; earlier versions of the diagram are in @HBL69 and @Hel70a.

In some ways this is like the inverse of the mosaic stacking example. There, reactions were balanced on Fe or Cu, and f~O2~ and pH were used as plotting variables. Here, the reactions are balanced on O~2~ and implicitly on H^+^ through the activity ratios with a~Fe^+2^~ and a~Cu^+^~, which are the plotting variables.

More common diagrams of this type are balanced on Si or Al. See demo(saturation) for an example in the H~2~O-CO~2~-CaO-MgO-SiO~2~ system.

Mosaic Combo

Instead of adding minerals with different metals by stacking mosaic diagrams, it may be possible to include two different metals in the basis species and formed species. The mosaic() and equilibrate() functions can be combined to balance on two different elements. The example here compares two methods applied to N-, C-, and N+C-bearing species because bimetallic aqueous species are not currently available in the OBIGT database. The total activities used here are modified from the example for sedimentary basin brines described by @Sho93, which is also the source of the thermodynamic parameters for acetamide.

  1. The "effective equilibrium constant" (K^eff^) method [@RBG_21] is used to calculate the activity of acetamide for a total activities of neutral and ionized species, i.e. ∑(ammonia and ammonium) and ∑(acetic acid and acetate).

  2. Using the mosaic combo method, the mosaic() command calculates equilibrium activities of NH~3~ and NH~4~^+^ for a given total activity of N in the basis species, and calculates the corresponding affinities of the formed species. Then, the equilibrate() command calculates equilibrium activities of the formed species for given total activity of C, and combines them with the activities of the changing basis species (NH~3~ and NH~4~^+^).

The mosaic combo method (solid black line) produces results equivalent to those of the K^eff^ method (dashed blue line).

The diagram shows the ionization of acetic acid and NH~3~ at different pHs. The predicted appearance of acetamide (CH~3~CONH~2~) is a consequence of the interaction between the N-bearing and C-bearing species, and is analogous to the formation of a bimetallic complex.

Thanks to Kirt Robinson for the feature request and test case used in this example.

Document History


Try the CHNOSZ package in your browser

Any scripts or data that you put into this service are public.

CHNOSZ documentation built on March 31, 2023, 7:54 p.m.