Metal speciation from stream to open ocean: modelling v. measurement

It has long been recognised that the bioavailability, transport and retention of metals in the natural environment depend very much on their chemical speciation, by which is meant the different chemical species that make up the total amount of metal in a given sample.^[1] Bioavailability underpins the role of metals as essential nutrients and governs toxicity, and is considered to depend on the concentration of the free metal ion, combined with those of other competitors including protons.^[2] The key process of metal partitioning between solutions and mobile or immobile solids depends on speciation in both phases.^[3] A well-known example of the importance of metal speciation to ecosystem function is the case of freshwater acidification,^[4] echoed in more recent concern about the acidifying effect of excess carbon dioxide on marine waters.^[5,6]

Chemical speciation can be determined analytically and, in principle at least, this is the most reliable approach. However, direct measurements require considerable effort, either in terms of sampling and laboratory analysis^[7] or in the development and maintenance of in situ techniques,^[8] and of course cannot provide predictions of future conditions. To address these limitations, therefore, chemical speciation modelling is an essential complementary activity, offering the possibility to address ‘what if’ questions, and explain current conditions and forecast future ones, at different spatial scales. It is clearly desirable that measurements and model predictions should agree, and here we explore how well they do so. The analysis considers systems assumed to be at or near to chemical equilibrium, and thereby amenable to approximate but practical analysis using the internally consistent rules of equilibrium chemistry that follow from the laws of thermodynamics.

In aquatic systems, a dissolved cationic metal is present as the free metal ion, its complexes with inorganic ligands such as chloride and bicarbonate ions, and its complexes with organic compounds, either identifiable ones such as acetate or siderophores, or the complex mixture of partial breakdown products of living matter. How metals bind to the inorganic and identifiable organic ligands is, or can be, well understood, but interactions with the ‘natural organic matter’, ‘dissolved organic matter’ (DOM) or ‘humic substances’ are often quantitatively more important as repositories of the cations. Despite decades of research on this complex organic material, it remains poorly defined chemically. However, we know that it contains weak-acid functional groups, including carboxyl and phenolic groups, and entities containing N and S, and these provide a range of binding sites for cations.^[9] Models have been developed that take into account this high degree of binding site heterogeneity, the most successful and widely used ones to date being the Windermere Humic Aqueous Model (WHAM)^[10] and the Non-Ideal Competitive Adsorption–Donnan (NICA–Donnan) model,^[11] both of which have been parameterised with laboratory data. These models take into account binding at different metal loadings and the effects of pH, competition by major cations and ionic strength, and so are suitable for field application. However, they have been parameterised with data for isolated natural organic matter (chiefly freshwater and soil humic and fulvic acids) obtained by experiments in which the analyst could impose the solution or suspension conditions and thereby optimise the measurement of the metal free ion or related variables. It is much more difficult to measure chemical speciation on natural samples, even more difficult to do so in situ, but this is what we really want to know.

Some efforts have been made to compare techniques applied to common samples,^[12,13] and to compare results with model predictions, but progress has been limited, not least by the effort required and the time limitations of an individual study. A less direct but more widely applicable approach is the comparison of analytical results with outputs from the same model, and this also provides a test of the model. We have previously performed three such studies focussing on data from fresh waters,^[14] estuarine and coastal waters^[15] and the open ocean,^[16] and comparing them with predictions obtained using WHAM6^[17] or WHAM7.^[18] We found that measurements and model predictions based only on metal binding by humic-type ligands (i.e. DOM) were in broad agreement, but in many cases differences were greater than could be explained by data uncertainty, and there was a tendency for the differences to be larger at low metal concentrations.

The work described here extends in two ways our previous efforts to compare observations and predictions of metal speciation in aquatic systems. First, by combining the results of the three previous studies, we obtained a large data set (~2000 samples) covering all types of surface water, including the deep ocean, and with some overlap of techniques applied to different types of water. Second, we improved the quantification of comparison between measurements and predictions. As before, we compared results in terms of the free-ion concentration [M] and DOM-bound metal (ν, mol g^–1), but we added their ratio ν/[M] (L g^–1), which is a kind of ‘local’ partition coefficient. We combined output statistics to create a single variable that quantifies the degree of agreement between measurement and prediction. Then we compared measured and modelled values by surface-water type, by metal and by measurement technique. This was an empirical comparison, with the simple aim of quantifying differences. We anticipated that analysis of the large combined data set would increase the possibility of identifying underlying trends and discrepancies, and therefore help to define future directions in efforts to understand and predict metal speciation in the field.

