solubility: Equilibrium Chemical Activities of Species

Description

Calculate chemical activities of aqueous species in equilibrium with a mineral or gas.

Usage

solubility(iaq, ..., in.terms.of = NULL, dissociate = FALSE, find.IS = FALSE)

Arguments

iaq

character (names) or numeric (species indices) of aqueous species

...

arguments for affinity or mosaic (i.e. plotting variables)

in.terms.of

character, express the total solubility in terms of moles of this species

dissociate

logical, does the mineral undergo a dissociation reaction?

find.IS

logical, find the equilibrium ionic strength by iteration?

Backward Compatibility

For compatibility with previous versions of the function, the iaq argument can be the output of affinity or mosaic for aqueous species. The examples for ionic strength and dissociation reactions were designed for this calling style.

In this case the (single) mineral or gas being dissolved is taken from the current basis species. Usually, the basis species should be set up so that the first basis species represents the substance being dissolved (a mineral such as CaCO or gas such as CO). This is treated as the conserved basis species, so it must be present in all of the formation reactions of the aqueous species.

The species should be defined to represent one set of aqueous species (including ions and/or neutral complexes) formed in solution, all involving the conserved basis species. For a dissociation reaction, the second basis species should be used to represent the counterion (cation or anion). Other variables (pH, ionic strength, activities of other basis species) should be defined in the call to affinity to make iaq.

Dissociation Reactions

The function perfoms some additional steps to calculate the solubility of something that dissociates (not just dissolves). For example, the dissolution of calcite (CaCO), involves the release of both calcium ions and different forms of carbonate in solution, depending on the pH. The equilibrium calculation must take account of the total activity of the shared ion (Ca), which is unknown at the start of the calculation. The solution is accomplished by recalculating the affinities, essentially working backward from the assumption that the dissociation didn't occur. The resulting activities correspond to equilibrium considering the system-wide activity of Ca.

A not recommended alternative is to set dissociate to a numeric value corresponding to the number of dissociated species (i.e. 2 for a 1:1 electrolyte). This setting indicates to calculate activities on a per-reaction basis, where each reaction has its own (independent) activity of Ca. That does not give a complete equilibrium in the system, but may be required to reproduce some published diagrams (see comment in the calcite example of demo("solubility")).

Ionic Strength

Set find.IS to TRUE to determine the final ionic strength due to dissolution of a substance in pure water. This works by calculating the ionic strength from the amounts of aqueous species formed, then re-running affinity with the calculated IS value. Note that for dissociation reactions, the ionic strength is calculated from both the ions present in the species definition and the counter ion, which should be the second basis species. The calculation is iterated until the ionic strength deviation at every point is lower than a preset tolerance (1e-4). Alternatively, speciation of counterions (e.g. ionized forms of carbonate or sulfate) can also be accomplished by using the mosaic function instead of affinity; this is used in the second example below.

Warning

This function has not been tested for systems that may form dimers or higher-order complexes (such as AuS). Except for relatively simple systems, even after careful refinement, the results from CHNOSZ, which considers chemical activities as the independent variables, will not match the results from speciation-solubility (or Gibbs energy minimization) codes, where the system is defined by its bulk composition.

Details

solubility calculates the activities of aqueous species in equilibrium with one or more minerals or gases. The minerals or gases should be loaded as the formed species, and the aqueous species (including ions and/or neutral complexes) that can be produced by dissolution should be listed in the iaq argument. The definitions of plotting variables should be provided in ..., which are passed as arguments to affinity, or to mosaic if the first one is named bases.

It must be possible to obtain a valid set of basis species by substituting each of the minerals or gases in the first position of the current basis defintion, and all of the aqueous species must include that basis species in their formation reactions. (This essentially means that all minerals, gases and aqueous species must share a common element, which is what the reactions are balanced on.)

For a single mineral or gas, the output of solubility can be used by diagram with type = NULL (the default) to plot the activities of the aqueous species or with type = "loga.balance" to plot the sum of activities of aqueous species, which corresponds to the solubility of the mineral or gas. This value corresponds to the total extent of dissolution of the mineral or gas; in.terms.of can be used to express this value in terms of another species or element. For example, for dissolution of gaseous S, in.terms.of = "S" gives the total amount of S in solution, which is twice the amount of S dissolved. Likewise, the solubility of corundum (AlO) can be expressed in terms of the moles of Al in solution (see the vignette ).

For multiple minerals, the function calculates the solubilities for each of the minerals separately; these are stored in the loga.equil element of the output. The overall minimum solubility among all the minerals at each point is stored in loga.balance. This corresponds to the total activity of dissolved species in equilibrium with the most stable mineral. In contrast to the situation for a single mineral or gas, diagram by default plots loga.balance; type = "loga.equil" can be used to plot the solubilities for the individual minerals or gases.

