## 3.8.2. Metal speciation models

Authors: Wilko Verweij

Reviewers: John Parsons, Stephen Lofts

Learning objectives:

You should be able to

- Understand the basics of speciation modeling
- Understand the factors determining speciation and how to calculate them
- Understand in which types of situations speciation modeling can be helpful

Keywords: speciation modeling, solubility, organic complexation

Introduction

Speciation models allow users to calculate the speciation of a solution rather than to measure it in a chemical way or to assess it indirectly using bioassays (see section 3.5). As a rule, speciation models take total concentrations as input and calculate species concentrations.

Speciation models use thermodynamic data about chemical equilibria to calculate the speciation. This data, expressed in free energy or as equilibrium constants, can be found in the literature. The term ‘constant’ is slightly misleading as equilibrium constants depend on the temperature and ionic strength of the solution. The ionic strength is calculated from the concentrations (C) and charges (Z) of ions in solution using the equation:

For many equilibria, no information is available to correct for temperature. To correct for ionic strength, many semi-empirical methods are available, none of which is perfect.

How these models work

For each equilibrium reaction, an equilibrium constant can be defined. For example, for the reaction

Cu2+ + 4 Cl- ⇌ CuCl42-

the equilibrium constant can be defined as

Consequently, when the concentrations of free Cu2+ and free Cl- are known, the concentration of CuCl42- can be easily calculated as:

[CuCl42-] = β * [Cu2+] * [Cl-]4

In fact, the concentrations of free Cu2+ and free Cl- are often NOT known, but what is known are the total concentrations of Cu and Cl in the system. In order to find the speciation, we need to set up a set of mass balance equations needs to be set up, for example:

[total Cu] = [free Cu2+] + [CuOH+] + [Cu(OH)2] + [Cu(OH)3-] (..) + [CuCl+] + [CuCl2] (..) etc.

[total Cl] = (..)

Each concentration of a complex is a function of the free concentrations of the ions that make it up. So we can say that if we know the concentrations of all the free ions, we can calculate the concentrations of all the complexes, and then we can calculate the total concentrations. A solution to the problem cannot be found by rearranging the mass balance equations, because they are non-linear. What a speciation model does is to repeatedly estimate the free ion concentrations, on each loop adjusting them so that the calculated total concentrations more closely match the known totals. When the calculated and known total concentrations all agree to within a defined precision, the speciation has been calculated. The critical part of the calculation is adjusting the free ion concentrations in a sensible and efficient way to find the solution as quickly as possible. Several more or less sophisticated methods are available to solve this, but usually a Newton-Raphson method is applied.

Influence of temperature and ionic strength