The measured speciation data, taken from the literature, were obtained by eight different techniques, briefly described below. They are each designed to measure (quasi-) equilibrium speciation in terms of conventional chemistry, based on the concentration of the free metal ion.

Cation-exchange, dynamic mode (CED) was introduced by Driscoll^[19] to determine Al speciation in acid to neutral fresh waters. It operates on the principle that most inorganic Al is cationic, and can be rapidly removed during passage through a column of cation-exchange resin (labile fraction), whereas Al complexed by organic matter remains unaffected (non-labile). A correction to account for the minor dissociation of the organic complex can be applied.^[20] The free-ion (Al³⁺) concentration is not measured, but calculated by equilibrium speciation modelling of the inorganic forms only.

Competing ligand methods (CL-ASV, CL-CSV) involve the titration of multiple subsamples (to which a buffer and varying concentrations of a metal spike have been added) with a strongly complexing ligand (CL) as a means of converting the electrochemically inert fraction of metal into a single, well-characterised complex that is reducible and thus can be subsequently measured by adsorptive stripping voltammetry.^[21] Anodic stripping voltammetry (ASV) uses a negative voltage for collection of the CL–metal complexes at the working electrode,^[22] whereas cathodic stripping voltammetry (CSV) uses a zero or positive value.^[23] Data from multiple experiments are then used to determine the conditional stability constant and complexation capacities of the organic ligand using non-linear curve fitting. Calculations require known side reaction coefficients for both the competing ligand and the inorganic metal.^[24]

The Donnan membrane technique (DMT) involves the separation of the in situ solution being analysed from an acceptor solution by a cation-exchange membrane.^[25] The acceptor solution contains Ca(NO₃)₂ and humic acid. The DMT cell is immersed into the river or lake for a set number of days, during which time free metal ions diffuse across the cation-exchange membrane and are complexed by the humic acid. After equilibration, metal concentrations in the acceptor solution are measured, and the free-ion concentration is computed by speciation modelling. If equilibrium is not achieved, the free-ion concentration is computed assuming diffusion-limited transport of the free metal ion across the cation-exchange membrane.

Differential pulse anodic stripping voltammetry (DPASV) uses the electrochemical technique used in CL-ASV but without the addition of a competing ligand. Samples are titrated with metal between voltammetric scans to estimate concentrations of labile metal, from which free-ion concentrations are calculated with inorganic side reaction coefficients, and then a linear transformation model is applied to determine conditional stability constants for one or two classes of ligands.^[26]

In the ion-exchange column technique (IET),^[27,28] a column of sulphonic acid-type resin is initially equilibrated with an electrolyte solution containing defined concentrations of Na, Mg, K and Ca. The resin is then calibrated by equilibrating it with solutions containing defined concentrations of the trace metal(s) of interest at a range of pH and a fixed ionic strength. A distribution coefficient, defining the ratio of resin-adsorbed metal to free metal ion, is calculated for the pH of the field water for which free-ion measurements are desired. The precalibrated column is then equilibrated with the field water of interest and the free metal ion concentrations computed using the distribution coefficient and the measured amounts of metal bound to the resin. Prior to equilibration with the column, the ionic strength of the field water is corrected to that of the calibration solutions using a multi-electrolyte stock containing Na, Mg, K and Ca nitrates.

In the ion-selective electrode (ISE) method, the electrode potential between an ion-selective electrode and a reference electrode is used for the direct determination of free copper ion concentration in a flow-through system.^[29]

Permeable liquid membrane (PLM) techniques^[30,31] are based on the measurement of metal flux, as generated by transport via a carrier molecule, through a hydrophobic membrane. The latter is sandwiched between a sample solution and a strip of receiving solution that contains a strong chelating agent to allow preconcentration. The metal of interest accumulates in the strip solution over time, the flux being directly related to either the concentration of free ion or to the free ion and labile metal complexes depending on the experimental setup.

