Surface complexation modeling in Variable charge SoilS : prediction of cadmium adSorption

intrinsic equilibrium constants for 22 representative brazilian oxisols were estimated from a cadmium adsorption experiment. Equilibrium constants were fitted to two surface complexation models: diffuse layer and constant capacitance. intrinsic equilibrium constants were optimized by fiteQl and by hand calculation using Visual minteQ in sweep mode, and excel spreadsheets. data from both models were incorporated into Visual minteQ. constants estimated by fiteQl and incorporated in Visual minteQ software failed to predict observed data accurately. however, fiteQl raw output data rendered good results when predicted values were directly compared with observed values, instead of incorporating the estimated constants into Visual minteQ. intrinsic equilibrium constants optimized by hand calculation and incorporated in Visual minteQ reliably predicted cd adsorption reactions on soil surfaces under changing environmental conditions.


introduction
Anthropogenic activities, such as mining or industrial activities, agricultural chemical applications -mainly phosphate and zinc micronutrients (Guilherme et al., 2014), and organic and inorganic residue application to soils may cause cadmium (Cd) concentration to build up over time.Cadmium is bound to the solid phase of soils as a result of surface precipitation and adsorption processes.Its mobility and destination in the soil environment are directly related to processes that are highly pH dependent (Lee et al., 1996).
Soil contamination by Cd is of great concern in terms of its entry into the food chain for it can be taken up by food plants in large amounts relative to its concentration in the soil.Furthermore, Cd accumulates over a lifetime in the body (ATSDR, 2012) and has adverse health effects, mainly in the form of renal dysfunction (Hooda, 2010).Environmental and human health risk assessment of metals depends to a great extent on modeling the destination and mobility of metals based on soil-liquid partitioning coefficients (Sauvé et al., 2000) that can be estimated locally.
Cadmium adsorption is mainly controlled by soil organic matter (Lee et al., 1996;Shi et al., 2007), pH (Sauvé et al., 2000), cation exchange capacity, clay content or type, surface area (Appel and Ma, 2002), surface charge (He et al., 2005), or a combination of those soil components.Different components in soils may contribute to Cd adsorption to different extents (Shi et al., 2007).Indeed, controlling the factors of Cd adsorption in highly weathered soils, such as Oxisols, are not clear (Alleoni and Camargo, 1994).Correlation coefficients between Cd adsorption and many soil properties for 17 Oxisols were inconclusive at pH 4.5, but highly significant at pH values of 5.5 and 6.5.Properties such as surface area, CEC, clay, kaolinite, and hematite contents were highly correlated (Pierangeli et al., 2003), emphasizing the importance of pH, surface area, and charge for Cd adsorption in these soils.Understanding mechanisms of metal adsorption in soils is important as these reactions control the strength of the metal-soil surface interactions (He et al., 2005).
Oxisols are characterized by B horizons with 1:1 clay minerals and variable amounts of Fe and Al hydrous oxides (Reatto et al., 2009a), exhibiting generally low surface charge density with predominant pH-dependent surface charge.Charge develops from protonation and deprotonation reactions of hydroxyl groups protruding from the surfaces of the minerals, and by deprotonation of acidic functional groups exposed on the surfaces of solid organic matter (Charlet and Sposito, 1987).
Little research has been done directly comparing surface charge with heavy metal sorption in variable charge systems (Appel and Ma, 2002).In this regard, surface complexation models have distinct advantages over empirical models (which provide soil-liquid partitioning coefficients for further modeling) because surface complexation models can be extrapolated to systems of different ionic strengths, pH 's, and component compositions (Koretsky, 2000;Bethke, 2008).However, application of the surface complexation modeling approach to natural systems, such as soils, is considerably more difficult, and rare, than application to simple mineral-water systems (Davis and Kent, 1990;Davis, 2001).Theory needs to describe hydrolysis and the mineral surfaces, account for electrical charge, and provide for mass balance on the sorbing sites.In addition, an internally consistent and sufficiently broad database of sorption reactions should accompany the theory (Bethke, 2008).Estimated intrinsic equilibrium constants must ultimately describe reliable results within chemical speciation software packages, such as Visual MINTEQ (Gustafsson, 2014), allowing end users to perform predictions and apply results in further technologies.
Intrinsic equilibrium constants estimated from Cd adsorption in soils, optimized by FITEQL and by hand calculation may be incorporated into Visual MINTEQ.Results may be validated as comparing observed results of Cd adsorption in four Oxisols with those predicted by Visual MINTEQ.
The overall objectives of this study were to estimate intrinsic equilibrium constants of Cd

