Environ. Sci. Technol. 2003, 37, 1553-1562
Modeling Complexometric Titrations of Natural Water Samples R O B E R T J . M . H U D S O N , * ,† EDEN L. RUE,‡ AND KENNETH W. BRULAND‡ Department of Natural Resources and Environmental Sciences, W-503 Turner Hall, 1102 South Goodwin Avenue, University of Illinois, Urbana, Illinois 61801, and Department of Ocean Sciences, University of California, Santa Cruz, California 95064
Complexometric titrations are the primary source of metal speciation data for aquatic systems, yet their interpretation in waters containing humic and fulvic acids remains problematic. In particular, the accuracy of inferred ambient free metal ion concentrations and parameters quantifying metal complexation by natural ligands has been challenged because of the difficulties inherent in calibrating common analytical methods and in modeling the diverse array of ligands present. This work tests and applies a new method of modeling titration data that combines calibration of analytical sensitivity (S) and estimation of concentrations and stability constants for discrete natural ligand classes ([Li]T and Ki) into a single step using nonlinear regression and a new analytical solution to the one-metal/twoligand equilibrium problem. When applied to jointly model data from multiple titrations conducted at different analytical windows, it yields accurate estimates of S, [Li]T, Ki, and [Cu2+] plus Monte Carlo-based estimates of the uncertainty in [Cu2+]. Jointly modeling titration data at lowand high-analytical windows leads to an efficient adaptation of the recently proposed “overload” approach to calibrating ACSV/CLE measurements. Application of the method to published data sets yields model results with greater accuracy and precision than originally obtained. The discrete ligand-class model is also re-parametrized, using humic and fulvic acids, L1 class (K1 ) 1013 M-1), and strong ligands (LS) with KS . K1 as “natural components”. This approach suggests that Cu complexation in NW Mediterranean Sea water can be well represented as 0.8 ( 0.3/0.2 mg humic equiv/L, 13 ( 1 nM L1, and 2.5 ( 0.1 nM LS with [Cu]T ) 3 nM. In coastal seawater from Narragansett Bay, RI, Cu speciation can be modeled as 0.6 ( 0.1 mg humic equiv/L and 22 ( 1 nM L1 or ∼12 nM L1 and ∼9 nM LS, with [Cu]T ) 13 nM. In both waters, the large excess (∼10 nM) of high-affinity, Cu-binding ligands over [Cu]T results in low equilibrium [Cu2+] of 10-14.5(0.2 M and 10-13.3(0.4 M, respectively.
Introduction Natural waters contain a diverse array of ligands that complex metal ions, strongly affecting both their biogeochemical * Corresponding author phone: (217)333-7641; fax: (217)244-3219; e-mail:
[email protected]. † University of Illinois, Urbana. ‡ University of California, Santa Cruz. 10.1021/es025751a CCC: $25.00 Published on Web 03/18/2003
2003 American Chemical Society
cycles and their ecotoxicological impacts (1-4). Major and trace inorganic ligands (5, 6), synthetic chelators (7), humic and fulvic acids (8), colloids (9), and microbial chelators (10) are among the most important of these complexing agents. While the coordination equilibria of major inorganic and synthetic ligands can be modeled accurately using critically reviewed stability constants (5, 11), those of natural organic and trace inorganic ligands cannot since their concentrations and stability constants are difficult to characterize. For this reason, accurately quantifying the speciation of metals and properties of relevant ligands in surface water systems requires the direct analysis of environmental samples. Making such speciation measurements is an expensive and painstaking process requiring trace metal clean procedures for collecting and processing samples (12) along with sensitive methods of analysis (13, 14). Furthermore, in many cases the ambient concentrations of natural ligands and of aquo metal ions are not directly measurable. Rather, they are derived from complexometric titration data by modeling, a two-step process comprised of (i) calibration of the analytical method and (ii) estimation of speciation model parameters. While analysts have long been aware of the difficulties inherent in calibration, it has been argued that the subsequent parameter estimates are biased largely at high metal additions, making published estimates of ambient metal ion concentrations and strong ligand parameters accurate (15). However, some have recently re-emphasized the biases inherent in conventional calibration methods (1618), particularly for water samples containing humic and fulvic acids, and in resultant speciation models. This critique augments a recent challenge to the meaningfulness of the distinctions made between classes of ligands (19). Thus, the interpretation of this extensive and still growing body of speciation data remains controversial. This paper presents a new method of analyzing complexometric titration data that addresses these important issues. Its first major element is a new analytical solution for the two-ligand problem used in conjunction with the largely neglected approach of simultaneously estimating all model parameterssthe analytical sensitivity, ligand concentrations, and complex stability constantssby nonlinear regression (20). Using this tool, we demonstrate that collinearity (21) can cause parameter estimates and derived ambient free metal ion concentrations obtained by modeling single titration curves to be much more uncertain and biased than is generally realized. The second key element of the method is jointly modeling the data from multiple titrations conducted at different analytical windows, which solves the collinearity problem and permits speciation techniques to be accurately calibrated. Application of the method to published data sets yields improved fits, consistently larger (yet plausible) analytical sensitivities, and meaningful detection of ligand classes. Finally, we illustrate how one may adapt the method to partition metal complexation into contributions from strong ligands and background humic and fulvic acids.
Background and Theory We focus here on two voltammetric techniques that are frequently employed in complexometric analyses of surface waters: (i) competitive ligand equilibration-adsorptive cathodic stripping voltammetry (CLE-ACSV or CSV), currently the most sensitive electrochemical method for measuring complexation of metals such as copper and iron, and (ii) anodic stripping voltammetry (ASV). The CSV method of van den Berg (22) requires that the analyst amend the water sample with a well-characterized, “added ligand” (AL) that VOL. 37, NO. 8, 2003 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
1553
forms planar metal complexes capable of adsorbing to the Hg electrode (23). ASV methods detect metal complexes with weak natural or added ligands that are labile toward reduction on the electrode (24). For both techniques, the same basic coordination equilibria of the aquo metal ion (Mz+) must be considered: KIN
Mz+ + XIN 798 M′ βALx
Mz+ + xAL 798 M(AL)x
(1)
Ki
Mz+ + Li 798 MLi Mz+ and its complexes with the major inorganic ligands (XIN), including OH-, Cl-, and CO32-, constitute the aggregate M′ species, while its complexes with AL make up M(AL)x. Herein, the term “reference complexes” (MREF) refers to M′ together with M(AL)x, while “reference ligands” refers to XIN and AL. One or more discrete ligand classes (Li) are used to represent the broad spectrum of remaining “natural ligands” being titrated. Classes for which both a total concentration ([Li]T) and a conditional stability constant (Ki) can be estimated from the data with enough precision to be considered statistically significant are assigned a numerical index, i.e., i ) 1, ..., 3. In some cases, a ligand class with only one significant parameter, termed “weak” (LW) or “strong” (LS) ligands, are fitted to data at high or low [Mz+], respectively. By convention, KS . K1 > K2 > K3 > KW. Standard mass law and conservation equations (15, 17, 22, 25) relate the concentration of each metal complex to its conditional stability constant (Ki or βx), the total ligand concentration, and [Mz+] (all in molar units). The ratio of the concentration of any metal complex to [Mz+], also known as its R or side reaction coefficient, is useful as a comparative measure of its abundance (22). The R coefficient of the complex MLi is
RMLi ≡ [MLi]/[Mz+] ) Ki[Li]T/(1 + Ki[Mz+])
(2)
It follows that regardless of the number of ligands present in a sample, the total quantity of metal complexed by natural ligands ([∑MLi]) equals [Mz+] times the sum of RMLi over all i. Thus, [∑MLi] is a single-valued function of [Mz+] and the parameters describing each ligand class:
[ΣMLi] ) ν([Mz+], [LS]T, Ki, [Li]T, RMLw)
(3)
This relationship lies behind graphs of [Mz+] versus the total concentration of metal in inorganic and natural ligand complexes, denoted [MNAT] (26) or [M]/T (16) ([M]/T ≡ [∑MLi] + [M′]). Because of the variety of ligands and competing metals present in aquatic systems, the form of the function ν is not known a priori, and for the most part, the parameters Ki and [Li]T are viewed as empirical descriptors of ν rather than constants for any particular ligands (27). Continuous distributions of Ki are also commonly used (28-30). The goal of the analyst is to determine ν by measuring [Mz+] and [∑MLi] over a range of [M]/T. [Mz+] is proportional to the measured peak current (I) in voltammetric methods:
[Mz+] ) I/(R′S)
(4)
where the sensitivity (S ≡ I/[MREF]) is a calibration factor specific to the analytical method and the sample being titrated. Because the total concentration of metal ([M]T ) [MREF] + [∑MLi]) generally is kept much below [XIN] and [AL]T throughout a titration, XIN and AL remain largely uncomplexed, making the R coefficient of the reference complexes 1554
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 37, NO. 8, 2003
FIGURE 1. ASV/Cu titration of Narragansett Bay seawater (24). Data (b). Curves for S2 (s) and S2\SREF (- -) models. Tangent lines with S equal to (i) SSIC (‚‚‚), (ii) SICI (- ‚‚ -), (iii) SS2 (s), or (iv) SREF (- -) and ∑[Li]T of (i) 55 nM from eq 6, (ii) 86 nM (S2\SICI), (iii) 102 nM (S2), and (iv) 120 nM (S2\SREF). (R′) or “analytical window” (31) nearly constant. The definitions of R′ and [M]T require that [∑MLi] ) [M]T - I/S. “External calibration” entails the assumption that S equals the slope (SREF) of an I vs [M]T plot measured in a reference medium. Since proper reference media, such as UV-oxidized water, have the same salinity, [AL]T, and pH as the sample, they contain identical concentrations of reference ligands but are devoid of Li, surfactants, metals, and other electroactive substances that can affect S (13). Surfactants in particular can cause S to be lower than SREF by a factor φ (g0):
S ) SREF(1 - φ)
(5)
As there is currently no basis for predicting φ in natural water samples, various internal calibration or standardization procedures to estimate S have been devised (15, 17, 18, 22). “Simple” internal calibration involves determining the sensitivity (SSIC) by linear regression of I versus [M]T at the highest metal additions (16), where the natural ligands appear to have been titrated on the basis of an essentially linear response (Figure 1):
I|High[M]T ≈ SSIC([M]T - [ΣLi]SIC T )
(6)
Note that the total concentration of natural ligands SIC (22), although its ([∑Li]SIC T ) must also be fitted to obtain S consistency with subsequent estimates of [∑Li]T is rarely checked. Titrations of natural ligands are normally incomplete in samples containing humic and fulvic acids despite their apparently linear response at high metal additions (15). Completely titrating these weak ligands can require attaining [M]T that are very much greater than ambient and result in [MREF] that are beyond the linear range of the analytical system (13). Consequently, the true S is often the slope of the tangent to the titration curve at a point that lies beyond the range of the data, causing significant bias in SSIC (15, 17, 18). To compensate, some analysts have graphically estimated S by fitting theoretical curves to titration data using SSIC and SREF as lower and upper bounds, respectively (24). We refer to
this method as “internal calibration by inspection” and to S estimated in this way as SICI. This approach has been formalized in a manual procedure that involves (i) calculating SSIC, (ii) fitting K1 and [L1]T using a linearization, (iii) generating synthetic data, and (iv) re-estimating S before iteratively repeating (ii) (18). In this recursive internal calibration, SRIC is essentially a third parameter that is fitted along with K and [L]T. While either approach may improve on simple internal calibration, they each have shortcomings and are still only reliable for nearly complete titrations, [M]T . [∑Li]T (18), which can be as difficult to identify as to conduct. A modified CSV internal calibration procedure that exploits differences between the equilibration times of the reactions in eq 1 has also been employed (13, 32). At a titration’s end, one further metal addition is made and I is measured immediately following the 60-s deposition period. If MLi formation is too slow to be significant during the deposition, then the fit of the additional datum to the simple internal calibration line could test the completeness of a titration. While attaining equilibrium with some natural ligands may take several hours as compared to seconds or less for reference ligands (13, 23), other recent CSV studies of fulvic and humic acids have found that equilibrium occurred within 30 s (16) and that significant weak complexation in natural water samples occurred within 5 min (32). Furthermore, since to date the test has been conducted by inspection (33), the usefulness of this “nonequilibrium internal calibration” (SNIC) is difficult to evaluate objectively. Now it can be shown that the relative negative bias (S) in the sensitivity estimated using any particular method of internal calibration (Sˆ ) that is attributable to untitrated ligands is
ˆ )/S ) RMLw/(R′ + RMLw) S ≡ (S - S
(7)
where RMLw is the side reaction coefficient of complexes that remain undetected and are by definition “weak” relative to the “detected” (parametrized) ligands in the analytical window (R′) under consideration. [Mz+] calculated using Sˆ in eq 4 will be biased high by an amount 1/(1 - S ). Similarly, since undetected weak complexes have by definition not been included, calculated [∑MLi] will be biased low (17) if S > 0. Because these concentrations define ν, negative bias in Sˆ ultimately results in fitted ligand parameters (Ki and [Li]T) that underestimate complexation by natural ligands. However, when Sˆ < SREF, as is commonly observed, one cannot determine what combination of φ > 0 and S > 0 is the cause on the basis of data from a single titration unless the Li are completely titrated. Solutions to this dilemma lie in exploiting the R′ dependence of S in eq 7. Since RMLw is defined relative to R′, S can be made to approach zero if R′ is made to greatly exceed the combined maximum RMLi for all ligands, i.e., ∑(Ki[Li]T). This is the essence of Kogut and Voelker’s (16) “overload” titration, in which a linear fit of high R′ data is used to estimate the sensitivity (SOV) and a concentration of strong ligands ([LS]OV T ), as in eq 6. A true overload titration must appear linear with SOV ) SREF and/or SOV ) SSIC for the titration at the next lower R′. Since parameters for the natural ligands of interest can only be estimated from titrations at R′ near RMLi (31), S at these lower R′ are derived from SOV by assuming that φ is independent of [AL]T. If correct, ratios of S at different [AL]T should exhibit the same ratios as SREF at the same [AL]T (16), leading to
S([AL]T) ) S0RAL
(8)
where RAL ≡ SREF([AL]T)/SREF and the subscript “0” denotes 0 S at the highest [AL]T used in a data set. For common AL,
there is a broad optimal range of [AL]T where RAL ≈ 1, although it may decrease at low [AL]T (13, 16). Combining SOV and measured RAL using eq 8 proved successful in modeling results of a CSV study of Cu complexation by Suwanee River fulvic (SRFA) and humic (SRHA) acids (16). In each of these approaches to modeling complexometric titration data, calibration and speciation modeling are handled in separate steps when performing the calculations. That is, one first estimates S in order to transform the observations (I) into the intermediate variables [Mz+] and [∑MLi] and then fits the parameters of the speciation model to these independent variables using one of several graphical and computational methods (17, 23). The use of SSIC and SOV maintains complete separation between the steps by calibrating with a subset of the data where a linear approximation approximates the limit of [M]T . [∑Li]T or R′ . RMLw, respectively, but suffers from bias to the extent that the limit is not reached. Recognizing the potential gain in accuracy over SSIC in relaxing the assumption of linearity lead some analysts to adopt SICI and SRIC, which allow results of the speciation modeling step to affect the calibration to a limited extent. Note that consistency between [∑Li]T, which must be estimated in eq 6, and the sum of the fitted [Li]T is not even required in some forms of calibration. However, in none of these methods can Sˆ and the ligand parameters interact fully. This can result in suboptimal fits to modeled data, to parameter sets that are not optimally coherent, and to underestimates of parameter uncertainties (see Results and Discussion). Objectively, the optimal internal calibration is obtained when one directly fits the actual observed variable (I) by simultaneously estimating S, Ki, and [Li]T using the complete model expressed in terms of the actual independent variables [M]T and [AL]T. While Shuman and Cromer (20) combined the calibration of their ASV technique and the fitting of parameters for a one-ligand Cu speciation model in modeling titrations of humic acid some time ago, their approach has not received much use since. The key reasons may be that (i) no mathematical tools for fitting multi-ligand models in this manner have been developed, (ii) even single-ligand models sometimes do not converge (22), and (iii) directly modeling I does not by itself solve the calibration dilemma. Herein we propose a method that addresses each of these deficiencies. First, even though it has recently been said to not be possible (17, 22), [Mz+] can be solved analytically in two-ligand/one-metal equilibrium systems (Table 1) by expressing the speciation model as a cubic equation, whose solution was first achieved by del Ferro of Bologna ca. 1500 (34). Second, we show that jointly modeling data from multiple titrations at different R′ can minimize the collinearity of S and ligand parameters, thereby permitting simultaneous calibration and speciation modeling with multiple ligands. This approach exploits the constraint of continuity in the underlying function ν relating [∑MLi] to [Mz+] through all titrations of a sample (eq 3). Thus, portions of titration curves at different R′ that overlap in [M]/T must be describable by the same set of speciation model parameters and the S appropriate to each. To our knowledge, the simultaneous calibration of S and estimation of speciation model parameters using titration data from multiple analytical windows has not been attempted previously, although single sets of ligand parameters have been fitted to multiple titrations using fixed values of SOV (16) or SSIC (26, 35).
Modeling and Data Modeling Methods. The basic model equations define I as a function of the independent variables [M]T and [AL]T (or R′) and the parameters given above. These include (i) analytical solutions of speciation models with one (20) or two (Table 1) natural ligand classes at constant R′ (nested VOL. 37, NO. 8, 2003 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
1555
TABLE 1. Analytical Solution for [Mz+] in a System Containing Two Ligands (Li) That Form 1:1 MLi Complexesa [Mz+]3 + a2[Mz+]2 + a1[Mz+] + a0 ) 0
a0 ≡ -
a1 ≡
a2 ≡
[M]T
K1 K2 R +
+ K2[L2]T + K1[L1]T - (K1 + K2)[M]T + R
K1 K2 R + [L1]T + [L2]T - [M]T R+
+
K2 + K1 K1K2
R+ ≡ R′ + RMLw (RMLw assumed constant) [M]T ≡ max ([M]T - [LS]T,0)
P≡
a1 a22 3 9
Q≡
a1a2 a0 a32 - 6 2 27
If (-Q2 > P3) then θ ≡ arcos(Q/x- P3) [M2+] ) 2x- P cos(θ/3) - a2/3 else
U( ≡ [Q ( xQ2 + P3]1/3 [Mz+] ) (U+ + U-) - a2/3 a
Solution to cubic equation adapted from ref 34.
inside a loop that iteratively corrects for the effects of M(AL)x formation on R′), (ii) eq 8, and (iii) eq 4. The SAS code and an approach to testing hypotheses concerning S are described in the Supporting Information. R′ values are computed using the following conditional stability constants (log K1 or β2) for added ligands: CuSA (9.57, 9.77) and Cu(SA)2 (14.57, 15.08) at pH/salinity combinations of 8.0/30.3 and 8.35/38, respectively (36); Cu(bzac)2 (10.7) (23). RCu(oxine)x (17.02) (13) and RCu′ are derived elsewhere (5,24). RSA values used in eq 8 are either directly measured values (Narragansett Bay) or calculated from the following empirical function (Mediterranean Sea): RSA ) {[SA]T(1 + 1.08[SA]T0)}/{[SA]T0(1 + 1.08[SA]T)} for [SA]T in µM (See Supporting Information). RBZAC is assumed to equal unity as in the original analysis. Models are named according to the number of analytical windows fitted and the structure of the speciation model employed. Models of single- and multi-window data sets are denoted “Sx” and “Mx”, respectively, where x ) “G” refers to “general discrete ligand class” and x ) “NC” refers to “natural components” speciation models. Estimates of S made using these models are denoted SSx or SMx. To provide more detail, “G” is usually replaced by the number of ligands employed that have defined Ki with “S” and/or “W” appended if LS or LW are included. Models that employ an Sˆ fixed by a separate calibration step are denoted S(M)x\Sˆ -models. Units of S are always nA (or equivalent measure of peak height)/nM for a constant deposition period. All models were estimated using the Levenburg-Marquadt nonlinear regression algorithm of the SAS/ETS version 8.02 (SAS Institute) “model” procedure, with a set of initial guesses specified for each parameter. To fit SG and MG 1556
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 37, NO. 8, 2003
models, we begin from a speciation model that includes L1 and L2. Next, LS and/or LW are added or L2 is removed as needed to attain convergence and to minimize the relative error in predictions of I, i.e., the RMSE normalized by the error degrees of freedom. Fitting logarithms of parameters (S0, Ki, and [Li]T) facilitates convergence and properly accommodates the large standard errors (SE) of parameter estimates that often occur. In such cases, untransformed parameters often have physically impossible estimated parameter distributions with a significant probability of negative values. Since the error in measurements of I increases with its magnitude (22), residuals were inversesquare-weighted either by (i) I for actual data or (ii) the standard deviation of the error for synthetic data (see Data Sources). The reported uncertainty in each model parameter is the standard error of the estimate (SE) computed by SAS and converted to arithmetic values for [Li]T and S. All parameter estimates are significantly different from an appropriate reference value, e.g., [Li]T > 0.1 nM, Ki > 1 M-1, and RMLw > RM′ (P < 0.05). A ligand class is considered “detected” when its parameter(s) are significant. Modeled [Mz+] are calculated from the best-fit parameters while their uncertainties are derived from SAS Monte Carlo simulations. Using the parameters and covariance matrix from the data fitting step, 104 simulations are run at each [M]T of interest with [AL]T ) 0. The 2.5th and 97.5th percentiles of the distribution of [Mz+] are reported as 95% confidence limits (CL). Data Sources. Original field data from previously published, complexometric Cu(II) titrations were provided by the analysts and are described fully in the Supporting Information. Synthetic CSV/Cu titrations of coastal seawater samples were generated using the “model A” parameters for SRHA (16) in a chemical equilibrium model (KINETEQL) (37) that uses the equilibrium constant approach (38). Titrations of WHAM fulvic acid were generated using model VI (28). Random error was introduced using the normal distribution functions in SAS with standard deviations defined as the maximum of 0.07 nM or 3% of [CuREF] (15, 17). The coefficient defining RCuLw in eq 9 was fitted using error-free data for model B SRHA (16) at [Cu]T ) 2-200 nM. At low metal additions, I is frequently not measurable. Maximum likelihood estimator methods can reduce the bias that either deleting these left-censored data or substituting 0.5× or 1× detection limits introduces into descriptive statistics and regression parameters, but they behave poorly for the small data sets typically used here (39). In preliminary tests, deletion yielded more reasonable parameter estimates than substitution and is adopted here. Developing a means of using the information contained in these data is an important subject for future research.
Results and Discussion Collinearity in Modeling Single Titration Curves. Kozelka and Bruland (24) have described in detail their analysis of an ASV/Cu titration of a coastal seawater sample (Figure 1). From the absence of ASV-labile Cu at ambient [Cu]T and the findings of CSV studies at nearby sites that the ambient Cu was strongly complexed (23), they reasoned that [L1]T ≈ 16 nM. Using SSIC ) 126.8 ( 0.1 and SREF ≈ 175 (est.) as constraints, they obtained SICI ) 142 by inspection. Because the Scatchard linearization of the data yielded a curved plot, they fitted parameters for a two-ligand speciation model: [L2]T ) 20 ( 2 nM, K2′ ) 108.8 ( 0.1; and [L3]T ) 54 ( 4 nM, K3′ ) 107.7 ( 0.05 (Ki′ are stability constants with respect to [Cu′]). In contrast, simultaneously calibrating S and fitting a two-ligand speciation model (model S2), yields parameter estimates of SS2 ) 160 ( 150/80, [L2]T ) 9 ( 14/6 nM, K2′ ) 109.4 ( 0.9, [L3]T ) 90 ( 200/60 nM, and K3′ ) 107.5 ( 0.4. These optimized parameters fit the data (I) with a much
FIGURE 2. S of (A) RMSE, (B) [Li]T, (C) Ki′, and (D) [Cu′] calculated at [Cu]T ) 26 nM for an ASV/Cu titration of coastal seawater (Figure 1). Results of fitting data at fixed S using speciation models with either L2 + L3 (s; dark and medium gray shading ( 1 SE), LS + L3 (- ‚‚ -; ‚‚‚ ( 1 SE), or L3 only (- -; light gray shading ( 1 SE). Single points with error bars are parameter estimates ( 1 SE from models S2, S1S, and S1 plotted at S ) SS2 (157), SS1S (136), and SS1 (91), respectively. For S2 model, log K2′ ) 10 ( 5.6 at S ) 120; at S < 120 model fails to converge. Original RMSE, [Li]T, Ki′, and [Cu′] (9) plotted at S ) SICI (142). In panel D, error bars indicate 95% CL from Monte Carlo simulations for each model; shaded range is 25th-75th percentiles; (-) median. x-y drop line indicates SREF (175). lower RMSE (12%) than the published SICI and ligand parameter values (30%) do using identical model equations. Although the new parameter estimates are not within the published confidence limits, the published estimates are encompassed by the new SE since they are at least 3-50-fold larger. Compensating interactions between model parameters, or collinearity, can amplify the standard errors of parameter estimates and cause nonlinear models to not converge. Collinearity commonly occurs in linear models when independent variables are correlated (21) but can also result from the structure of nonlinear models. That such interactions between S and [Li]T are possible in this context is suggested by the range of x-intercepts of the tangent lines in Figure 1. Collinearity diagnostics (21) for model S2 indicate strong interactionssa condition number of 65 and explained proportions of variation of 0.70 to 0.98sbetween S and [L2]T, K3, and [L3]T in order of increasing severity. Collinearity in this model is also apparent in the correlations of fitted speciation model parameters with S (Figure 2B,C)snegative for the Ki′ and positive for both [Li]T. As a consequence, perturbations in S affect the RMSE very little over the range 120 < S < 200 due to compensatory changes in ligand parameters, making it difficult to calibrate S. Collinearity also amplifies the SE of estimated parameters (Figure 2). Consider that the SE for [L3]T is (17/14 nM with collinearity eliminated by fixing S at SS2 (model S2\160) versus (200/60 nM when S is a free parameter (model S2). The (log) SE of the less strongly collinear parameters are lower by 29-37% at fixed S. These results also illustrate how calibrating first and then fitting speciation model parameters at fixed S can artificially decouple collinear parameters, causing their uncertainty (SE) to be underestimated. In this case, the collinearity and SE are lower for models with fewer parameters (Figure 2). Replacing L2 with LS (one fewer parameter) in model S1S gives SS1S ) 140 ( 40 (RMSE ) 11%) while dropping L2 entirely yields SS1 ) 90 ( 11 (RMSE ) 15%). While still large, the SE of SS1S is much lower because of reduced collinearity, as confirmed by SAS diagnostics. Model S1S’s RMSE is improved over S2 because it has a comparable absolute fit with one less parameter. Note that the remaining parameter estimates are similar to the S2-values (Figure 2B,C) with smaller SE: [L3]T ) 71 (
46/28 nM, K3′ ) 107.6 ( 0.2; [LS]T ) 5.7 ( 1/0.8 nM. Use of model S1 is not recommended because of the degraded fit. The optimized speciation model parameters obtained at various S yield [Cu′] that vary inversely with S (Figure 2D), as suggested by eq 4. The influence of model structure is clear since the overly simple S1 model predicts systematically higher [Cu′] while the two better-fitting models yield similar [Cu′] at the same S and [Cu]T. Finally, [Cu′] computed using best-fit parameters and 95% CL from Monte Carlo simulations are depicted as well. The best combination of low RMSE (Figure 2A) and low uncertainty is obtained with model S1S, confirming that two ligand classes can be detected in these data but not fully parametrized as published. Modeling Multi-Window Data Sets. To demonstrate that jointly modeling data from multiple titrations conducted at different analytical windows resolves the calibration dilemma, we generated synthetic titrations of 3 mg/L SRHA in seawater at pH 8 (S ) 35) using recently reported constants for a twoligand model (16) at R′ ) 102.7, 103.1, 103.5, and 104.1. Although the last 4-5 points in each curve appear linear, the [Cu]T range of 2-200 nM (n ) 18) is only sufficient to half-titrate L2. S2 models of these titrations are affected by collinearity (Figure 3A). Some do not converge, e.g., R′ ) 103.1, while those that do yield inaccurate estimates of parameters, e.g., SS2 ) 0.23 ( 0.64/0.17 at R′ ) 102.7 and SS2 ) 0.65 ( 0.09/0.08 at R′ ) 104.1 rather than the true S of unity (Figure 3). Note that [L2]T is biased low at R′ ) 104.1, a known problem in higher R′ titrations, which results in high calculated [Cu2+] when the stronger ligand classes are titrated (23) (Figure 3D). Both sets of model S2 parameters yield inaccurate and highly uncertain estimates of [Cu2+] as well. In contrast, jointly modeling the data from these two titrations (M2-model) yields SM2 ) 0.99 ( 0.02 and estimates of ligand parameters not significantly different from the true values. The derived [Cu2+] are both accurate and quite precise, e.g., 3.46 ( 0.3 × 10-13 M and 1.11 ( 0.06 × 10-11 M at [Cu]/T ) 10 and 100 nM, respectively (true values are 3.32 × 10-13 and 1.11 × 10-11 M). Moreover, model M2 converges for any pair of these titrations, even including titrations for which S2 does not. As long as the two log R′ values differ by more than 0.4, the SM2 are well-calibrated, i.e., between 0.967 and 1.02 with a mean SE of 0.03. VOL. 37, NO. 8, 2003 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
1557
FIGURE 3. S of (A) RMSE, (B) [Li]T, (C) Ki, and (D) [Cu2+] at [Cu]/T ) 10 and 100 nM. Results of fitting two-ligand speciation model at various fixed S (S2\S or M2\S models) to synthetic CSV/Cu titrations of 3 mg/L Suwanee River humic acid (RAL ) 1) at r′ ) 102.7 (‚‚‚), r′ ) 104.1 (- -), and both r′ jointly (s). Optimal values for S2 model fits to r′ ) 102.7 (4) and 104.1 (]) data shown at SS2 ) 0.23 and 0.65 and for M2 model fit to both jointly (0) and SM2 ) 1.00. Error bars in panel D are Monte Carlo CL. True values (9) at S ) 1. Closer examination of the results of fitting titrations at R′ ) 102.7 and 104.1 illustrates why this occurs (Figure 3). While the RMSE of fits to individual curves are insensitive to S, the RMSE of a joint fit using one set of parameter values (model M2) has a distinct minimum at S ) 0.99 that is nearly equal to the error inherent in the data (3%) (Figure 3). This minimum corresponds to the S at which the [Cu2+] calculated from the two sets of I - [Cu]T data agree most closely at equal [Cu]/T so that a single set of ligand parameters can describe both accurately. In these data (Figure 3D), calculated [Cu2+] at [Cu]/T ) 100 nM agree only at the true S, while at [Cu]/T ) 10 nM, [Cu2+] agree for S > 1. It follows that [L1]T and K1 estimated from the individual titration curves closely agree over a wide range in S, while the optimal values of L2 parameters, particularly [L2]T, agree only at S near unity. Thus, the greater accuracy and precision in the M2 model must result from decoupling S from the most strongly collinear parameter, [L2]T. To apply the overload approach to this sample, we generated additional titration data at R′ ) 105.85 (SSIC ) 0.957 ( 0.008), which overloaded 1 mg/L humic acid in seawater (16), R′ ) 106.18 (SSIC ) 0.970 ( 0.009), and 106.52 (SSIC ) 0.967 ( 0.005). The latter was taken to be SOV with [LS]T ) 0. Although this SOV is quite accurate, proving that it corresponds to overload conditions (assuming that φ > 0) requires obtaining a minimum of two high-R′ titrations that yield relatively little information about natural ligands (other than LS). Fitting model M2\SOV to the R′ ) 102.7 and 104.1 data jointly yields [Cu2+] ) 3.70 ( 0.16 × 10-13 M and 1.16 ( 0.03 × 10-11 at the same [Cu]/T as above. Because of the 3% bias in SOV and the underestimated CL, the true values do not fall within the CL of the model result despite using two more titrations than model M2. This example illustrates the disadvantage of decoupling calibration and speciation modeling. Since the speciation model accounts for the imperceptible curvature in the overload titration, jointly fitting it and any titration with R′ e 104.1 yields a value of SM2 that has a similar precision as SOV and is more accurate (results not shown). The M2-model also yields larger and more accurate estimates of the confidence limits for [Cu2+] than obtained with the separate overload titration, since the uncertainty in SOV is not propagated. Jointly modeling (M2) any two titrations with one R′ g 104.1 and the other R′ e 104.1 also yielded SM2 between 0.99 and 1.01 with correct ligand parameters, except for the pair with the greatest difference in R′ (SM2 ) 1.05 for R′ ) 103.1 and 106.52). These results suggest a more efficient experimental design for speciation studies: conduct necessary titrations 1558
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 37, NO. 8, 2003
FIGURE 4. Parameter estimates for synthetic CSV/Cu titrations of WHAM fulvic acid (5 mg/L). Estimates of K1 (0), [L1]T (9), K2 (O), and [L2]T (b) are for S2 models of single titrations at each r′ and for an M2SW model applied to all titrations simultaneously (ALL). Connected, small filled circles are SS2 and SM2SW (ALL); rCuLw (4) and [LS]T (1) for the M2SW fit. True S is unity. at R′ in the range of interest and a single titration at high R′. Modeled jointly, they should yield more accurate results than possible by estimating SOV independently using one or two fewer titrations. Since the synthetic data used above were generated and fitted using a discrete ligand-class model, one might question whether modeling multiple titrations works for samples with broader distributions of ligand parameters. Synthetic multiwindow titrations of 5 mg/L WHAM fulvic acid, which has a broad distribution of Ki centered at two main values, were generated over a range of [Cu]T from 1 to 200 nM (n ) 22) (28). Fits to individual titrations yield SS2 between 0.73 and 0.87 (Figure 4). As expected (23, 31), the Ki increase in parallel with R′ while [Li]T decrease. Jointly fitting all curves yielded estimates of SM2 with [Cu]T. g Different from S
included. This nearly equals the inherent error of the synthetic data using only 7 parameters (S, 2Ki, 2[Li]T, [LS]T, and RCuLw) rather than the total of 17 employed when each titration is fitted separately. These results also illustrate how one can demonstrate detection of new ligand classes: changes in RMSE will be noticeable as ligands are added, and statistical tests will confirm the significance of the parameters. Multi-Window Titrations of Natural Water Samples. Three published CSV studies of copper complexation in marine systems provide excellent sets of field data for testing the proposed method (Table 2). Rue and Bruland (in ref 23) titrated a filtered and frozen sample collected within a phytoplankton bloom in Narragansett Bay, RI, at five analytical windows using the added ligand salicylaldoxime (SA). Moffett (in ref 23) titrated a fresh portion of the same sample at three windows using benzoyl acetone (bzac) as the added ligand. Finally, Campos and van den Berg’s (13) early study of Cu complexation in NW Mediterranean Sea water employed one oxine and three SA windows. Originally, each titration in these studies was analyzed individually, using some form of internal calibration and one- or two-ligand speciation models. Since consistency between titrations was not a constraint, the resulting model curves are discontinuous on [Cu]/T plots (Figures 5A-7A). Furthermore, for several of the curves the high-[Cu]/T points are aligned vertically, as is characteristic of the “linear” points used in simple internal calibrations (17). Finally, the overall RMSE of I values
predicted using the original model parameters are quite larges50, 24, and 43%, respectively (Table 2)sdue in part to over-parametrization. For example, the overall RMSE for the Narragansett SA study is 50% since 22 parameterss5 × (SSIC, [∑Li]TSIC, [L1]T, K1) + 1 × ([L2]T, K2)swere fitted to 29 data points, leaving only 7 error degrees of freedom. Simultaneous calibration and speciation modeling using generic discreet ligand-classes (MG models) yield improved fits of these multi-window data sets. The RMSE decrease to 12-19% (Table 2) and the [Cu2+] - [Cu]/T plots (Figures 5B7B) are visually superior to those from the original analyses. SIC As expected, the calibrated SMG 0 both exceed S0 and remain below SREF (where known). Although overload conditions 0 SIC were not attained in any study, the ratios of SMG do 0 :S0 decrease from 2.0 to 1.1 with increasing R′ as predicted from eq 7. All of these considerations confirm that the SMG are an improvement over the original Sˆ . Note that surfactants clearly did affect the SA titrations, since φ is 0.16 ( 0.06 and ∼0.3 in the two SA data sets. The assumption that φ is independent of [AL]T (eq 8) is validated if not quite proven by the quality of these fits. The MG models yield ligand parameters that for the most part seem similar to the original estimates (Table 2). All three data sets have concentrations of strong ligands (Ki g 1012.6) that exceed ambient [Cu]T. However, the ambient [Cu2+] derived here are 3-6-fold lower than the original estimates. VOL. 37, NO. 8, 2003 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
1559
FIGURE 5. CSV/Cu/SA titrations of Narragansett Bay water (23). [SA]T ) 55.5 (9), 27.5 (3), 11 (1), 2.75 (O), and 1.1 µM (b). (A) Original analysis. (B) M2 model (L1 + L2) fit to all five titrations jointly (s) with 95% CL (- -). Model curve from Figure 6B (gray). (C) Natural component model curves for HA only (- -), HA + LS (- ‚‚ -), and HA + LS + L1 (s). Model curve from Figure 6C (gray). Plots that use [Cu2+] - [Cu]/T coordinates actually introduce considerable uncertainty in both x and y coordinates of data when S is uncertain. Since this effect is systematic across all data points and is accounted for in the models, we have chosen to depict the resulting uncertainty in model results rather than use bi-directional error bars, which suggest independent errors for each point.
FIGURE 6. CSV/Cu/bzac titrations of Narragansett Bay water (23). [bzac]T ) 500 (1), 250 (O), and 100 µM (b). (A) Original analysis. (B) M1SW model (L1 + LS + LW) fit to all three titrations jointly. Lines as in Figure 5B. Model curve from Figure 5B (gray). (C) Natural component model curves as in Figure 5C, except without LS. Model curve from Figure 5C (gray).
FIGURE 7. CSV/Cu/SA/oxine titrations of Mediterranean Seawater (23). [SA]T ) 10 (1), 2 (O), and 1 µM (b), [oxine]T ) 1 µM (3). (A) Original analysis. (B) L1 + L2 model fit to all three SA titrations jointly. Lines as in Figure 5B. Oxine data calculated using optimized S. (C) Natural component model curves as in Figure 5C. The two lowest [Cu]T points in the 1 and 2 µM [SA]T titrations (panel A) have reported I < 0.1 nA (detection limit) and are excluded from panels B and C. 1560
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 37, NO. 8, 2003
SIC This is caused both by the fact that SMG 0 > S0 and by our use of the highest R′ titration data in calculating [Cu2+], which permits ligand-classes with higher Ki to be characterized (31). Note that the original analyses of some of these data yielded [L1]T at or below ambient [Cu]T (21) and could therefore not be used. The [Cu2+] derived here are lower than the 10-12.210-12.8 M recently reported for pristine coastal U. S. surface waters (40) and are comparable to or lower than the 10-13.1 M reported in the central northeast Pacific Ocean (2). The accuracy of the results obtained with the MG models is supported by the good agreement that exists between the modeled [Cu2+]T - [Cu]/T curves for the SA and bzac data sets from the Narragansett Bay sample. Above [Cu]/T ) 15 nM (Figures 5B and 6B), the model curve for the bzac data mostly lies within the 95% CL of the SA model and vice versa. Potential causes of the differences that do exist include the difference in pH of 0.1 unit, the effects of freezing the water sample used in the SA study on colloidal ligands, and possible discrepancies in the Cu(AL)x stability constants. The most important difference may lie in the sequence of reagent additions, which caused [Cu2+] to approach equilibrium from below (high [Cu]/T) in the SA titrations and from above (low [Cu]/T) in the bzac titrations. Note that (i) unlike the other data sets neither of the two models of the SA data (Figure 5B,C) match any of the individual curves in their entirety very well, although they are individually very coherent and (ii) the highest R′ curve is quite poorly described. It may be worth investigating whether portions of this data set were not at equilibrium. The Mediterranean Sea data set (13) also contains one titration conducted using a second AL, oxine (Figure 7A). While we cannot estimate an independent set of parameters from a single curve, the SM2(oxine) estimated by finding the value that permits them to agree most closely with the SA model (Figure 7B) is 2.2 times the original SSIC(oxine) (RMSE ) 20%, r 2 ) 0.80 for the oxine data). The oxine and SA titration data can also be jointly fitted using model M2, yielding a total RMSE of 14% without changing model parameters significantly. An nearly identical increase in SMNC(oxine) is obtained with the MNC model below (RMSE ) 17%, r 2 ) 0.86). Natural Components. Empirically describing the function ν is essential for estimating ambient [Cu2+] but does not serve as a basis for comparing the contributions of different natural biogeochemical components to metal complexation in natural waters (19). For example, given the strong Cu complexation observed in the above data sets, it is of great interest to determine whether high-affinity Cu chelators, such as those produced by cyanobacteria or sulfide clusters or background humic and fulvic acids in a sample, are the major sources of ligands (8, 16, 17). To provide a better means of testing such hypotheses than reporting [L1]T and [L2]T, we suggest an alternate parametrization of the discrete ligandclass model. This “natural components model” expresses the concentration of ligands in each class i as the sum of contributions from humic acid present at a concentration [HA]NC and any remaining ligands [Li]NC T :
[Li]T ) [HA]NCni + [Li]NC T
(9)
Following others (16, 27), we represent complexation by humic acids using fixed stability constants (K1 ) 1013; K2 ) 1011.5) and constant quantities of humic ligands in each class (n1 ) 3.4, n2 ) 139 nmol/mg HA). The weak complexation parameter is assumed to be proportional to the concentration of humic acidsRCuLw ) 102.942[HA]NCswhile [LS]NC T is independent of [HA]NC. To apply this model, the concentrations NC of the natural componentss[HA]NC, [L1]NC T , and [LS]T sare estimated rather than the ligand parameters of the MG model.
The natural component models (MNC) fit the field data sets nearly as well as the general ligand-class (MG) models, except in the case of the Narragansett Bay SA data (Table 2; Figures 5C-7C), and the fitted SMNC are not significantly different from the SMG. Although the calculated [Cu2+] in all cases are slightly greater than the MG model values, they have overlapping confidence limits. These analyses of complexometric titrations of water collected during the spring phytoplankton bloom in Narragansett Bay and from the Mediterranean Sea suggest that high-affinity chelators do regulate [Cu2+] in these marine systems (10). The concentraNC tions of these ligands ([L1]NC T + [LS]T ) are higher than the ambient [Cu]T by at least 10 nM in all cases. In each case, ambient [Cu]T > [LS]T but a large excess of L1 depresses [Cu2+]. Interestingly, while the MNC model does not fit the Narragansett Bay SA data well, the MNC models of the SA and bzac data sets agree even more closely at [Cu]/T > 15 nM than the MG models do (Figures 5C and 6C) and the modeled concentrations of high-affinity ligands are similar; [L1]NC T + [LS]NC T ) 21 ( 3 nM for the SA titrations and 22 ( 1 nM for the bzac data. The higher R′ of the SA titrations permit the distinction between L1 and LS classes to be made but do not quantify complexation by humic acid as precisely as the lower R′ bzac data. A recent CSV analysis of filtered water (salinity ) 28) gathered during wintertime from another site in Narragansett Bay (41) found Cu-binding ligands equivalent to 0.27 mg humic equiv/L and 4-8 nM inert LS but no L1 class ligands. Of course, this natural components approach requires a valid model of complexation by humic and fulvic substances. Kogut and Voelker’s SRHA constants (16) suit our purposes very well since they were (i) obtained using CSV in seawater and (ii) cover the appropriate range of [Cu2+]. The differences between the composition of the medium in which the experimental constants were obtained (pH 8.2; S ) 35) and the samples from Narragansett Bay (pH 8.0, S ) 30.3) and the Mediterranean Sea (pH 8.35, S ) 38) should exert only modest effects on the complexation of Cu by humic and fulvic acid according to WHAM calculations (28). At any given [Cu2+], humic-associated CuLi should differ from the reference conditions by a relatively constant amount. This bias of -13% at pH 8 (S ) 30.3) and +15% at pH 8.35 (S ) 38) acts to narrow the differences between the [HA]NC of the samples estimated here (Table 2). Thus, we interpret the [HA]NC parameter as a measure of the intensity of humic complexation of Cu by an equivalent concentration of SRHA rather than an actual concentration of humic acid in a sample. Ultimately validated models for humic/fulvic acids that account for changes in medium compositionspH, salinity, and competing metalssshould be incorporated into the modeling of natural components. Other useful adaptations of this approach will likely involve developing correlations between field measurements of the physicochemical properties of NOM (42) and quantities of other natural ligands (4) and speciation model parameters derived from modeling complexometric titrations of natural water samples.
Acknowledgments We are grateful to L. Campos, P. Kozelka, J. Moffett, and C. van den Berg for providing data and to G. Gertner and S. Wente for helpful advice on statistical issues. Suggestions by C. van den Berg, M. Kogut, C. W. Boast, and anonymous reviewers helped improve the manuscript. This material is based upon work supported by the Cooperative State Research, Education and Extension Service, U.S. Department of Agriculture under Project 399 (R.J.M.H.) and by the National Science Foundation under Grants OCE-0137085 and CEBIC-34460011 (K.W.B. and E.L.R.). VOL. 37, NO. 8, 2003 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
1561
Supporting Information Available Raw data, RSA values, and SAS code including 5 tables and 1 figure. This material is available free of charge via the Internet at http://pubs.acs.org.
Literature Cited (1) Sunda, W. G.; Guillard, R. R. L. J. Mar. Res. 1976, 34, 511-529. (2) Bruland, K. W.; Donat, J. R.; Hutchins, D. A. Limnol. Oceanogr. 1991, 36, 1555-1577. (3) Morel, F. M. M.; Hering, J. G. Principles and Applications of Aquatic Chemistry; John Wiley & Sons: New York, 1993. (4) Campbell, P. G.; Errecalde, O.; Fortin, C.; Hiriart-Baer, W.; Vigneault, B. Comp. Biochem. Physiol., Part C: Toxicol. Pharmacol. 2002, 133, 189-206. (5) Byrne, R. H.; Kump, L. R.; Cantrell, K. J. Mar. Chem. 1988, 25, 163-181. (6) Al-Farawati, R.; van den Berg, C. M. G. Mar. Chem. 1999, 63, 331-352. (7) Breault, R. F.; Colman, J. A.; Aiken, G. R.; McKnight, D. Environ. Sci. Technol. 1996, 30, 3477-3486. (8) Xue, H.; Sigg, L. Aquat. Geochem. 1999, 5, 313-335. (9) Wells, M. L.; Kozelka, P. B.; Bruland, K. W. Mar. Chem. 1998, 62, 203-217. (10) Moffett, J. W.; Brand, L. E. Limnol. Oceanogr. 1996, 41, 388393. (11) Martell, A. E.; Smith, R. M.; Motekaitis, R. J. NIST Standard Reference Database 46 Version 6.0; NIST Standard Reference Data: Gaithersburg, MD, 2001. (12) Bruland, K. W.; Franks, R. P.; Knauer, G. A.; Martin, J. M. Anal. Chim. Acta 1979, 105, 223-245. (13) Campos, M. L. A. M.; van den Berg, C. M. G. Anal. Chim. Acta 1994, 284, 481-496. (14) Coale, K. H.; Bruland, K. W. Limnol. Oceanogr. 1988, 33, 10841101. (15) Miller, L. A.; Bruland, K. W. Anal. Chim. Acta 1997, 343, 161181. (16) Kogut, M. B.; Voelker, B. M. Environ. Sci. Technol. 2001, 35, 1149-1156. (17) Voelker, B. M.; Kogut, M. B. Mar. Chem. 2001, 74, 303-318. (18) Turoczy, N. J.; Sherwood, J. E. Anal. Chim. Acta 1997, 354, 1521. (19) Town, R. M.; Filella, M. Limnol. Oceanogr. 2000, 45, 13411357. (20) Shuman, M. S.; Cromer, J. L. Environ. Sci. Technol. 1979, 13, 543-545. (21) Belsley, D. A. Conditioning Diagnostics: Collinearity and Weak Data in Regression; John Wiley & Sons: New York, 1991.
1562
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 37, NO. 8, 2003
(22) van den Berg, C. M. G. Mar. Chem. 1984, 15, 1-18. (23) Bruland, K. W.; Rue, E. L.; Donat, J. R.; Skrabal, S. A.; Moffett, J. W. Anal. Chim. Acta 2000, 405, 99-113. (24) Kozelka, P. B.; Bruland, K. W. Mar. Chem. 1998, 60, 267-282. (25) Xue, H.; Sigg, L. In Environmental Electrochemistry; Taillefert, M., Rozan, T. F., Eds.; ACS Smposium Series 811; American Chemical Society: Washington, DC, 1999; pp 313-335. (26) Xue, H.; Sunda, W. G. Environ. Sci. Technol. 1997, 31, 19021909. (27) Westall, J. C.; Jones, J. D.; Turner, G. D.; Zachara, J. M. Environ. Sci. Technol. 1995, 29, 951-959. (28) Tipping, E. Aquatic Geochem. 1998, 4, 3-48. (29) Kinniburgh, D. G.; Milne, C. J.; Benedetti, M. F.; Pinheiro, J. P.; Filius, J.; Koopal, L. K.; Van Riemsdijk, W. H. Environ. Sci. Technol. 1996, 30, 1687-1698. (30) Fish, W.; Dzombak, D. A.; Morel, F. M. M. Environ. Sci. Technol. 1986, 20, 676-683. (31) van den Berg, C. M. G.; Nimmo, M.; Daly, P.; Turner, D. R. Anal. Chim. Acta 1990, 232, 149-159. (32) Leal, M. F. C.; van den Berg, C. M. G. Aquat. Geochem. 1998, 4, 49-75. (33) van den Berg, C. M. G. University of Liverpool, personal communication, 2002. (34) Gellert, W.; Ku ¨ stner, H.; Hellwich, M.; Ka¨stner, H. The VNR Concise Encyclopedia of Mathematics; Van Nostrand Reinhold: New York, 1975. (35) Moffett, J. W.; Zika, R. G. Mar. Chem. 1987, 21, 301-313. (36) Kogut, M. B. Ph.D. Thesis, Massachusetts Institute of Technology, 2002. (37) Giambalvo, E. M.S. Thesis, University of California at Santa Cruz, 1997. (38) Westall, J. C. MICROQL. I. A Chemical Equilibrium Program in BASIC: Department of Chemistry, Oregon State University: 1986. (39) Clarke, J. U. Environ. Sci. Technol. 1998, 32, 177-183. (40) Twiss, M. R.; Moffett, J. W. Environ. Sci. Technol. 2002, 36, 10611068. (41) Kogut, M. B.; Voelker, B. M. Environ. Sci. Technol. 2003, 37, 509-518. (42) Ross, A. R. S.; Ikonomou, M. G.; Orians, K. J. Anal. Chim. Acta 2000, 411, 91-102.
Received for review April 26, 2002. Revised manuscript received December 6, 2002. Accepted December 18, 2002. ES025751A