The data set of measured speciation variables, together with the necessary additional information required as model inputs (pH, concentrations of major ions and DOM), was initially compiled from results for 2088 different samples. However, we restricted our analysis here to results that were realistic and complete, i.e. only accepting samples for which the values of ν were positive. This meant rejecting a total of 66 samples (3 %). The data set is summarised in Tables 1 and S1, with further details in references.^[14–16]

Table 1.  Numbers of samples used in the analysis
FW, fresh water; E&C, estuarine and coastal waters; OO, open ocean. See text for key to methods and abbreviation definitions

We used WHAM^[10] incorporating Humic Ion-Binding Model VII^[18] to perform the speciation calculations; previously, the open-ocean calculations^[16] were done using Model VI.^[17] Models VI and VII use structured formulations of discrete, chemically plausible binding sites for protons in humic and fulvic acids (HA, FA), in order to allow the creation of regular arrays of bidentate and tridentate binding sites for metals. Metal aquo ions (Al³⁺, Cu²⁺, Cd²⁺, etc.) and their first hydrolysis products (AlOH²⁺, CuOH⁺, CdOH⁺, etc.) compete with each other, and with protons, for binding. The same intrinsic equilibrium constant for binding to carboxyl or type A groups (K_MA) is assumed to apply to the aquo ion and its first hydrolysis product. The constant for binding to weaker acid groups (K_MB) is related to K_MA, and the contributions of rarer ‘soft’ ligand atoms are factored in using a correlation with equilibrium constants for metal binding by NH₃.^[9,17] The intrinsic equilibrium constants are modified by empirical electrostatic terms that take into account the attractive or repulsive interactions between ions and the charged macromolecule. WHAM constants are derived from the results of numerous studies of proton and metal binding by isolated humic substances, together with linear free-energy relationships.^[18] To make WHAM7, the humic ion-binding model is combined with an inorganic speciation model, the species list and constants for which were given by Tipping.^[10] The inorganic reactions in this database are restricted to monomeric complexes of metals.

Temperature effects on reactions between inorganic species are taken into account using published or estimated enthalpy data, but in the absence of sufficient experimental information, reactions involving humic substances are assumed to be independent of temperature.

The effects of ionic strength on the inorganic reactions are taken into account using the extended Debye–Hückel equation, which relates the charge on an ion to an activity coefficient at a specified ionic strength. This approach is generally accepted for the low ionic strengths of fresh waters,^[32] but in the modelling of high-ionic-strength systems, such as estuarine and marine systems, it is usually recommended to use the ‘mean salt method’ for calculating activity coefficients, rather than the methods normally applied to lower-ionic-strength fresh water systems. The ion-pairing model described by Millero and Schreiber^[33] incorporates the mean salt method with Pitzer’s equations. We previously^[15] compared the results from ion-pairing with the extended Debye–Hückel (EDH) equation, and found that for the great majority of the ions considered in the present work, differences in the free-ion activities calculated using the two methods were within ±16 %. An exception was nitrate, where the free ion activity was 24 % higher when the EDH method was employed. These differences are small in comparison with the variations when comparing modelled and measured values of trace metal free-ion concentrations, and so the use of the extended Debye–Hückel equation for all waters is justified, and consistent with several other studies that have used component-independent relationships for activities in marine systems.^[34,35]

To run the model, we used measured total solute concentrations and pH. For all the freshwater samples, dissolved organic carbon (DOC) concentrations were available. For estuarine and coastal waters, some missing DOC concentrations had to be estimated from other studies carried out at the same location,^[15] and for open ocean samples, [DOC] was taken to be the value at the nearest location for which data were available.^[16]

The binding activity of DOM was estimated by assuming DOM to be 50 % carbon, and that 65 % of the DOM behaves like isolated FA whereas the rest is inert.^[36,37] For example, a DOC concentration of 5 mg L^–1 corresponds to [DOM] of 10 mg L^–1, and so the concentration of FA for modelling is 6.5 mg L^–1. All the predictions reported here were made using only humic-type DOM, i.e. no additional calculations were performed assuming the presence of anthropogenic ligands, e.g. ethylenediamine tetraacetic acid (EDTA).