Examples

Run this code

# NOT RUN {
## EXAMPLE 1

# Calculate solubility of a single substance:
# Gaseous SO2 with a given fugacity
# Define basis species (any S-bearing basis species should be first)
basis(c("sulfur", "oxygen", "H2O", "H+"))
basis("pH", 6)
# Load the substances (minerals or gases) to be dissolved
species("sulfur dioxide", -20)
# List the formed aqueous species
# We can use retrieve() to identify all the possible aqueous species
iaq <- retrieve("S", c("O", "H"), "aq")
# Place arguments for affinity() after the first argument of solubility()
s1 <- solubility(iaq, O2 = c(-56, -46), T = 125, in.terms.of = "S")

# Calculate overall solubility for multiple substances:
# Gaseous S2 and SO2 with a given fugacity
basis(c("sulfur", "oxygen", "H2O", "H+"))
basis("pH", 6)
species(c("S2", "sulfur dioxide"), -20)
s2 <- solubility(iaq, O2 = c(-56, -46), T = 125, in.terms.of = "S")

# Make expressions for legends
S_ <- expr.species("SO2", "gas", -20, TRUE)
pH_ <- quote(pH == 6)
T_ <- lT(125)
lexpr <- lex(S_, pH_, T_)

# Make diagrams from the results of solubility calculations
layout(matrix(c(1, 3, 2, 3), nrow = 2))
# Logarithm of activity of aqueous species in equilibrium with SO2(gas)
diagram(s1, ylim = c(-15, 0))
diagram(s1, type = "loga.balance", col = 3, lwd = 3, add = TRUE)
legend("topright", legend = lexpr, bty = "n")
# Logarithm of concentration (parts per million) of aqueous species
sppm <- convert(s1, "logppm")
diagram(sppm, ylim = c(-10, 5))
diagram(sppm, type = "loga.balance", col = 3, lwd = 3, add = TRUE)
legend("topright", legend = lexpr, bty = "n")
par(xpd = NA)
text(-58, 6.5, paste("Solubility of gaseous SO2 (green line) is",
  "sum of concentrations of aqueous species"), cex = 1.2, font = 2)
par(xpd = FALSE)

# Show overall (minimum) solubility of multiple gases
diagram(s2, col = 4, lwd = 3)
# Show solubilities of individual gases
names <- info(species()$ispecies)$formula
diagram(s2, type = "loga.equil", names = names, add = TRUE)
title("Minimum solubility (blue line) corresponds to the most stable gas")
layout(matrix(1))

## EXAMPLE 2

## Two ways to calculate pH-dependent solubility of calcite
## with ionic strength determination
## Method 1: CO2 and carbonate species as formed species
basis(c("CO2", "Ca+2", "H2O", "O2", "H+"))
species("calcite")
iaq <- info(c("CO2", "HCO3-", "CO3-2"))
# Ionic strength calculations don't converge below around pH = 3
sa0 <- solubility(iaq, pH = c(4, 14), dissociate = TRUE)
saI <- solubility(iaq, pH = c(4, 14), dissociate = TRUE, find.IS = TRUE)
## Method 2: CO2 and carbonate species as basis species
basis(c("Ca+2", "CO2", "H2O", "O2", "H+"))
species("calcite")
iaq <- info("Ca+2")
bases <- c("CO2", "HCO3-", "CO3-2")
sm0 <- solubility(iaq, bases = bases, pH = c(4, 14), dissociate = TRUE)
smI <- solubility(iaq, bases = bases, pH = c(4, 14), dissociate = TRUE, find.IS = TRUE)
## Plot the results
plot(0, 0, xlab="pH", ylab="solubility, log mol", xlim = c(4, 14), ylim = c(-5, 2))
# Method 1 with/without ionic strength
lines(saI$vals[[1]], saI$loga.balance, lwd = 5, col = "lightblue")
lines(sa0$vals[[1]], sa0$loga.balance, lwd = 5, col = "pink")
# Method 2 with/without ionic strength
lines(smI$vals[[1]], smI$loga.balance, lty = 2)
lines(sm0$vals[[1]], sm0$loga.balance, lty = 2)
legend("topright", c("I = 0", "I = calculated", "mosaic method"),
       col = c("pink", "lightblue", "black"), lwd = c(5, 5, 1), lty = c(1, 1, 2))
title(main = "Solubility of calcite: Ionic strength and mosaic method")
# The two methods give nearly equivalent results
stopifnot(all.equal(sa0$loga.balance, sm0$loga.balance))
stopifnot(all.equal(saI$loga.balance, smI$loga.balance, tolerance = 0.003))
## NOTE: the second method (using mosaic) is slower, but is
## more flexible; e.g. complexes with Ca+2 could be included
# }

Run the code above in your browser using DataLab