In fact the explanation above is too simple. Equilbrium constants are valid under specific conditions for temperature and ionic strength (for example the standard conditions of 25oC and [endif]--> and need to be converted to the temperature and ionic strength of the system for which speciated is being calculated. It is possible to adapt the equilibrium constants for non-standard temperatures, but this requires knowledge of heat capacity (ΔH) data of each equilibrium. That knowledge is often not available. Constants can be converted from 25°C to other temperatures using the Van ‘t Hoff-equation:

where K1 and K2 are the constants, T1 and T2 the temperatures, Δ H is the enthalpy of a reaction and R is the gas constant.

Equilbrium constants are also valid for one specific value of ionic strength. For conversion from one value of ionic strength to another, many different approaches may be used. This conversion is quite important, because already at relatively low ionic strengths, deviations from ideality become significant, and the activity of a species starts to deviate from its concentration. Hence, the intrinsic, or thermodynamic, equilibrium constants (i.e. constants at a hypothetical ionic strength of zero) are no longer valid and the activity a of ions at non-zero ionic strength needs to be calculated from the concentration and the activity coefficient:

a = γ * c

where γ is the activity coefficient (dimensionless; sometimes also called f) and c is the concentration; a and c are in mol/liter.

The first solution to calculate activity coefficients for non-zero ionic strength was proposed by Debye and Hückel in 1923. The Debye-Hückel theory assumes ions are point charges so it does not take into account the volume that these ions occupy nor the volume of the shell of ligands and/or water molecules around them. The Debye-Hückel gives good approximations, up to circa 0.01 M for a 1:1-electrolyte, but only up to circa 0.001 M for a 2:2-electrolyte. When the ionic strength exceeds these values, the activity coefficients that the Debye-Hückel approximation predicts deviate significantly from experimental values. Many environmental applications require conversions for higher ionic strengths making the Debye-Hückel-equation insufficient. To overcome this problem, many researchers have suggested other methods, like the extended Debye-Hückel-equation, the Güntelberg-equation and the Davies-equation, but also the Bromley-equation, the Pitzer-equation and the Specific Ion Interaction Theory (SIT).

Many programs use the Davies-equation, which calculates activity coefficients γ as follows:

where z is the charge of the species and I the ionic strength. Sometimes 0.2 instead of 0.3 is used. Basically all these approaches take the Debye-Hückel-equation as a starting point, and add one or more terms to correct for deviations at higher ionic strengths. Although many of these methods are able to predict the activity of ions fairly well, they are in fact mainly empirical extensions without a solid theoretical basis.

Solubility

Most salts have a limited solubility; in several cases the solubility is also important under conditions that occur in the environment. For instance, for CaCO3 the solubility product is 10-8.48, which means that when [Ca2+] * [CO32-] > 10-8.48, CaCO3 will precipitate, until [Ca2+] * [CO32-] = 10-8.48. But it also works the other way around: if solid CaCO3 is present in a solution where [Ca2+] * [CO32-] < 10-8.48 (note the ‘<’-sign), solid CaCO3 will dissolve, until [Ca2+] * [CO32-] = 10-8.48. Note that the Ca and CO3 in the formula here refer to free ions. For example, a 10-13 M solution of Ag2S will lead to precipitation of Ag2S. The free concentrations of Ag and S are 6.5*10-15 M and 1.8*10-22 M resp. (which corresponds with the solubility product of 10-50.12, but the dissolved concentrations of Ag and S are 7.1*10-15 M and 3.6*10-15 M resp., so for S seven order of magnitude higher. This is caused by the formation of S-complexes with protons (HS- and H2S (aq)) and to a lesser extent with Ag.

Complexation by organic matter

Complexation with Dissolved Organic Carbon (DOC) is different from inorganic complexation or complexation with well-defined compounds such as acetate or NTA. The reasons for that difference are as follows.

- DOC is very heterogeneous; DOC isolated at two sites may be very different (not to mention the difficulty of selecting isolation procedures).
- Complexation with DOC generally shows a continous range of equilibrium constants, due to chemical and steric differences in neighbouring groups.
- Increased cation binding and/or the ionic strength of the solution change electrostatic interactions among the functional groups in DOC-molecules, which influences the equilibrium constants.
- In addition, changing electrostatic interactions may cause conformational changes of the molecules.

Among the most popular models to assess organic complexation are Model V (1992), VI (1998) and VII (2011), also known as WHAM, written by Tipping and co-authors (Tipping & Hurley, 1992; Tipping, 1994, 1998; Tipping, Lofts & Sonke, 2011). All these models assume that two types of binding occur: specific binding and accumulation in the diffuse double layer. Specific binding is the formation of a chemical bond between an ion and a functional group (or groups) on the organic molecule. Diffuse double layer accumulation is the accumulation of ions of opposite electrical charge adjacent to the molecule, without formation of a chemical bond (the electrical charge is usually negative, so the ions that accumulate are cations).

For specific binding, all these models distinguish fulvic acids (FA) and humic acids (HA) which are treated separately. These two classes of DOC are typically the most abundant components of natural organic matter in the environment – in surface freshwaters, the fulvic acids are typically the most abundant. For each class, eight different discrete binding sites are used in the model. The sites have a range of acid-base properties. Metals bind to these sites, either to one site alone (monodentate), to two sites (bidentate) or, starting with Model VI, to three (tridentate). A fraction of the sites is allowed to form bidentate complexes. Starting with Model VI, for each bidentate and tridentate group three sub-groups are assumed to be present – this further increases the range of metal binding strengths.

Binding constants depend on ionic strength and electrostatic interactions. Conditional constants are calculated in the same way in Model V, VI and VII, as follows:

where:

- Z is the charge of the organic acid (in moles per gram organic matter);
- w is calculated by :

where:

- P is a constant term (different for FA and HA, and different for each model);
- I is the ionic strength.

Therefore, the conditional constant depends on the charge on the organic acids as well as on the ionic strength. For the binding of metals, the calculation of the conditional constant occurs in a similar way.

The diffuse double layer is usually negatively charged, so it is usually populated by cations, in order to maintain electric neutrality. Calculations for the diffuse double layer are the same for Model V, Model VI and Model VII. The volume of the diffuse double layer is calculated separately for each type of acid, as follows:

where:

- NAv is Avogadro's number;
- M is the molecular weight of the acid;
- r is the radius of the molecule (0.8 nm for fulvic acids, 1.72 for humic acids);
- κ is the Debye-Hückel parameter, which is dependent on ionic strength.

Simply applying this formula in situations of low ionic strength and high content of organic acid would lead to artifacts (where volume of diffuse layer can be calculated to be more than 1 liter/liter). Therefore, some "tricks" are implemented to limit the volume of the diffuse double layer to 25% of the total.

In case the acid has a negative charge (as it has in most cases), positive and neutral species are allowed to enter the diffuse double layer, just enough to make the diffuse double layer electrically neutral. When the acid has a positive charge, negative and neutral species are present.

The concentration of species in the diffuse double layer is calculated by assuming that the concentration of that species in the diffuse double layer depends on the concentration in the bulk solution and the charge.

In formula:

where R is calculated iteratively, to ensure the diffuse double layer is electrically neutral.

Applications

Speciation models can be used for many purposes. Basically, two groups of applications can be distinguished. The first group consists of applications meant to understand the chemical behaviour of any system. The second group focuses on bioavailability.

Chemical behaviour; laboratory situations

Speciation models can be helpful in understanding chemical behaviour in either laboratory situations or field situations. For instance, if you want to add EDTA to a solution to prevent metals from precipitation, the choice of the EDTA-substance also determines the pH of the final solution. Figure 1 shows the pH of a 1 mM solution of EDTA for five different EDTA-salts. This shows that if you want to end up with a near neutral solution, the best choice is to add EDTA as the Na3HEDTA-salt. Adding a different salt requires adding either acid or base, or more buffer capacity, which in turn will influence the chemical behaviour of the solution.

Figure 1. pH of a 1 mM EDTA-solution for different EDTA-salts. Data obtained using speciation program CHEAQS Next.

If you have field measurements of redox potential, speciation models can help to predict whether iron will be present as Fe(II) or Fe(III), which is important because Fe(II) behaves quite different chemically than Fe(III) and also has a quite different bioavailability. The same holds for other elements that undergo redox equilibria like N, S, Cu or Mn.

Phase reactions can be predicted with speciation models, for example the dissolution of carbonate due to the gas solution reaction of CO2. Another example is the speciation in Dutch Standard Water (DSW), a frequently used test medium for ecotoxicological experiments, which is oversaturated with respect to CaCO3 and therefore displays a part of Ca as a precipitate. The fraction that precipitates is very small (less than 2% of the Ca) so it seems unimportant at first glance, but the precipitate induces a pH-shift of 0.22, a factor of almost two in the concentration of free H+.

Many metals are amphoteric and have therefore a minimum solubility at a moderate pH, while dissolving more at both higher and lower pH-values. This can easily be seen in the case of Al: Figure 2 shows the concentration of dissolved Al as a function of pH (note log-scale for Y-axis). Around pH of 6.2, the solubility is at its minimum. At higher and lower pH-values, the solubility is (much) higher.

Figure 2. Soluble Al as function of pH. Data obtained using speciation program CHEAQS Next.

Speciation models can also help to understand differences in the growth of organisms or adverse effects on organisms, in different chemical solutions. For example, Figure 3 shows that changes in speciation of boron can be expected only between roughly pH 8 and 10.5, so when you observe a biological difference between pH 7 and 8, it is not likely that boron is the cause. Copper on the other hand (see Figure 4) does display differences in speciation between pH 7 and 8 so is a more likely cause of different biological behaviour.

Figure 3. Speciation of B as function of pH; concentration of boron was 1x10-6 M. At higher concentrations, complexes with 2, 3, 4 or 5 B-ions can be formed at significant concentrations. Data obtained using speciation program CHEAQS Next.

Figure 4. Speciation of copper(II) as a function of pH; concentration of copper was 3x10-8 M. At higher concentrations, complexes with 2 or 3 Cu(II)-ions can be formed at significant concentrations. Data obtained using speciation program CHEAQS Next.

Chemical behaviour: field situations

In field situations, the chemistry is usually much more complex than under laboratory conditions. Decomposition of organisms (including plants) results in a huge variety of organic compounds like fulvic acids, humic acids, proteins, amino acids, carbohydrates, etc. Many of these compounds interact strongly with cations, some also with anions or uncharged molecules. In addition, metals easily adsorb to clay and sand particles that are found everywhere in nature. To make it more complex, suspended matter can contain a high content of organic material which is also capable of binding cations.

For complexation by fulvic and humic acids, Tipping and co-workers have developed a unifying model (Tipping & Hurley, 1992; Tipping, 1994, 1998; Tipping, Lofts & Sonke, 2011). The most recent version, WHAM 7 (Tipping, Lofts & Sonke, 2011), is able to predict cation complexation by fulvic acids and humic acids over a wide range of chemical circumstances, despite the large difference in composition of these acids. This model is now incorporated in several speciation programs.

Suspended matter may be of organic or of inorganic character. Inorganic matter usually consists of (hydr)oxides of metals, such as Mn, Fe, Al, Si or Ti, and clay minerals. In practice, the (hydr)oxides and clays occur together, but the mutual proportions may differ dramatically depending on the source. Since the chemical properties of these metal (hydr)oxides clays are quite different, there is a huge variation in the chemical properties of inorganic suspended matter in different places and different times. As a consequence, modeling interactions between dissolved constituents and suspended inorganic matter is challenging. Only by measuring some properties of suspended inorganic matter, can modeling be applied successfully. For suspended organic matter, the variation in properties is also large and modelling is challenging.

Bioavailability

Speciation models are useful in understanding and assessing the bioavailability of metals and other elements in test media. Test media often contain substances like EDTA to keep metals in solution. EDTA-complexes in general are not bioavailable, so in addition to keeping metals in solution they also change their bioavailability. Models can calculate the speciation and help you to assess what is actually happening in a test medium. An often forgotten aspect is the influence of CO2. CO2 from the ambient atmosphere can enter a solution or carbonate in solution (if in excess over the equilibrium concentration) can escape to the atmosphere. The degree to which this exchange takes place, influences the pH of the solution as well as the amount of carbonate that stays in solution (carbonates are often poorly soluble).

Similarly, in field situations models can help to understand the bioavailability of elements. As stated above, the influence of DOC can nowadays be assessed properly in many situations, the influence of suspended matter remains more difficult to assess. Nevertheless models can deliver insights in seconds that otherwise can be obtained only with great difficulty.

Models

There are many speciation programs available and several of them are freely available. Usually they take a set of total concentrations as input, plus information about parameters such as pH, redox, concentration of organic carbon etc. Then the programs calculate the speciation and present them to the user. The equations cannot be solved analytically, so an iterative procedure is required. Although different numerical approaches are used, most programs construct a set of non-linear mass balance equations and solve them by simple or advanced mathematics. A complication in this procedure is that the equilibrium constants depend on the ionic strength of the solution, and that this ionic strength can only be calculated when the speciation is known. The same holds for the precipitation of solids. The procedure is shown in Figure 5.

Figure 5. Typical flow diagram of a speciation program.

Limitations

For modeling speciation, thermodynamic data is needed for all relevant equilibrium reactions. For many equilibria, this information is available, but not for all. This hampers the usefulness of speciation modeling. In addition, there can be large variations in the thermodynamic values found in the literature, resulting in uncertainty about the correct value. A factor of 10 between the highest and lowest values found is not an exception. This of course influences the reliability of speciation calculations. For many equilibria, the thermodynamic data is only available for the standard temperature of 25°C and no information is available to assess the data at other temperatures, although the effect of temperature can be quite strong. Also ionic strength has a high impact on equilibrium ‘constants’; there are many methods available to correct for the effect of ionic strength, but most of them are at best semi-empirical. Simonin (2017) recently proposed a method with a solid theoretical basis; however, the data required for his method are available only for a few complexes so far.

More fundamentally, you should realize that speciation programs typically calculate the equilibrium situation, while some reactions are very slow and, more inportant, nature is in fact a very dynamic system and therefore never in equilibrium. If a system is close to equilibrium, speciation programs can often make a good assessment of the actual situation, but the more dynamic a system, the more care you should take in believing the programs’ results. Nevertheless it is good to realise that a chemical system will always move towards the equilibrium situation, while organisms may move them away from equilibrium. Phototrophs are able to move a system away from its equilibrium situation whereas decomposers and heterotrophs generally help to move a system towards its equilibrium state.

References

Simonin , J.-P. (2017). Thermodynamic consistency in the modeling of speciation in self-complexing electrolytes. Ind. Eng. Chem. Res. 56, 9721-9733.

Tipping, E., Hurley, M.A. (1992). A unifying model of cation binding by humic substances. Geochimica et Cosmochimica Acta 56, 3627 - 3641.

Tipping, E. (1994). WHAM - A chemical equilibrium model and computer code for waters, sediments, and soils incorporating a discrete site/electrostatic model of ion-binding by humic substances. Computers & Geosciences 20, 973 - 1023.

Tipping, E. (1998). Humic Ion-Binding Model VI: An Improved Description of the Interactions of Protons and Metal Ions with Humic Substances. Aquatic Geochemistry 4, 3 - 48.

Tipping, E., Lofts, S., Sonke, J.E. (2011). Humic Ion-Binding Model VII: a revised parameterisation of cation-binding by humic substances. Environmental Chemistry 8, 228 - 235.

Further reading

Stumm, W., Morgan, J.J. (1981). Aquatic chemistry. John Wiley & Sons, New York.

More, F.M.M., Hering, J.G. (1993). Principles and Applications of Aquatic Chemistry. John Wiley & Sons, New York.