For fresh waters from the field, we estimated truly dissolved Fe^III concentrations with the empirical equation of Lofts et al.,^[38] suitably modified for Humic Binding Model VII.^[14] For samples where Fe was not the metal of interest with respect to speciation, the dissolved Fe^III concentration was set to 10 nmol L^–1 for estuarine and coastal sites and to 1 nmol L^–1 for open-ocean sites. These procedures were adopted to take into account the significant competition by Fe^III for binding at the stronger sites.^[39]

We used three variables to compare measurements and model predictions. The ‘master’ variable is the free metal-ion concentration (Al³⁺, Cu²⁺, etc) denoted by [M] with units of moles per litre. The variable ν is used to quantify metal associated with DOM, with unit of moles per gram. A combination of these, ν/[M], with units of litres per gram, is useful because it permits comparisons of both strongly organically bound metals, for which the great majority of the metal is bound so that ν is well predicted, and the concern is about [M], and weakly bound ones, for which [M] is well predicted and the concern is about ν.

The expected behaviour of ν/[M] for the highly heterogeneous binding sites of DOM can be understood with reference to the theoretical calculations of Fig. 1. If DOM possessed only one type of binding site, with no electrostatic interactions, binding behaviour would be as shown in the upper three panels of Fig. 1. At low [M], ν increases such that the log–log plot is linear with a slope of unity, and log (ν/[M]) is constant, and equal to the equilibrium constant, when plotted against either log [M] or log ν. When the metal’s occupation of the sites passes ~10 % of the total (log ν ~ –3.3 in Fig. 1), log (ν/[M]) starts to decrease, and an increase in [M] causes a relatively smaller increase in ν. In other words, the affinity of the sites for metals effectively declines, owing to the fact that they are occupied already. We therefore can think of ν/[M] as a ‘local’ partition coefficient. The same behaviour is observed in the presence of competing metals (or the proton), which act to alter the effective equilibrium constant for binding.^[9] For DOM with heterogeneous binding sites, the relationships are less straightforward (Fig. 1, lower panels). The slope of log ν against log [M] is no longer unity, because the system passes through sites with decreasing affinity, and the different metals display marked differences with respect to binding site heterogeneity, from the low heterogeneity of Zn to the high of Hg. Nonetheless, ν/[M] continues to operate as a ‘local’ partition coefficient, providing a measure of the propensity of the DOM to bind more metal.

Fig. 1.  Explanation of ν/[M], where ν is bound metal per unit mass of dissolved organic matter (DOM) and [M] is the free metal ion concentration. The upper three panels show theoretical results for a single site ligand with log K (L mol^–1) = 3 (full line) or 6 (dashed line). The DOM has 5 × 10^–3 mol g^–1 of metal-binding sites for the metal, so that the maximum log ν value is –2.3. The lower panels are predictions with WHAM7 for the binding of four metals to DOM at I = 0.1 M and pH 7.

We made comparisons of logarithmic values of the variables, (1) in order to be able to cover the large ranges, (2) because comparisons are better represented by ratios of observed to calculated values, rather than absolute differences, which are biased towards high values, and (3) because the data are noisy.

Measured and model-predicted values were compared by plotting log values and deriving regression slopes, forced through (0,0), and values of r². Ideally the slopes and r² values for all three variables would be close to unity. However, their values in practice inevitably depend on the ranges of the variables, and so they do not tell the full story. Therefore, we used as the main criteria of agreement (1) the root-mean-squared deviations (RMSD) between measured and predicted values, and (2) the average deviations (AD) between measured and predicted. Ideally, both of these should be zero. Because appreciation and discussion of these four measures (each applying to three variables) would be involved and complex, we used a single measure based on RMSD and AD. We refer to this as ΣRA, defined as:

where the vertical lines indicate absolute values. Division by 6 simply reduces the total value so that it is more like the individual ones. Ideally, ΣRA should be zero. Because deviations are based on differences in log values, values of 0.5, 1 and 2 indicate overall agreement to approximately a factor of three, to one order of magnitude and to two orders of magnitude respectively.

