Aqueous REE Speciation Calculations and Plotting in CHNOSZ

[lycunddemon]

Dear Prof. Dick,

I am a postdoctoral from Guangzhou Institute of Geochemistry, Chinese Academy of Sciences.
I found CHNOSZ when reading some journals, it is very powerful.
I am trying to use CHNOSZ to draw a figure of concentration comparison of REE complexes containing different ligands (like the picture below). I tried for a long time but failed, and I didn't find any similar demos.
Could you help me to draw such a picture?

taking from: Haas J R, Shock E L, Sassani D C. Rare earth elements in hydrothermal systems: estimates of standard partial molal thermodynamic properties of aqueous complexes of the rare earth elements at high pressures and temperatures[J]. Geochimica et cosmochimica acta, 1995, 59(21): 4329-4350.

Best wishes,
Yancheng Lou

[jedick]

Hi Yancheng, thanks for your question. I'll mention that CHNOSZ can't exactly recreate this diagram because it doesn't have a speciation model with mass-balance capabilities. What CHNOSZ can do is calculate relative abundances of species that all have a single element in common (lanthanum in this case), where those species may form by complexation with ligands whose distribution depends on variables like pH and temperature and with activity coefficients calculated for given ionic strength.

Here are a few lines of code to get started:

# Start with default settings
reset()

# Setup plot file
png("La_diagram.png", width = 8, height = 6, units = "in", res = 300)
par(mfrow = c(2, 2))

# Assign conditions
IS <- 1  # ionic strength (mol/kg)
T <- 300
P <- "Psat"
pH <- c(2, 7)

# Define basis species and set concentrations
basis(c("La+3", "Cl-", "F-", "H2O", "H+", "O2"))
basis("Cl-", 0)  # 1 M Cl-
basis("F-", -3)  # 0.001 M F-

The ionic strength and concentration values for Cl^- and F^- come from the legend of the figure you referenced by Haas et al. I write "concentration" instead of "activity" because the former is implied by setting non-zero ionic strength. To make a first attempt at a diagram, I tried this:

# Get aqueous species with La
iaq <- retrieve("La", ligands = c("H", "O", "F", "Cl"), state = "aq")
species(iaq)
a <- affinity(pH = pH, T = T, P = P, IS = IS)
diagram(a, alpha = TRUE, ylim = c(0, 1))
label.figure("A")

The result is plot A below. The upper two curves at high pH are similar to the reference diagram, but the predominance of LaF⁺² at lower pH doesn't look right. Note that anhydrous species' formulas are used in the database, so LaO₂H and LaO₂^- are equivalent to La(OH)₃ and La(OH)₄^-.

What seems to be happening is that HF(aq) forms at low pH and reduces the F^- available for complexation with La⁺³. The pKa for HF is about 5.6 at these conditions, calculated like this:

- subcrt(c("HF", "H+", "F-"), c(-1, 1, 1), T = T, P = P, IS = IS)$out$logK

Unlike HF(aq), the formation of HCl(aq) doesn't affect the availability of Cl^- very much over this pH range, because HCl is a strong acid. We could consider other aqueous species in the database that have Cl together with H and/or O ...

iCl <- retrieve("Cl", ligands = c("H", "O"), state = "aq")
dput(names(iCl))
# c("Cl-", "ClO3-", "ClO4-", "ClO-", "HClO", "ClO2-", "HClO2", "ClO2", "Cl2", "HCl")

... but unless conditions are very oxidizing, most of these aren't stable species. So let's just stick with the neutral acid species and their conjugate bases (which we started with).

The next calculation with the mosaic() function distributes the total concentration of F (defined after the basis() call above) between F^- and HF, and likewise for Cl. To be more exact, the total concentrations of F and Cl should also count the formation of complexes with La, but the relatively low concentration of La (1 ppb as used by Haas et al.) means that the values defined here aren't too much in error.

# Consider acid-base reactions
bases <- list(c("F-", "HF"), c("Cl-", "HCl"))
m <- mosaic(bases, pH = pH, T = T, P = P, IS = 1)
diagram(m$A.species, alpha = TRUE, ylim = c(0, 1))
label.figure("B")

This plot (B) is starting to look more like the reference diagram. A notable difference at low pH is the low La⁺³ concentration calculated with CHNOSZ. One guess is that there is some species (possibly a mineral) that takes up Cl and because of mass balance lowers the concentration of Cl^-. Since we haven't accounted for this, the relative abundance of La-Cl complexes is overestimated.

Another consideration is activity coefficients. Changing from the default (B-dot) to bgamma (b_γ) model increases the concentration of La⁺³ by a noticeable amount, but still doesn't get to the level seen in the reference diagram.

# Change activity coefficient model
thermo("opt$nonideal" = "bgamma")
m <- mosaic(bases, pH = pH, T = T, P = P, IS = 1)
diagram(m$A.species, alpha = TRUE, ylim = c(0, 1))
label.figure("C")
dev.off()

That's probably about as close as CHNOSZ can get to the diagram you sent. Besides figuring out what's keeping La⁺³ relatively stable at low pH, another thing to consider is how the total concentration of dissolved F changes with pH. In this example it has been set to a constant, but it would be affected by the solubility of fluorite, which possibly could be calculated with CHNOSZ or another way.

I hope that helps!

[lycunddemon]

Hi Dick,

Thank you very much, I learned a lot from your answer!

Best wishes,
Yancheng