Soil samples and properties
Soil samples from the 0.00-0.20 m depth from 20 sites were collected and characterized for mineralogy, texture, organic carbon, and cation exchange capacity at pH 7.0 by Rein (2008); additional soils from horizons A and Bw from three sites were collected and characterized by Reatto et al. (2007Reatto et al. ( , 2008Reatto et al. ( , 2009a,b),b).The characterization data, displayed in tables 1 and 2, were used in the present study.Site location was based on survey reports to provide a representative sample of pristine Cerrado (tropical savanna) Oxisols from Brazil.Consideration was given to geographic distribution, variations in texture, mineralogical composition of the clay fraction, and parent material.All sampled sites were under native vegetation, and had never been cropped or amended.Determination of surface area was performed as described by Cerato and Lutenegger (2002), with modifications.Soil was passed through a 2 mm sieve, oven dried for 12 h at 120 ºC to remove water, and then placed (1 g) in an aluminum tare.Three mL of ethylene glycol monoethyl ether (EGME) were mixed with the soil.Tares were placed in a desiccator and a 635 mm Hg vacuum was applied.One hundred and 10 grams of oven-dried CaCl 2 was used as a desiccant.After 18 and 24 h, the desiccator was opened and the tares were weighed.If the mass of the sample varied more than 0.001 g from the first time to the second, the sample was placed once more in the desiccator and evacuated again for 4 h more.The procedure was repeated successively until the tares reached constant weight.Soil surface area (a) was calculated as A Wa = Ws 2.86 × 10 -4 , where Wa is the weight of EGME retained by the sample in g; Ws is the oven dry weight of soil; and 0.000286 is the weight of EGME required to form a monomolecular layer on a square meter of surface (m 2 g -1 ).Reference clays from the Source Clay Mineral Repository at the University of Missouri at Columbia, Missouri, USA, were used to validate the method.

adsoption isotherms
Cadmium adsorption envelopes (amount of Cd adsorbed as a function of solution pH at a fixed total Cd concentration) were determined for all soil samples.One g of soil was placed in polypropylene centrifuge tubes and brought to equilibrium with 15 mL of a 0.1 mol L -1 NaCl.Tubes were agitated at 200 rpm for 12 h in a reciprocating shaker.Acid or base were added, and the operation was repeated several times, up to 72 h total time, to adjust initial suspension pH values to 4.5, 5.0, 5.5, 6.0, and 6.5.Each treatment had three replications.Five mL of 0.36 mmol L -1 Cd solution in 0.1 mol L -1 NaCl was mixed with the soil suspensions.Final Cd concentration in the soil suspension was 0.09 mmol L -1 .Soil suspensions were shaken at 200 rpm for 24 h, and the pH of suspensions at the end of this procedure was recorded.Soil suspensions were centrifuged at 8000 rpm for 10 min, and the supernatant was collected for Cd analysis by inductive coupled plasma with mass spectrometry (ICP-MS).
A review of the theory and assumptions of the constant capacitance model, as well as of the diffuse layer model was presented in Goldberg (1992).The soil will be considered as an integrated whole with an average surface functional group (expressed as SOH).The surface protonation and deprotonation of the soil are expressed as equations 1 and 2: The adsorption of Cd on the soil surface is expressed by the following surface complexation reaction (Equation 3).

SOH Cd
SOCd H Eq. 3 Other reactions with Cd and Cl -were also included in the modeling, e.g., CdCl + and CdCl 2 species played an important role in solution equilibria (equations summarized in table 3), even though CdCl + was not considered as a surface-bound species.
Intrinsic equilibrium constant expressions for the surface complexation reactions are: Eq. 6 where exp(± FΨ/RT) is derived from the Boltzmann equation and used to adjust for the electrostatic properties of the charged surfaces; F is the Faraday constant (C mol c -1 ); Ψ is the surface potential (V); R is the molar gas constant (J mol -1 K -1 ); T is the absolute temperature (K); and square brackets indicate concentrations (mol L -1 ).
Optimal best-fit values for the surface site acidity constants, log K a1 int and log K a2 int , as well as surface site density for Oxisols from the Brazilian Cerrado, were calculated from a previous study (Table 1).These constants were fixed when the metal binding constants were calculated from the batch adsorption data.
The Davies equation (Davies, 1938)and a revision of the dissociation constants of some sulphates Journal of the Chemical Society (Resumed and the specific interaction theory (Sukhno and Buzko, 2004), both used for activity correction in Visual MINTEQ, were applied to adjust activity coefficients in extrapolation of equilibrium constants in i = 0.1 mol L -1 (Table 3).Equilibrium constants of dissolved species extrapolated to i = 0.1 mol L -1 by the Davies equation (Equation 7) were used to optimize log K SoCd int by FITEQL 4.0 (Herbelin and Westall, 1999), and were used in further modeling through Visual MINTEQ.There is no generally accepted convention for treating activity coefficients of surface species (Richter et al., 2005).In general, there is little or no change in the pH dependence of adsorption with ionic strength (between 0.001-1.0mol L -1 ) for specifically adsorbed ions like Cu 2+ , Pb 2+ , Ni 2+ , and Cd 2+ (Hayes and Leckie, 1987); therefore, adopted protonation/dissociation constants were obtained from a previous study in 0.1 mol L -1 NaCl.The Davies equation is: where the subscript i refers to each of the reactants and products in the reaction; z i is the ionic charge of each reactant or product; and i is the ionic strength reported for the experimental data.
Once computed, the activity coefficients were used in the following relationship to correct the equilibrium constants to I = 0.1 mol L -1 : Eq. 8 where v represents the reactant or product stoichiometric coefficient.
The dominant Cd solution species in groundwater at pH values greater than 8.2 is CdCO 3 (aq).Precipitation with carbonate is increasingly important in systems with a pH greater than 8 (USEPA, 1999).However, Brazilian Oxisols are rather acidic soils; therefore, the effect of partial CO 2 pressure, and its related compounds, in simulated suspensions was not included in Visual MINTEQ or FITEQL calculations.
The FITEQL program was used to optimize surface complexation modeling of the experimental Cd adsorption.The program was set to allow ionic strength and activity corrections.FITEQL uses the Davies equation to correct ionic activity (Equation 7).Obtained equilibrium constants were averaged using the weighting method of Dzombak and Morel (1990), in which the weighting factor w i is defined as: ) Eq. 9 where (σ log K )i is the standard deviation of log K calculated by FITEQL for the i th data set.The best estimate for log K is then calculated as: Eq. 10 Another method to optimize log K SoCd int values was employed.Using observed data as a guide parameter, simulations with Visual MINTEQ were used to optimize by trial and error (systematic variations of assumed values of log K SoCd int ) using the method of least squares regression (with a log K SoCd int precision of 0.01).This method seeks to minimize the sum of the squared errors (SSE) between observed and calculated values of the dependent variable, in this case the sorbed Cd concentration, S (Bolster and Hornberger, 2007): Eq. 11 where SSE is the objective function to be minimized; N is the number of observations; Si is the i th measured value of the dependent variable; and Si is the i th model-predicted value of the dependent variable.
The constant capacitance model is very insensitive to values of capacity density (C 1 ) (Goldberg, 1995).The choice of the C 1 value is arbitrary, and it is recommend to use the best fit values (~1.0 F m -2 ) (Hayes et al., 1991).We chose to use the C 1 value of 1.06 F m -2 (derived from Al oxides) (Westall and Hohl, 1980), which is usually used for soil modeling (Goldberg et al., 2000).

Validation
The validation of the surface complexation models obtained with adsorption of Cd in 16 soil samples collected by Rein (2008) and six by Reatto et al. (2008) was carried out using four soil samples from different localities (Sete Lagoas, São Gotardo, Uberaba, and Itumbiara), collected in a transect by Rein (2008).These four soil samples were not used during model prediction.The set up for the validation experiments was the same as that used for prediction.All model parameters have been maintained unvaried; therefore, the results presented are a test of the ability of the model to simulate Cd adsorption in these four soils.Observed vs. predicted plots of Cd adsorption values from this series of soils were compared using the root mean square error (RMSE) and a dimensionless statistic which directly relates model predictions to observed data such as modeling efficiency (EF) (Mayer and Butler, 1993): Eq. 12 where y i represents observed values; ŷ i , simulated values; and n, the number of pairs.

EF ŷ y ŷ y
Eq. 13 where y ¯ represents the observed mean.
Intrinsic equilibrium constants estimated by FITEQL were 0.29 and 0.38 units smaller for DLM and CCM, respectively, than those calculated by hand (Table 4).Confidence intervals were also narrower for the hand calculation group than those calculated by FITEQL.Very similar statistical results (standard deviation and confidence intervals) were found for DLM and CCM.Intrinsic equilibrium constants optimized by FITEQL and hand calculation differed significantly (Table 4).A possible explanation for the difference between log K SoCd int results from FITEQL and from hand calculation using Visual MINTEQ may be related to FITEQL numerical instability or convergence problems (Villegas-Jiménez and Mucci, 2009).Speciation constants for major reactions used in FITEQL exactly matched those in Visual MINTEQ (Table 3).Output data from FITEQL, however, was able to predict the observed data with precision (Figure 1), and preliminary tests optimizing constants without taking into account background ionic strength (0.1 mol L -1 NaCl; in both software programs, FITEQL and Visual MINTEQ) were satisfactory, and predicted the observed results perfectly, but they did not correspond to reality.When optimized intrinsic equilibrium constants estimated by FITEQL were added to the Visual MINTEQ database including background ionic strength data, the program failed to accurately predict observed results of Cd adsorption (Figures 2 and 3).Some authors use site density to ultimately optimize curve fitting (Davis and .997 -13.782 -13.796 (1) Suspension concentration = 50 g L -1 ; site concentration (Nt; mol kg -1 ) may be converted to site density (Ns; sites nm -2 ) by the following expression: Nt = (SA × Ns × CS × 10 18 )/ NA, where SA is the surface area (m 2 g -1 ); CS is g of soil L -1 ; and NA is Avogadro's number; (2) 298.15K; (3) Data from Martell and Smith (2003); (4) i: ionic strength in mol L -1 ; ( 5) Davies ( 1938); (6)  Sukhno and Buzko (2004).Surface complexation modeling in Variable charge SoilS: prediction of cadmium... Kent, 1990) because as site density increases, the adsorption edge shifts to lower pH values, whereas with decreasing site density, the adsorption edge shifts to higher pH values, allowing the modelist to adjust curve fitting "manually".Using changes in site density to adjust curve fitting was not our intent in the present study.Rather, we searched for a hand calculation method that could provide more accurate results when using Visual MINTEQ, and that could be used as a reference, since as it was optimized accessing by Visual MINTEQ data directly (Table 4).
Intrinsic equilibrium constants optimized by hand calculation, considering a 0.1 mol L -1 NaCl background, was incorporated in the Visual MINTEQ database.The results showed good precision in reproducing the observed values (Figures 2 and 3).Using this approach, the diffuse layer and constant capacitance models produced very similar results.In fact, differences between models are insignificant (Davis and Kent, 1990) when models have matching pH and ionic strength.A sample from an Oxisol from Buritis, MG, based on estimated parameters from hand calculation, exhibited values of 0.14 and 0.15 for RSME for Cd adsorption, and 0.97 and 0.96 for EF in regard to DLM and CCM, respectively, showing very similar prediction values between models (Figures 2 and  3).When parameters estimated from FITEQL were added to Visual MINTEQ, RSME was high, and EF was low for both DLM and CCM (Figures 2  and 3).The latter statistical results clearly show disagreement among FITEQL optimization of log K SoCd int added into Visual MINTEQ and observed data (Figures 2 and 3).

Validation
Cadmium adsorption in four Oxisol samples collected in Sete Lagoas, MG; São Gotardo, MG; Uberaba, MG; and Itumbiara, GO, was modeled using parameters estimated by hand calculation in Visual MINTEQ (Table 4, Figure 4).Parameters derived from adsorption obtained at several fixed pH values were incorporated in the DLM to reproduce pH edge curves.Only DLM using optimization by hand calculation was shown in figure 4 since DLM and CCM showed equivalent results (Figure 1).Another reason for showing only DLM using optimization by hand calculation in figure 4  CI of mean (4) 0.18 0.17 (1) Suspension concentration = 50 g L -1 ; Capacitance of inner Helmholtz layer = 1.06 F m -2 ; Site density = 0.4295 and 0.4444 sites nm -2 for DLM and CCM, respectively; (2) WSOS is the weighted sum of squares of the residuals and DF is the degrees of freedom; (3) SD: standard deviation from log K values; (4) CI of mean = 95 % confidence interval; (5) SSE: sum of the squared errors.4).Although it is interesting to know the likely range of a true value using log K SoCd int , it is important to show that each value of log K SoCd int is related to its soil surface area.Therefore, the confidence interval established using only one value of soil surface area (from the sample under analysis) shows narrowed confidence intervals.Taking soil surface area into account, calculated confidence limits would probably show larger intervals.Therefore, another model of goodness of fit statistics was assessed, for purposes of validation, using RMSE and EF (Equations 12 and 13).In the end, cadmium adsorption was modeled satisfactorily in all four soils, according to validation statistics (Figure 5).
-    With the use of log K SoCd int estimated by hand calculation, predicted Cd adsorption vs. observed values rendered lower RSME and higher EF than the log K SoCd int estimated by FITEQL (when running these constants under Visual MINTEQ).
figure 5. observed vs. predicted cadmium adsorption as a function of ph in four oxisols from Sete lagoas, mg; São gotardo, mg; uberaba, mg; and itumbiara, go.intrinsic equilibrium constants optimized by hand calculation using Visual minteQ were used as input values in Visual minteQ for cd prediction using the diffuse layer model; suspension concentration = 50 g l -1 ; 0.1 mol l -1 nacl; log K SOCd int = -1.95.
log K SoCd int , as well as all other parameters involved, such as surface acidity constants and site concentration, are ready to be incorporated in chemical equilibrium software.Diffuse layer and constant capacitance model parameters, when incorporated in chemical equilibrium software, allow prediction of Cd adsorption in Oxisols with high efficiency and accuracy.Although FITEQL output rendered significant predicted vs. observed values of Cd adsorption, when the optimized log K SoCd int was incorporated into Visual MINTEQ, predicted Cd adsorption vs. observed values rendered high RSME and low EF, showing discrepancies between programs.

table 4 . average intrinsic equilibrium constants for cd adsorption onto 20 oxisols from the brazilian cerrado using the diffuse layer model (dlm) and the constant capacitance model (ccm), estimated by fiteQl 4.0 and by hand calculation using Visual minteQ (1)
was that FITEQL optimized