The WHAM model incorporating Humic Ion Binding Model VII provides a consistent means of calculating the speciation of dissolved metals in different surface waters, taking into account competition effects, including the effects of pH and ionic strength. The equilibrium and other constants used in the model are based on many laboratory measurements with isolated humic substances,^[18] and therefore have not been calibrated with any field data. The factor of 65 % used to convert from DOM to FA (see above) is based on different field data from those considered here. Therefore, our speciation predictions are truly independent of the measured speciation variables, and the comparisons therefore provide an unbiased evaluation of the agreement or otherwise between measurements and predictions. In essence, we are taking information that many workers have obtained over many years from laboratory experiments, as encapsulated within WHAM, and comparing it with field observations.

The results cover the three major types of surface water, i.e. fresh, estuarine–coastal and the open ocean (Table 1). Therefore, they include wide ranges of ionic strength and of the concentrations of major metals (Mg, Ca) that compete significantly for binding to DOM. The DOM itself must vary in its composition and properties, given its different sources (terrestrial plants, freshwater and marine phytoplankton, microbes, etc). From Table S1, we see that whereas the estuarine–coastal and open-ocean samples refer to quite narrow ranges of pH, the range for fresh waters is wider, although acid systems are poorly represented, except in the Al data set. We have data for nine different metals, and these reflect research interest together with suitability for analysis. Some metals are better represented than others; for example, Al only appears in fresh waters, whereas most of the Fe measurements are for the open ocean. Copper is the most-analysed metal (35 % of all results) and covers all three of the surface-water types. To test model predictions and compare different analytical techniques against the model, the ideal data set would provide even coverage of conditions, a wide range of metals with similar amounts of data for each aquatic environment, and analysis by several methods on the same or similar samples common samples. Clearly, the present data set does not fully meet these characteristics, but it is large and wide-ranging and definitely useful for a meta-analysis.

All data are plotted in Fig. 2, and statistical comparisons for the different types of surface water and the combined data set are given in Table 2. The results in Fig. 2 show strong correlations between predicted and measured values (high R²), which suggests that at the large scale, the general chemistry of metals in natural waters follows the expected trend from the laboratory-based model data. The regression depends considerably on data for Hg and Fe, both of which are very strongly bound (high ν/[M]), whereas the other metals show less variation at this scale.

Fig. 2.  Measured and predicted variables (ν, [M], ν/[M], see text for explanation) for all metals in fresh waters (circles), estuarine and coastal waters (triangles) and the open ocean (squares). The 1 : 1 lines are shown, together with one-log-unit ranges.

Table 2.  Summary of correlations and errors for all metals
The slopes were obtained by plotting log predicted v. log measured values, forced through the origin; RMSD, root-mean-squared error in residuals between measured and predicted values; AD, average residual (measured – predicted); ΣRA is defined by Eqn 1. Variables: [M], free ion concentration; ν, moles of metal bound per gram of dissolved organic matter

Regression slope, RMSD and AD values show the same pattern for each surface-water type and for the combined data set (Table 2). The log [M] slopes are all fairly close to unity, but consistently too low, whereas the log ν slopes are consistently slightly greater than unity, and the log ν/[M] slopes are consistently less than unity. The AD values follow the pattern of the slope values, which means that predictions of log [M] tend to be greater, and those of log ν and log ν/[M] smaller, than the measured values.

Comparing the three types of surface water by the variable ΣRA, it is found that the best agreement between predictions and measurements is found for fresh waters (ΣRA = 0.52), followed by estuarine–coastal systems (ΣRA = 0.77) and then the open ocean (ΣRA = 1.24). The differences reflect increasing noise in the data rather than a systematic change in the relationship between measured and predicted values. This may simply reflect the greater analytical challenge of making speciation measurements in marine and marine-influenced waters.

For fresh waters, Lofts and Tipping^[14] calculated how much error could be associated with the model predictions on the basis of uncertainties in input data, model parameters and variations in DOM binding properties. They then added uncertainty from the analytical precision of the speciation measurement. Only for Al, analysed by the CED method, could differences between measurements and predictions be accounted for by these factors, and it was concluded that there must be additional reasons for discrepancies in predicted and measured values. The same conclusion is likely to apply to results for estuarine and coastal systems and the open ocean, given that they show worse agreements than fresh waters. One possible reason for differences greater than those attributable to data uncertainty and DOM variations emerged in our analysis of data from estuarine and coastal waters^[15] and is the presence of anthropogenic ligands, such as EDTA, which are found in the receiving waters of many industrial areas and can be extremely persistent in wastewater treatment plants and natural waters.^[40] Better agreement with observed speciation was obtained when reasonable concentrations of EDTA were included in the model inputs (see also Baken et al.^[41] and Ahmed et al.^[42]). However, enhanced organic complexation due to anthropogenic ligands will apply only to a few of the samples considered here. Another possibility is that some natural waters contain high-affinity ligands that are not represented in humic-type DOM, as suggested in particular for marine waters.^[43,44] Nonetheless, it is likely that in many or most cases, differences between predicted and measured speciation variables arise from faults in the measurement techniques, or in the model, or in both.

Fig. 3 compares results for the nine different metals for which we have data. The results for Al are good, helped by its relatively high concentrations, compared to those of the other metals, and the low pH values, which lead to readily measurable and predictable values of [M], ν and ν/[M]. In contrast, for Fe^III total concentrations are rather low, pH values neutral or higher, and the samples are from estuarine and coastal and open ocean waters, and so [M] and ν values are much lower than those for Al, and measurement and prediction are harder. The predictions of [M] tend to be too high, more so at lower measured values. Values of ν for Fe^III are well predicted because a very high proportion of truly dissolved Fe (i.e. not including oxide colloids) in natural waters is bound to organic matter, but because the [M] values are too high, values of ν/[M] are too low.

Fig. 3.  Measured and predicted variables (ν, [M], ν/[M]. see text for explanation) for individual metals. The 1 : 1 lines are shown, together with one-log-unit ranges.

These trends for Fe^III are also seen for Cu. For Ni, Zn, Cd and Pb, agreements of both log [M] and log ν are good at higher concentrations, but predicted log [M] is too high and predicted log ν and log ν/[M] too low at lower concentrations. Agreements for Co^II are especially poor, perhaps because much of the metal was present as Co^III in marine samples.^[16] The variable ν/[M] is poorly predicted for Cd in the open ocean.

The binding of Hg by natural organic matter is extremely strong^[45] and consequently, there is almost perfect agreement between measured and predicted values of log ν (Fig. 3). This is the one case where predicted log [M] values are consistently lower than the measured ones, and consequently the model predicts higher values of log ν/[M]. Both measurements and predictions yield very high values of log ν/[M], of the order of 20, which is the highest of all the metals, ~5 log units greater than the values for the second most strongly bound metal, Fe^III.

A plot of the predicted log ν/[M] averaged over all data for each metal v. the average of the measured values (Fig. 4) shows a strong relationship, which indicates that the order of reactivity of the metals, and their consequent binding to organic matter in natural waters follow the expectations of laboratory experiments, as encapsulated in WHAM7. (Only Co deviates strongly from the trend, as discussed above.) This reinforces the more general agreements shown in Fig. 2, confirming the basic assumption of our modelling approach, that the chemical speciation of metals in natural surface waters arises from equilibrium metal–ligand interactions involving humic-type organic matter.

Fig. 4.  Measured v. predicted values of log ν/[M] (see text for explanation), averaged over all results for each metal. The 1 : 1 line is shown. Standard deviations scarcely exceed the areas of the points.

Fig. 5 summarises results obtained by different measurement methods. Table 3 compares results, including the variable ΣRA. The lowest ΣRA is produced by CED, because this method was used only for Al in fresh waters, providing favourable conditions for measurement, as discussed above. The next lowest ΣRA is given by PLM, but only on the basis of 19 points and 3 metals. Furthermore, the statistical variables for this method benefit from our removal of impossible results, i.e. negative values of ν.

Fig. 5.  Measured and predicted variables (ν, [M], ν/[M], see text for explanation) for different methods in fresh waters (circles), estuarine and coastal waters (triangles) and the open ocean (squares). The 1 : 1 lines are shown, together with one-log-unit ranges.

Table 3.  Summary of statistical data for different measurement methods
See text for key to abbreviations of measurement methods, [M] and ν as in Table 2

The CL-ASV, DMT and IET methods give similar values of ΣRA (Table 3). DMT and IET both benefit statistically from the removal of impossible results, and neither produces any systematic deviations; the error terms are mainly due to noise. Predicted and measured results obtained by the CL-ASV method are in fair agreement, although there is a tendency for predicted free-ion values to be higher than measured ones, whereas the opposite applies to ν and ν/[M].

The DPASV method has been used widely in all three types of surface water; there is better agreement for fresh waters, but in estuarine and coastal samples, and more so in those from the open ocean, the model calculates [M] too high and ν and ν/[M] too low.

The most commonly applied technique (34 % of all analyses considered here) has been CL-CSV. The model predictions of [M] tend to be either in agreement or too high (Fig. 3), predictions of ν and ν/[M] either in agreement or too low. Consequently, the results overall are biased, leading to a high value of ΣRA.

The ISE results only apply to estuaries, and the results show very poor agreement between measurements and predictions. Because dissolved Cu is largely in the organically bound state, agreement is good for log ν but the predicted values of log [M] are much lower than the measured values and this leads to poor agreement in log ν/[M]. Reasons for this were discussed by Stockdale et al.^[15] and the two most likely possibilities were first that the ‘dissolved’ Cu included significant amounts of colloidal material^[46] and so the input concentration for modelling is too high, as well as the measured value of ν, and second that anthropogenic ligands were present in the samples, leading to more binding of the metal than would occur with only humic-type DOM.

The competitive ligand approach (i.e. the CL-ASV and CL-CSV methods) has been criticised on the grounds that equilibration is unlikely,^[47] which would cause underestimation of [M], and overestimation of ν and ν/[M]; van Leeuwen and Town^[47] showed that errors of several orders of magnitude in ν/[M] could result. Therefore, it is instructive to compare the results obtained by CL-ASV and CL-CSV with those obtained by the other six techniques (Table 4). Taken overall, the non-CL methods yield unbiased results, the average differences between measured and predicted values of log [M], log ν and log ν/[M] being close to zero. This applies whether or not the ISE data are included, but the standard deviations are appreciably reduced when the ISE results are omitted. In contrast, for the CL methods, the average difference in log [M] is –1.5, i.e. the measured values are lower than predicted ones by a factor of 30, the difference in log ν is 0.5 (factor of 3) and the difference in log ν/[M] is 2.0 (factor of 100). Only the few results for Hg (Fig. 3) show the opposite trend. On this basis, results from the CL methods differ from those obtained by other methods.

Table 4.  Average differences between measurements and predictions, comparing competitive ligand (CL) methods and all other methods
If the value is positive, the measured value exceeds the predicted one and vice versa, [M] and ν as in Table 2

A key aspect of DOM in its interactions with metals is binding site heterogeneity, and for many natural waters, especially seawater, where metal concentrations are low, the strongest binding sites are of primary importance. Therefore, it is worthwhile to consider how differences between measured and predicted values vary with the extent to which the strong and weaker sites are occupied, which can be most straightforwardly done by plotting the differences against log ν (Fig. 6). When all the data are considered, we find significant trends with log ν, so that as log ν becomes increasingly negative, measured log [M] values become more negative than the predicted ones, whereas the opposite is true for log ν and log ν/[M]. But if only data for log ν > –6 are considered (right-hand panels of Fig. 6), although statistically significant variations occur, they are much smaller than at low values of ν, a finding that fits with results of Ahmed et al.^[48] For the data at log ν > –6 (838 samples), the value of ΣRA is quite small, at 0.34. Thus, there is little systematic bias and fair agreement between measured and predicted speciation variables at log ν > –6, whereas at lower loadings of the DOM, metal binding is stronger than predicted.

Fig. 6.  Differences between measured and predicted logarithmic variables (ν, [M], ν/[M], see text for explanation), based on measurements made by methods other than CL-ASV, CL-CSV and ISE (see text for explanation), plotted against measured values of ν (moles of metal bound per unit mass of dissolved organic matter). Left panels show all data (n = 1125), right panels data for ν > 10^–6 mol g^–1. All the slopes are significant (P < 0.001).

Model parameters for binding at the strong sites have been derived from laboratory studies in which the bias was towards higher loadings of the organic matter, i.e. high values of ν, largely because the determination of binding at high metal concentrations is easier and therefore within the capabilities of more researchers. Nonetheless, appreciable data at low values of ν were used by Tipping et al.^[18] in parameterising Humic Ion-Binding Model VII. Of the 220 data sets employed, 41 contained values <–6 and 16 contained values <–8, the most negative value being –15. These cover the ranges found by measurement (Figs 3, 5, 6), and so it does not appear that the trends in Fig. 6 simply arise because data at low values of ν were not used in the model parameterisation. Indeed, the successful fitting of the laboratory data^[18] means that the model’s formulation of strong binding sites present at low densities is adequate, i.e. there is presently no evidence for stronger binding sites in isolated humic substances that might explain the deviations between measurements and predictions in the field data analysed here. Therefore, if the measurements are accurate and unbiased, the trends in Fig. 6 would indicate the presence in field samples of ligands with stronger affinities than binding sites in isolated humic materials. These might include, for example, the metal-specific ligands thought to be produced by marine phytoplankton.^[43,44]

The variable ν is central to the WHAM-F_Tox model of metal mixture toxicity,^[49] in which the observed parallelism between WHAM-calculated binding of metals by HA and metal accumulation by aquatic organisms^[50,51] is exploited to obtain a measure of toxic exposure in acute laboratory experiments. Toxicity of the accumulated metal is quantified by the variable F_Tox defined as Σα_iν_i where the α_i are toxicity coefficients for the different metals. The highest value of α_i derived so far is 2.1 × 10⁶ g mol^–1 (for cadmium effects on trout),^[52] and toxic effects occur only when F_Tox exceeds ~2.0. Therefore, the minimum value of ν for toxicity is ~10^–6 mol g^–1. For less toxic metals, the required values of ν are much higher, of the order of 10^–3 mol g^–1 for Zn for example. Therefore, the range of ν over which measurements and model predictions are in reasonable agreement, i.e. >10^–6 mol g^–1 (Fig. 6), corresponds well to the range over which metal toxic effects are observed, implying that the WHAM7 model is appropriate for toxicity evaluation and prediction.

The comparisons of measured and predicted speciation variables presented here suggest that if we are to move closer to the reliable measurement and prediction of metal speciation in surface waters, three principal issues need to be resolved. The first is the apparent systematic difference between measurements made by the competitive ligand equilibrium techniques and those made by other methods. As noted above, there may theoretical reasons for this.^[47] A way forward would be a concerted effort to compare techniques, building on previous work and ideas,^[13,53–55] but covering a wider and more representative range of conditions and metals.

The second issue arises from our finding that model predictions differ most from measured values (obtained using techniques other than the CL-based ones) only at low values of ν (<10^–6 mol g^–1). The solution conditions under which such low loadings of the DOM occur are, of course, those that present the greatest analytical challenges, and therefore will produce the greatest uncertainties. Research is required to evaluate the analytical methodologies at low values of ν, for both natural water samples and laboratory-prepared solutions using isolated humic materials. The key requirement is to ascertain whether or not natural waters possess ligands with higher affinities for metals than humic-type DOM.

Third, the possible involvement in metal speciation of synthetic ligands such as EDTA in anthropogenically affected waters deserves more attention. By measuring their concentrations in surface-water samples for which speciation is determined analytically, and including them in the speciation calculations, the contributions of natural DOM and the synthetic ligands to metal speciation could be established.

The emphasis needs to be on analytical measurements, because ultimately, these are the basis of speciation knowledge, and further model developments based on existing data are unlikely to improve understanding or predictive capabilities. An alternative, more pragmatic, approach would be to accept the present field analytical data and use them to modify WHAM7, for example by creating additional strong binding sites, or altering competition effects, to make the results at low ν agree better. In our opinion, such modifications would be premature. We suggest that the best way forward is through coordinated analytical efforts and speciation method comparison, supported by speciation modelling.

A summary of data used in the analysis is available in Table S1 and the source references are available in Table S2 (see http://www.publish.csiro.au/?act=view_file&file_id=EN15111_AC.pdf).