Assessment of Modeling Uncertainties Using a Multistart Optimization

Article

Subscriber access provided by FONDREN LIBRARY, RICE UNIVERSITY

Assessment of modeling uncertainties using a multi-start optimization tool for surface complexation equilibrium parameters (MUSE) Nefeli Bompoti, Maria Chrysochoou, and Michael L. Machesky ACS Earth Space Chem., Just Accepted Manuscript • DOI: 10.1021/ acsearthspacechem.8b00125 • Publication Date (Web): 27 Dec 2018 Downloaded from http://pubs.acs.org on January 2, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Earth and Space Chemistry

Assessment of modeling uncertainties using a multi-start optimization tool for surface complexation equilibrium parameters (MUSE) Nefeli Maria Bompoti*,†, Maria Chrysochoou† and Michael L Machesky‡

†

Department of Civil and Environmental Engineering, University of Connecticut, Storrs, CT USA ‡

University of Illinois, Illinois State Water Survey, Champaign, IL USA

ABSTRACT The MUSE algorithm has been developed to optimize the fitting of thermodynamic constants for surface complexation modeling (SCM). Although there is a plethora of software to perform data fitting and determine intrinsic equilibrium constants, the algorithms used are highly dependent on initial values and choice of parameters. This limits their transferability to model other systems, for example, reactive transport processes. With this in mind, a hybridized optimization approach, based on a multi–start algorithm combined with a local optimizer, has been developed to allow the simultaneous optimization of SCM parameters and to assess the sensitivity of these parameters to changes in the model assumptions. In this study, the CD–MUSIC formalism with a Basic Stern electrostatic model is utilized to model chromate adsorption on ferrihydrite, although the MUSE algorithm can be applied to any adsorption data set and be implemented in any model formulation. This study offers two innovative components to the inverse SCM modeling approach: a) determination of the true global optimum by performing multiple minimizations of the mean squared error between the simulated and observed data using a large number of initial starting points, and b) quantitative simulation of spectroscopic pH–dependent profiles for two

ACS Paragon Plus Environment

ACS Earth and Space Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

chromate surface complexes. We demonstrate that when MUSE is implemented to determine chromate log Ks, their dependence on other adjustable parameters such as specific surface area (SSA) and capacitance is relatively small (i.e., less than one unit difference for chromate log Ks on ferrihydrite) and can be accounted by mathematical functions determined through the MUSE algorithm. The robustness of the algorithm is demonstrated in the absence of the spectroscopy data as well, with traditional batch tests yielding similar thermodynamic constants as the spectroscopic profiles. Keywords: MUSE, Adsorption, Surface complexation modeling, Optimization, CD - MUSIC, Iron oxides, Chromate.

INTRODUCTION Surface complexation models (SCMs) provide a thermodynamic framework to describe adsorption processes, and can, in principle, replace empirical distribution factors for fate and transport modeling of contaminants.1–5 Several SCM formulations exist, arising from combinations of different pK models (1-pK or 2-pK), single or multisite expressions (charge distribution multisite complexation (CD–MUSIC)),6 and various electrostatic models (constant capacitance (CCM),7 Basic Stern (BS), diffuse layer model (DLM),8,9 and triple layer model (TLM)).4,5 The drawback of using SCMs is the high degree of parameterization, even when pure sorbate-sorbent systems are considered. Parameters related to experimental conditions can be measured, or otherwise constrained, while others are either calculated, assumed or fitted, limiting the transferability of SCMs between systems. Parameters, other than those experimentally determined, are tied to the type of SCM since the underlying assumptions that simulate the solid/liquid interface of each SCM are different.


Page 2 of 42

Page 3 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Specifically, measurable parameters include the solid-to-liquid ratio, ionic strength and solution composition, and the specific surface area (SSA). While the popular N2-BET10,11 method provides accurate SSA values for many materials,12 the difficulty of measuring the SSA on poorly crystallized oxy-hydroxides, such as ferrihydrite, is well established.3,13,14 The uncertainty in SSA values is such that Villalobos and Antelo13 treated it as a normalization parameter in their ferrihydrite study. Site densities (Ns) are also required and constant values such as 3 sites/nm2 or 2.315 sites/nm2 are often used in modeling studies. Crystallographic considerations tied to microscopy tools, such as transmission electron microscopy (TEM) and atomic force microscopy (AFM),16–18 are being increasingly employed to examine individual crystal faces which can exhibit different site densities and binding affinities. Depending on the choice of the electrostatic model, the description of the surface potential may require no electrostatic parameters (DLM), or one to two parameters to describe capacitance (CCM, BS, and TLM), plus additional parameters to describe charge distribution, i.e. the CD factors in the CD–MUSIC model. Capacitance values are treated in various ways; adjusted/ fitted19 to potentiometric titration data, predicted by the radius of the electrolyte cations in TLM,20 predicted as a linear combination of the capacitances of different electric double layer structures,21 obtained from classical molecular dynamics simulations,22 or even determined indirectly using radioisotopic methods.23,24 It is also possible to adopt a Generalized Composite model, coupled with a CCM or a TLM,25 or as often used, with a non-electrostatic model.26 However, this causes equilibrium constants to become more heavily dependent on experimental conditions and site densities.26 Typically, any ligand-specific SCM requires at a minimum four equilibrium constants: protonation (log KH+), electrolyte (log KC+, log KA-), and specific ion/contaminant. If more than



one oxygen type (singly- vs. triply- coordinated oxygens on different locations on the surface) and more than one surface species are considered, the number of necessary thermodynamic constants increases. The surface complexes formed are nowadays usually determined from spectroscopy insights (X-ray absorption near edge structure (XANES), attenuated total reflectance Fourier transform infrared spectroscopy (ATR-FTIR), extended X-Ray Absorption fine Structure (EXAFS)),27–31 mostly in a qualitative manner, although notable exceptions include studies that incorporated quantitative analysis of EXAFS,32,33 nuclear magnetic resonance data (NMR),34 and time resolved laser fluorescence spectroscopy (TRLFS)35 spectra. Protonation equilibrium constants attributed to particular surface oxygen groups may be calculated, as in the case of the MUSIC model,36 or fitted to experimental titration data. Electrolyte and specific ion equilibrium constants are usually fitted to describe potentiometric titrations and macroscopic adsorption data, respectively. A number of fitting algorithms are used to determine surface complexation parameters, including FITEQL,37 GEOSURF,38 ECOSAT-FIT,39,40 and most recently MINFIT.41 These programs allow for optimization of only a single set of parameters during the fitting exercise. More sophisticated fitting routines exist such as PEST,42 which is a generic parameter estimation program. PEST operates as an extension of other geochemical modeling software, but except for coupling with PHREEQC or ORCHESTRA,43 is not in widespread use for surface complexation modeling. Dependence on certain types of software often constraints the modeler to certain types of SCMs and/or a limited choice of optimizable parameters. A user-friendly model for surface complexation parameter optimization to accommodate multiple types of SCMs and numbers of optimizable parameters is yet to be developed.


Page 4 of 42

Page 5 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Several studies have demonstrated the dependence of equilibrium constants and capacitance on site density,44–48 capacitance on specific surface area,44 and equilibrium constants on capacitance.45 Hayes et al.46 showed a decreasing trend of electrolyte equilibrium constants with increasing site densities in various SCMs (CCM, DLM, and TLM), and Hwang and Lenhart44 further illustrated a decreasing trend of capacitance with increasing SSA. Sverjensky49,50 proposed to determine equilibrium constants independent of the solid sorbent properties by normalizing the molar-based equilibrium constants for site density, surface area and solid concentration. However, this approach does not address other interdependencies, such as the correlation of equilibrium constants with capacitance values. Capacitance values are often treated as floating parameters and fitted arbitrarily, leading to non–unique fitting solutions.46. It follows from this analysis that the simultaneous optimization of multiple parameters is necessary to enhance the robustness of SCMs. The goodness of fit to experimental data is often judged subjectively, and not upon specific error criteria. There is little discussion in the literature about global vs local optimization, and non– uniqueness of solutions, with the exception of Jiménez and Mucci51 who employed an evolutionary genetic algorithm to detect global optima. The objective of this study is to address the uncertainty imposed by SCM parameters using a global optimization tool; the MUlti-start optimization algorithm for Surface complexation Equilibrium parameters (MUSE). The MUSE algorithm has been developed to optimize adsorption parameters, such as (but not limited to) equilibrium constants, capacitance values, and SSA. We illustrate its utility by fitting chromate adsorption data to ferrihydrite with a 1-pK Basic Stern CD-MUSIC SCM. The adsorption data include spectroscopic input in the form of multivariate curve resolution with alternating least



squares profiles (MCR–ALS) derived from ATR–FTIR spectra, which further demonstrates the robustness of MUSE. MATERIALS AND METHODS

Figure 1. Flow chart of MUSE algorithm. Experimental data. Two types of data were employed for purposes of model fitting: traditional batch adsorption data collected in this study, and concentration profiles for the monodentate and bidentate species of chromate on ferrihydrite previously presented in Johnston and Chrysochoou.52 Detailed experimental methods for both types of adsorption data are provided in Supporting Information (Figure S1). The utilization of the ATR spectra for modeling purposes required the conversion of the absorbance units to surface coverage units, i.e. mol of adsorbed Cr (VI) per g or m2 of mineral surface. Other studies performed this conversion using Beer’s law,53,54 which requires an estimate for the molar extinction coefficient (scattering coefficient) of the sorbed species. In this study, we utilized MCR–ALS relative distributions of surface complexes based on direct experimental data obtained in the flow cell, as described in Kabengi et al.55 Specifically, the amount of solid on the mineral film deposited on the ATR diamond was


Page 6 of 42

Page 7 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


estimated, and the experimental solid to liquid ratio was adjusted to yield a surface coverage comparable to the batch experiments and flow through experiments collected under similar conditions.55 More details about the conversion of spectroscopic data to surface coverage are given in Supporting Information. Model Description. The overall flowchart of the MUSE algorithm in shown in Figure 1. As with any other modeling approach, it is necessary to choose the formulation of surface reactions and the electrostatics, and here a 1-pK CD-MUSIC model with a Basic Stern layer was adopted. This was preferred because the 1-pK formalism reduces the number of protonation log Ks needed and the Basic Stern layer requires only one capacitance value for fitting. An adjusted capacitance value for spherical particles (1.15 F/m2)56 was utilized for a constant capacitance scenario, whereas appropriate calculations were incorporated into the model when the capacitance value was optimized.56. The algorithm can be modified to accommodate any adsorption data set, surface reaction formulation, and electrostatic model, with limitations imposed mainly by the number of fitting parameters. The fitting problem in this study was to determine the log K values for a ligand with two surface species, i.e. chromate on ferrihydrite. Three parameters are fitted simultaneously, the two log Ks for the surface species and the capacitance, to examine interdependencies between the parameters. Theoretically, the capacitance value depends on the dielectric constant of the solution near the surface and on the radius of adsorbing ions. The dielectric constant is, in turn, a function of the mineral properties and the concentration of charged ions near the surface. Titration curves are, therefore, ideally used to determine the capacitance of a particular mineral and this will vary even for a single surface because of differences in morphology; for example, different ferrihydrite datasets yielded capacitance values between 0.7 and 1.36 F/m2 in a previous



study.57 This capacitance value should then be carried forward to model specific ion adsorption of different ligands. In practice, however, most published adsorption datasets do not have corresponding charging data available and models typically assume capacitance values from the literature. In addition, because capacitance is difficult to determine experimentally, we chose to treat it as a fitted parameter in this study, although we also provide a model scenario which adopts a suggested capacitance value for ferrihydrite from the literature. The algorithm allows for the optimization of more than three parameters at the expense of computing time; depending on the dataset, this could also lead to overfitting and increase the possibility of non-unique solutions. Reducing the number of fitted parameters ultimately requires fixing others, i.e. the conundrum of optimization and transferability to other systems is not completely resolved. However, it is proposed that it can be mitigated. In this study, the choice of fixed versus fitted parameters was made with end users of SCMs in mind, i.e. users of geochemical software (MINTEQ, PHREEQC and others) who will adopt a model from the literature and adjust the parameters to describe their data in micro- and macroscale applications. Typically, adoption of existing SCMs relies on two steps: -

Choosing a database of reactions for a particular surface, which includes protonation, electrolyte and ligand reactions with the associated log Ks.

-

Entering system-specific parameters, including the solid concentration, ionic strength and composition, specific surface area and capacitance(s). Depending on the software, it is possible to choose site densities as well, or these are automatically adopted based on the suggested models within the software.

Addressing the experimental error of adsorption data remains a challenging issue in both experimental and modeling efforts. An adsorption dataset with fully-quantified experimental


Page 8 of 42

Page 9 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


errors due do different sources of variability, such as pH and adsorbate measurements or presence of competitive ions, is currently lacking. Although the uncertainty of experimental error is not directly addressed in this study, the uncertainty related to input parameters is addressed indirectly through sensitivity analysis. The choice of log KH+ values can itself be an optimization problem, if they are determined by fitting surface charge data. While calculated values based on theory have been proposed for iron oxides,19,58 calculation requires the consideration of multiple types of surface oxygens. For example, 11 different sites with log KH+ values ranging from 3 to 12.4 have been proposed for ferrihydrite.58 It is generally impractical to use a large number of variable sites to describe surface charge; Hiemstra et al.36 showed that two sites were adequate to describe the surface charge of goethite, while Bompoti et al.57 used three sites to describe the surface charge of several FH datasets with variable points of zero net proton charge (PZNPCs). The model developed by Bompoti et al.57 for fresh FH with three surface sites is adopted in this study (Table S1), including the log Ks for electrolytes. Briefly, the model consists of one singly coordinated (SC) hydroxyl group and one triply coordinated (TC) group located on the (1-10) and (1-11) planes, and one TC group on the (001) and (00-1) faces. The doubly–coordinated surface groups are not considered in the model since their pK is out of the normal environmental pH range and therefore are considered non-reactive.56 The site densities are taken as the sum of site densities proposed by Hiemstra58 and protonation log Ks are taken to be 8 for the (1-10) and (1-11) faces, and 9.5 for the (001) and (00-1) faces. The purpose of this approach was to allow flexibility in the model to capture differences in charging values on a mass and surface area basis, as well as differences in the PZNPC, which ranges between 8.0 and 8.7.57 These variations appear to be driven by the size of FH crystals, which is influenced by the precipitation method as well as



aging during and after preparation.14 The proposed surface model was able to capture the charging behavior of both fresh and aged FH by increasing the contribution of the basal planes. In this study, we used the site densities for the fresh FH model since the chromate pH edge experiments in this study were performed with fresh FH. Table S2 shows the surface complexation reactions, CD factors and equilibrium constants. CD factors may also be treated as variable parameters and will impact the results for log K values (Figure S2). However, in this study they were based on the structure of the adsorbed ligand, as in many other studies,59 and were fixed at -0.5 and -1.0 for the monodentate and bidentate chromate surface species, respectively. Optimization Algorithm. The model was built with custom-made Mathematica™ notebooks, versions of which have been previously used to simulate the rutile,60–63 magnetite,64 ferrihydrite57 and goethite65 interfaces. They have also been extended to fit adsorption data to 290° C.60,66 The multi-start optimization routine was integrated in the original SCM model to optimize multiple parameters using as criteria the mean squared error (MSE) and Model Selection Criterion (MSC), with larger MSC values indicating a better fit.60,64 Both MSE and MSC depend on the experimental data, specifically the amount and variability of available data, and number of fitting parameters for the latter; thus, they are used to compare different fits for a single data set. The MUSE algorithm also supports the calculation of fit statistics, such as standard deviation and correlations among the parameters. A minimization function was applied, tied to a Nelder-Mead downhill simplex method for optimization. Compared to the other direct search methods available in Mathematica (differential evolution, simulated annealing and random search), Nelder-Mead was free of convergence problems and required less computing time. Nelder–Mead is a nonlinear, derivative-free adaptive


Page 10 of 42

Page 11 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


method applicable for non-linear problems and continuous variables.67,68 For each iteration, a 3dimensional space is formed to minimize the objective function by performing tetrahedral reflections. However, similar to other minimization techniques, Nelder-Mead is dependent on the choice of starting points and boundary conditions, and is often trapped in a local minimum around the initial guesses. To address this problem, the Nelder-Mead method was hybridized with a multi-start algorithm, an approach that has shown promise in other fields for global optimization problems.69–71 Briefly, a random number generator was implemented to produce a 3-dimensional matrix of starting points for the local optimization function. The dimensions of the matrix depend on the number of fitted parameters and their constraints, which was three or four in this case. For each initial point, Nelder-Mead provided an optimal solution; comparison of the solutions using the MSE as selection criterion was implemented to identify the “best local minimum”, which is taken to correspond to the global minimum. From this definition, it follows that finding the true global minimum depends on the number of initial points and parameter boundaries, as defined by the user. RESULTS AND DISCUSSION Algorithm performance. Figure 2a illustrates how the number of starting points influences the MSE value for the spectroscopy dataset, using two different SSA values as constants within the algorithm. The MSE decreases with an increasing number of starting points and plateaus at a minimum value, indicating that the MUSE algorithm successfully detected the global optimum within the selected boundaries for the three parameters. Compared to local optimizers that start with one initial point, the multi – start algorithm provided a solution with an order of magnitude lower MSE (for SSA 600 m2/g, MSE: 2.36*10-13 for one starting point, compared to 2.51*10-14



for starting points > 200, as shown in Figure 2a). Figure 2b plots the evolution in the optimized parameters, showing that the identified global minimum had a much lower capacitance and higher log K values compared to the initial solution with one starting point. The number of points for which the algorithm identified a better minimum varied for different datasets and was even different for the same dataset with different SSA values (Figure 2a). Generally, however, it was observed that starting with more than 200 points was usually sufficient for all datasets when three parameters were fitted simultaneously. The results presented henceforth were all obtained using 500 initial points. Figure S3 and S4 show the distribution of the 500 points used by the algorithm and the local optimization process for three iterations for the spectroscopy dataset, for which the parameter boundaries were chosen to be relatively narrow (5-15 for log K monodentate, 14-23 for log K bidentate and 0.3-3 F/m2 for the capacitance). In this case, the availability of the spectroscopy profiles enabled narrowing of the boundaries for the two log K values, as these were constrained by the respective species distribution. For the batch tests, the boundaries had to be expanded, given that there was no physical constraint for the two log K values, and were 0 to 22 for monodentate and 10 to 30 for bidentate. The boundaries for the capacitance were set to 0.3 – 3 F/m2. This range is wider than the physically plausible values for ferrihydrite (0.7-1.36 F/m2) as shown in a previous study.57 The purpose of optimizing the capacitance values over a wider range was to better evaluate the performance of the algorithm. However, a solution was considered acceptable only for capacitance values within the plausible range for ferrihydrite. The expansion of the parameter boundaries has two implications: a) the 500 points cover a wider domain and are thus less likely to identify the global minimum; and, b) more equivalent solutions


Page 12 of 42

Page 13 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


emerge. These issues will be further illustrated in the discussion of the model results for the spectroscopy and batch datasets.

Figure 2. MUSE results for the FH spectroscopy dataset using two different SSAs: a) MSE as a function of the number of starting points in the multi – start optimization algorithm, and b) extracted parameters (Capacitance (triangles), monodentate (circles) and bidentate (squares) log K values as a function of the number of starting points.

Spectroscopy profiles. Figure 3a shows the model fits for the concentration profiles obtained by MCR-ALS for two SSA values. Along with Table S3, Figure 3 also includes the results of a sensitivity analysis performed for the three fitted parameters as a function of SSA. The extraction of these relationships between the four parameters was enabled by the MUSE algorithm. Other algorithms may be trapped in local minima, away from the optimum solution, providing nonunique thermodynamic parameters. For instance, if the capacitance value is not minimized simultaneously with the equilibrium constants, the dependence of those parameters on SSA is not easily observed.



Figure 3. Adsorption profiles produced from ATR spectroscopy and MCR-ALS (a). Adsorption data (points) and model fits for two SSA values (dashed lines) are shown for total chromate adsorbed and individual surface complexes; corresponding fitted log Ks (b), MSE (c), and capacitance values (d) as a function of SSA. Trend lines represent the relationships for log Ks (c) and Capacitance (d) with SSA: log K MD: -0.664 ln (SSA) + 14.85, log K MD: -0.717 ln (SSA) + 22.215, Capacitance: 108.54 (SSA)-0.736.

This dependence was investigated by conducting a sensitivity analysis with the SSA ranging between 200 and 1200 m2/g; this extends the upper and lower limits of plausible SSAs for fresh FH as reported in the literature. Assuming individual spherical particles, the theoretical SSA of a 1 nm particle would be 1,680 m2/g13, however we did not include such high SSA values since no previous study has reported those. The TEM analysis of the FH used in this study showed particles of 3.4 nm median diameter, which corresponds to a theoretical SSA of 480 m2/g.


Page 14 of 42

Page 15 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Hiemstra et al.72 estimated a SSA of 610 m2/g for a mean particle diameter of 2.5 nm. Aging is also associated with water content decrease and results in a decrease in the available surface area.72 In fact, the uncertainty of SSA is such that a decrease of ~ 100 m2/g was estimated after a 2 hour aging.73 The results of the sensitivity analysis clearly demonstrate a dependence of all three fitted parameters on SSA, and specifically a decrease with increasing SSA. Both log Ks varied in a consistent manner with SSA, i.e. the difference between the bidentate and monodentate log K was approximately constant (Δlog Ks: 6.98 - 7.2), which is necessary for the fitted distribution of the two complexes to closely match the spectroscopically determined distribution. Both equilibrium constants exhibited a range of approximately one log units across the adopted SSA range, decreasing proportionally to the logarithm of SSA (Figure 3b). The bidentate complex exhibits a slightly higher range (1.33 versus 1.1) to account for the double site occupancy. The dependency of log Ks on SSA has implications on the number of sites that are occupied by chromate. By increasing the available SSA, the total site concentration is also increasing, meaning that the ratio of Cr–binding sites to singly coordinated sites decreases. For this spectroscopic dataset, the Cr–binding sites account for 31% of the total sites when SSA is 350 m2/g, while the ratio decreases to 17% for SSA 650 m2/g. A surface area of 69.4 m2/g is the theoretical limiting SSA for the spectroscopic dataset; that is where all the available sites would be occupied by chromate. The inversely proportional change of the capacitance with SSA (Figure 3d) indicates that, mathematically, the capacitance value can account for the majority of the SSA uncertainty when treated as floating parameter. Modeling studies of hematite44 and goethite74 surface charge have shown strong correlations between SSA and capacitance, suggesting higher proton capacity of the lower SSA minerals due presumably to their greater



surface roughness.74 When the surface area increases, the surface loading decreases to account for the same amount of adsorbed chromate, and therefore, the proton co-adsorption is also decreased to maintain a constant proton exchange ratio. This change is reflected in the electrostatics and decrease of capacitance values, as shown in Figure 3d. The SSA and log K values can be further constrained if we consider the plausible range of capacitance values for ferrihydrite; as indicated previously, modelling a wide range of available surface charge data with the same model adopted here yielded capacitance values between 0.7 and 1.36 F/m2.57 As shown in Figure 3d, this range constrains SSA values for the ferrihydrite used in this study between 370 and 850 m2/g and reduces the variability in log K to 0.5 units (10.4 to 10.9) for MD and 0.6 (17.3 to 17.9) BD. For the theoretical value of SSA ~ 600 m2/g the capacitance value is 0.9 F/m2 and the log K MD 10.59 and log K BD 17.57. The sorbent dependency of molar-based equilibrium constants has been previously discussed in a thermodynamic framework by Kulik,75 as have the implications of the molar-based mass action expressions to surface properties.76 The MUSE algorithm accounts for surface activity by incorporating a surface mole fraction scale for surface species, while, in a similar approach, Sverjensky49 proposed the “site occupancy” standard state to account for different sorbent properties. In this model, the equilibrium constants are expressed as a surface mole fraction, in terms of moles/ m2. However, we observe that there is still some dependency on the surface area, although the range of the obtained log Ks is relatively narrow for common SSAs. Wang and Giammar76 have also mentioned findings similar to Kulik, suggesting that even when accounting for surface activities, “denticity” effects on the thermodynamic constants are not consistently addressed. While the discussion of sorbent dependency on intrinsic equilibrium constants has


Page 16 of 42

Page 17 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


focused on multidentate species,76–78 we show that the monodentate thermodynamic constant is also affected by differences in SSA. A benefit of utilizing the spectroscopy profiles is that additional MSEs can be calculated and used as constraints, i.e. the MSEs for the total adsorption envelope and for the individual profiles of the monodentate and bidentate species, which are also shown in Figure 3c. The MSE of the adsorption envelope alone decreased asymptotically with increasing SSA and no true minimum was identified over the range of SSAs tested. Conversely, the individual species profiles both exhibited a clear minimum around 250-350 m2/g and adding all three MSEs the optimal solution emerged at SSA 300 m2/g, with a capacitance of 1.69 F/m2, log KMD 11.01 and log KBD 18.08 (Table S3). The MUSE algorithm was tested for optimization of all four parameters simultaneously, including SSA. The optimal solution coincided with the results of the sensitivity analysis when running 4,000 starting points (SSA 290 m2/g, capacitance of 1.74 F/m2, log KMD 11.03 and log KBD 18.10); having 1,000 initial points yielded an optimal solution at (SSA 327 m2/g, capacitance of 1.55 F/m2, log KMD 10.95 and log KBD 18), i.e. 1000 starting points were not sufficient to find the true minimum. Thus, optimizing more parameters requires a greater number of starting points and increases the risk to model overfitting. However, it is notable that the lower number of points influences the optimal capacitance and SSA values some, but the impact on the optimized log Ks is 0.1 units, i.e. minimal. The optimum solution in terms of the total MSE does not correspond to exactly the BET measured surface area of 347 m2/g, or the theoretical surface area of 600 m2/g that is often used in FH modeling studies. While 290 m2/g is within the range of SSA previously reported BET measured values in the literature,57 it is at the low end, while the optimal capacitance is outside the range of values previously found for ferrihydrite surface charge curves. In other words, the



optimum fit in terms of MSE does not yield a solution that lies within expected values for SSA and capacitance. This begs the question: is the particular dataset an outlier or is the MUSE algorithm yielding values that are not physically plausible? Given that this data is derived from spectroscopic data that have experimental conditions different from traditional batch tests and rely on a number of assumptions to transform the ATR profiles into molar-based concentrations, it is not unlikely that the dataset is atypical. This is the “cost” associated with relying on spectroscopic techniques to obtain independently constrained values for two surface species. The two plausible starting points for SSA (BET and theoretical) also yield solutions with very good fits; in fact, every solution within the typical SSA range for FH (300-650 m2/g) yields a satisfactory fit for both the spectroscopic and the overall adsorption profiles, and capacitance values that are within the range determined by fitting titration data. At this point, there is a tradeoff: a modeler can either pick the optimal solution indicated by the algorithm, or pick one of the solutions that use SSA and capacitance values within the acceptable range of values. And still, the theoretical question remains, which set of log Ks is “true”? Clearly, the strict answer is that “true” log Ks are difficult to determine using any modeling approach. All log K values are a function of the SSA and capacitance values adopted, as well of other parameters that were treated as constants here (e.g., see Figure S2 for effect of CD values). An approach that would best address this issue is to implement functions within geochemical models that will correct log K values for differences in SSA, e.g. in this case the logarithmic functions fitted in Figure 3, as well as calculate or adopt optimal capacitances, which will relieve the end user of the necessity to pick an arbitrary capacitance value. When surface charge data are available, capacitance values can be estimated by fitting those data first. However, charging curve data do not always accompany adsorption data. Parameters that can be determined or


Page 18 of 42

Page 19 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


constrained experimentally should be used as model input and not optimized. Moreover, the functional relationships must be determined for each ligand and surface, which also implies that the functions determined using one specific dataset (the spectroscopic profiles in this study) are transferrable to any other dataset, which is questionable. Another possibility to deal with this uncertainty for modeling studies is to modify the way that models are fitted, reported and transferred into modeling codes. Instead of reporting ligand log Ks as constants, we propose to report average log Ks for a range of SSAs, along with the uncertainty for that range. For FH, SSA typically ranges between 300 and 650 m2/g, considering both BET measurements and theoretical values. The calculated average log Ks for this range are 10.8±0.3 and 17.8±0.3 for monodentate and bidentate chromate surface complexes, respectively. The fluctuations in the log K values are small, and other uncertainties are likely larger, e.g. the experimental errors. An end user may then describe their own dataset using the average log Ks (or choose to modify them within the given range if needed) and fit the SSA and capacitance values using the MUSE algorithm constraining them within the plausible range for each. If a user wants to completely avoid fitting any parameter for forward modeling, they may adopt the values suggested for a specific SSA. For instance, the end user can adopt the theoretical SSA value of 600 m2/g and log Ks of 10.59 and 17.57 for MD and BD, respectively. The capacitance value can be either the optimized value of 0.9 F/ m2 for SSA 600 m2/g or adjusted in the range of expected capacitance values 0.7 to 1.36 F/m2.57 In this case, the capacitance is treated as a floating parameter that is used to account for differences in experimental systems, while keeping the thermodynamic constants and a consistent SSA values. In case the BET value is used as SSA, the log Ks have to be adjusted accordingly. An evaluation of the suggested scenarios for a macroscopic dataset is given below.



Application to Macroscopic data. A conventional batch pH-envelope for FH was also fitted using the MUSE algorithm. The SSA was kept constant in this case, at the theoretical value of 600 m2/g. Several modeling scenarios were performed using the parameters extracted previously. In addition, a modeling scenario with optimized equilibrium constants evaluated the variability of extracted parameters and the model fitting capability without the spectroscopic insights. Since quantifiable spectroscopy data are not available for many ligands, it useful to evaluate the performance of the MUSE algorithm under less constrained conditions. Capacitance values were constrained in the range of plausible values for ferrihydrite (0.7-1.36 F/m2). The following three optimization scenarios were performed: Model with optimized MD and BD log K values: Employment of monodentate and bidentate surface complexes with optimization of both log Ks and capacitance. This scenario was evaluated with and without the consideration of competitive CO2 adsorption, with the CO2 model and results described in Supporting Information. Model Spectroscopy: In this scenario, both surface complexes and log Ks were adopted from the spectroscopic analysis, optimizing only the capacitance value. Specifically, the average calculated log Ks associated with an SSA 600 m2/g were employed as described previously (10.59 and 17.57 for monodentate and bidentate, respectively). Forward model: A simulation of the dataset was also performed using the log Ks and capacitance as extracted from the spectroscopic analysis for SSA 600 m2/g. The only difference with the spectroscopy model is that the capacitance is kept constant. This model was run twice, with a capacitance value of 0.9 F/m2, optimized by MUSE and one with the suggested capacitance of 1.15 F/m2.56


Page 20 of 42

Page 21 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The simulated adsorption data and extracted parameters for the above scenarios are given in Figure 4 and Table 1, respectively. Additional scenarios using the measured BET value of 347.2 m2/g are given in the SI (Figure S5 and Table S4). Utilizing only a visual assessment of the model fits, all three scenarios provide acceptable fits to the experimental data, except in the pH region 7.0-8.0, in which they fail to predict the sharp decrease in adsorption observed in the data. This could only be captured by including CO2 in the model, and both the model with fitted log Ks and the model using the average spectroscopic log Ks provided similar fits (Figure S6 and Table S5). For the model with log Ks optimized, a unique solution emerged where both species are suppressed compared to the spectroscopic values and the capacitance value was optimized to the upper limit of 1.36 F/m2. When the model was optimized at SSA 347.5 m2/g, several equivalent solutions emerged, in terms of both the MSE and the MSC. For both SSA values, the optimal solutions for the two species were outside the one-standard deviation range of the average values determined from the sensitivity analysis. For SSA 600, the MD value was within 2σ and the BD value was within 5σ of the average, both at the lower end. For SSA 347.5, the MD value was several log K units lower and the BD value was higher than the average by 3σ for all equivalent solutions. Clearly, optimization of the log K values for this macroscopic dataset does not yield optimal solutions within the range of the spectroscopic dataset for the two most plausible SSA values. This is an indication that we cannot rely on a single spectroscopic dataset to extrapolate species distribution for batch datasets that have very different experimental conditions.

That said, the two scenarios that used the average values of spectroscopy yielded fits that were fairly satisfactory, if not optimal. The scenario with the adopted values from the spectroscopic



analysis yielded a capacitance value of 0.81 F/m2 and gave a good fit, with the MSC decreasing by 0.3 units, compared to the optimized parameters. A lower MSC is expected for a lower number of optimized parameters, i.e. the more optimized parameters, the better the fitting of the data, but the lower the predictive power of the model. The predictive model (forward scenario) gave an acceptable fit as well (MSC 3.08), only for the MUSE-optimized capacitance of 0.9 F/m2. The scenario with the suggested capacitance value of 1.15 F/m2 yielded a poor fit (MSC 2.58) and overestimated chromate adsorption for the specific dataset. Thus, it may be concluded that although the optimized MUSE values were not the optimal, taken as a whole (the two log K values, SSA and associated optimal capacitance), gave a satisfactory fit with adequate description of the data.

The distribution of the two species are very different in the scenarios using optimized log Ks and the spectroscopic log Ks (Figure S7). When both log Ks are optimized the monodentate species is favored compared to bidentate (Figure S7a). As such, the approach of using the spectroscopic log K values based on SSA 600 m2/g to fit or simulate the dataset provides an equally successful fit to the data, since the spectroscopy data is also matched satisfactorily (Figure S7b and S7c). However, this assumes that the quantitative species distribution in the ATR and batch datasets is similar, despite the difference in experimental conditions.


Page 22 of 42

Page 23 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4. Simulated macroscopic adsorption data for FH at ionic strength 0.01 M and SSA 600 m2/g. Data (points) and model simulation (dashed lines) are shown for adsorbed chromate.

An additional scenario considering an outer-sphere complex instead of a monodentate species was evaluated. Although, spectroscopy does not clearly support the presence of an outer-sphere complex, DFT calculations on FH have shown more energetically favored outer-sphere complexes compared to monodentate complexes, while inner-sphere bidentate complexes are predicted to be the most energetically favored.55. For that reason, a model scenario including an outer-sphere complex instead of a monodentate complex was examined, yielding an equivalent fit to the model with optimized both monodentate and bidentate log Ks (MSC 3.18). The optimized equilibrium constant for the outer–sphere complex was log Kouter-sphere 2.01, while the bidentate equilibrium constant remained unchanged at log KBD 17.86. Model simulation for the outer sphere case is shown in Figure S8. Clearly, the monodentate inner-sphere and the outersphere complexes are equivalent from a pure modeling perspective. A difficulty in connecting the models developed using spectroscopy and batch data is that spectroscopic analyses are performed under very different experimental conditions. ATR-FTIR uses a very small mass of mineral and a very high volume of solution. Ultimately, it is up to the end user to decide the priority, i.e. choose the best model in terms of fit alone, or choose the best model in terms of fit



and agreement with spectroscopy. The application of the MUSE algorithm for chromate on ferrihydrite showed that both approaches can yield fits that are highly satisfactory, if not completely equivalent. Additional macroscopic data under different conditions would help to constrain the model, and those will be discussed in a subsequent study. Table 1. Optimized parameters and their standard deviation for batch chromate pH – envelope dataset under different modeling scenarios described in text. Number of C Model parameters log K MD log K BD MSC MSE (F/m2) optimized Optimized log Ks* 3 10.24 16.45 1.36 3.67 7.604E-10 Spectroscopy** 1 10.59 17.57 0.81 3.39 1.011E-09 Forward model*** 0 10.59 17.57 0.91 3.08 1.369E-09 Forward model*** 0 10.59 17.57 1.15 2.58 2.27E-09 *Log K MD, log K BD, and C were optimized in this scenario **Log K MD and log K BD were adopted from spectroscopic analysis for SSA 600 m2/g ***All parameters were adopted from spectroscopic analysis for SSA 600 m2/g CONCLUSIONS SCMs have been developed and applied to interpret and rationalize experimental adsorption data for the past 50 years. However, the transferability of specific SCMs to other experimental data sets and ultimately to reactive transport modeling is constrained by several lingering problems, including the fact that SCMs have a large number of interdependent parameters. While increased parametrization enhances SCM flexibility, it also increases the uncertainty in thermodynamic constants and end users of geochemical models lack the tools to choose and implement a selfconsistent model. Parameters such as those related to electrostatics are not easily constrained, while others (site densities and specific surface area) can be constrained by physical models or measurements. Ultimately, thermodynamic constants are dependent on all these values, a fact that end users of geochemical modeling software cannot easily account for. An experienced user


Page 24 of 42

Page 25 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


can inform the SCM with good initial guesses and find parameters very close to the global optimum ones, but not all users are sufficiently experienced. The MUSE algorithm addresses this experience deficit by allowing the simultaneous optimization of several parameters and enabling the determination of a global minimum that is not constrained by the initial guesses for the parameter values. It can be applied to any adsorption data set, and any type of SCM and choice of parameters, including equilibrium constants, capacitance values and sorbent properties. We assert that when MUSE is implemented to determine adsorption log Ks, their dependence on other adjustable parameters such as SSA and capacitance will be relatively small, and certainly smaller compared to traditional single-point optimization. Thus, the MUSE algorithm facilitates parameter optimization and reduces parameter interdependence which is critical for adapting SCMs to diverse data sets and ultimately reactive transport models. The utility of the MUSE algorithm was demonstrated by fitting a set of batch and spectroscopic adsorption data for chromate adsorption by ferrihydrite with a CD-MUSIC SCM. Mean square error values for model fits decreased by an order of magnitude, and best-fit parameter values (capacitance, and monodentate and bidentate log K values) changed significantly when multiple starting parameter values were utilized within the MUSE algorithm. Moreover, the MCR-ALS profiles derived from ATR-FTIR spectra, which provided the relative distribution between chromate adsorbed as monodentate and bidentate species, were fit directly, thereby constraining the fit log Ks. Such molecular-level insight will be extremely useful in developing a consistent set of equilibrium constants for interpreting other adsorption data sets and reactive transport modeling purposes.



ASSOCIATED CONTENT Supporting Information. Experimental techniques for mineral synthesis, batch adsorption Phenvelopes and ATR–FTIR conversion to surface coverage. Detailed description of the Surface Complexation model, including the FH surface structure and surface complexation reactions, as well as the MUSE optimization process, modeling scenarios including evaluation of the model under different SSA, carbonate competitions and presence of outer-sphere complexes, speciation profiles of inner–sphere complexes for all modeling scenarios, and CD factors and SSA sensitivity analysis results. AUTHOR INFORMATION Corresponding Author *Phone: +1860 771 8519, E-mail: [email protected] Notes The authors declare no competing financial interest. ACKNOWLEDGMENT This work has been supported by the National Science Foundation Environmental Chemical Sciences program under Award Number CHE-1404643 entitled “Collaborative Research: Toward a unified model for ferrihydrite nanoparticles behavior in the environment: a multipronged investigation of surface structure and reactivity”. The TEM work was supported by the award of an FEI Graduate Fellowship.


Page 26 of 42

Page 27 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


REFERENCES (1) Goldberg, S.; Criscenti, L. J.; Turner, D. R.; Davis, J. A.; Cantrell, K. J. Adsorption– Desorption Processes in Subsurface Reactive Transport Modeling. Vadose Zone J. 2007, 6 (3), 407. (2) Stumm, W.; Kummert, R.; Sigg, L. A Ligand-Exchange Model for the Adsorption of Inorganic and Organic-Ligands at Hydrous Oxide Interfaces. Croat. Chem. Acta 1980, 53 (2), 291–312. (3) Dzombak, D. A.; Morel, F. M. M. Surface Complexation Modeling: Hydrous Ferric Oxide; John Wiley & Sons, 1990. (4) Yates, D. E.; Levine, S.; Healy, T. W. Site-Binding Model of the Electrical Double Layer at the Oxide/Water Interface. J. Chem. Soc. Faraday Trans. 1 Phys. Chem. Condens. Phases 1974, 70 (0), 1807. (5) Davis, J. A.; James, R. O.; Leckie, J. O. Surface Ionization and Complexation at the Oxide/Water Interface. J. Colloid Interface Sci. 1978, 63 (3), 480–499. (6) Hiemstra, T.; Van Riemsdijk, W. H. A Surface Structural Approach to Ion Adsorption: The Charge Distribution (CD) Model. J. Colloid Interface Sci. 1996, 179 (2), 488–508. (7) Stumm, W.; Hohl, H.; Dalang, F. Interaction of Metal-Ions with Hydrous Oxide Surfaces. Croat. Chem. Acta 1976, 48 (4), 491–504. (8) Stumm, W.; Huang, C. P.; Jenkins, S. R. Specific Chemical Interaction Affecting Stability of Dispersed Systems. Croat. Chem. Acta 1970, 42 (2), 223. (9) Huang, C.-P.; Stumm, W. Specific Adsorption of Cations on Hydrous γ-Al2O3. J. Colloid Interface Sci. 1973, 43 (2), 409–420.



(10) Brunauer, S.; Emmett, P. H.; Teller, E. Adsorption of Gases in Multimolecular Layers. J. Am. Chem. Soc. 1938, 60 (2), 309–319. (11) Gregg, S. J.; Sing, K. S. W.; Salzberg, H. W. Adsorption Surface Area and Porosity. J. Electrochem. Soc. 1967, 114 (11), 279C-279C. (12) Everett, D. H.; Parfitt, G. D.; Sing, K. S. W.; Wilson, R. The SCI/IUPAC/NPL Project on Surface Area Standards. J. Appl. Chem. Biotechnol. 1974, 24 (4–5), 199–219. (13) Villalobos, M.; Antelo, J. A unified structural model for ferrihydrite: Proton charge, electrolyte binding, and arsenate adsorption. Revista Internacional de Contaminación Ambiental 2011, 27 (2). (14) Villacís-García, M.; Ugalde-Arzate, M.; Vaca-Escobar, K.; Villalobos, M.; Zanella, R.; Martínez-Villegas, N. Laboratory Synthesis of Goethite and Ferrihydrite of Controlled Particle Sizes. Bol. Soc. Geológica Mex. 2015, 67 (3), 433. (15) Davis, J. A.; Kent, D. B. Surface Complexation Modeling in Aqueous Geochemistry. Rev. Mineral. Geochem. 1990, 23 (1), 177–260. (16) Gaboriaud, F.; Ehrhardt, J.-J. Effects of Different Crystal Faces on the Surface Charge of Colloidal Goethite (α-FeOOH) Particles: An Experimental and Modeling Study. Geochim. Cosmochim. Acta 2003, 67 (5), 967–983. (17) Villalobos, M.; Cheney, M. A.; Alcaraz-Cienfuegos, J. Goethite Surface Reactivity: II. A Microscopic Site-Density Model That Describes Its Surface Area-Normalized Variability. J. Colloid Interface Sci. 2009, 336 (2), 412–422. (18) Livi, K. J. T.; Villalobos, M.; Leary, R.; Varela, M.; Barnard, J.; Villacís-García, M.; Zanella, R.; Goodridge, A.; Midgley, P. Crystal Face Distributions and Surface Site


Page 28 of 42

Page 29 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Densities of Two Synthetic Goethites: Implications for Adsorption Capacities as a Function of Particle Size. Langmuir 2017, 33 (36), 8924–893. (19) Venema, P.; Hiemstra, T.; Weidler, P. G.; van Riemsdijk, W. H. Intrinsic Proton Affinity of Reactive Surface Groups of Metal (Hydr)Oxides: Application to Iron (Hydr)Oxides. J. Colloid Interface Sci. 1998, 198 (2), 282–295. (20) Sverjensky, D. A. Interpretation and Prediction of Triple-Layer Model Capacitances and the Structure of the Oxide-Electrolyte-Water Interface. Geochim. Cosmochim. Acta 2001, 65 (21), 3643–3655. (21) Boily, J.-F. The Variable Capacitance Model: A Strategy for Treating Contrasting ChargeNeutralizing Capabilities of Counterions at the Mineral/Water Interface. Langmuir 2014, 30 (8), 2009–2018. (22) Parez, S.; Předota, M.; Machesky, M. Dielectric Properties of Water at Rutile and Graphite Surfaces: Effect of Molecular Structure. J. Phys. Chem. C 2014, 118 (9), 4818–4834. (23) Sprycha, R. Electrical Double Layer at Alumina/Electrolyte Interface: I. Surface Charge and Zeta Potential. J. Colloid Interface Sci. 1989, 127 (1), 1–11. (24) Sprycha, R. Electrical Double Layer at Alumina/Electrolyte Interface: II. Adsorption of Supporting Electrolyte Ions. J. Colloid Interface Sci. 1989, 127 (1), 12–25. (25) Goldberg, S. Application of Surface Complexation Models to Anion Adsorption by Natural Materials. Environ. Toxicol. Chem. 2014, 33 (10), 2172–2180. (26) Stoliker, D. L.; Kent, D. B.; Zachara, J. M. Quantifying Differences in the Impact of Variable Chemistry on Equilibrium Uranium(VI) Adsorption Properties of Aquifer Sediments. Environ. Sci. Technol. 2011, 45 (20), 8733–8740.



(27) Tiberg, C.; Sjöstedt, C.; Persson, I.; Gustafsson, J. P. Phosphate Effects on Copper(II) and Lead(II) Sorption to Ferrihydrite. Geochim. Cosmochim. Acta 2013, 120, 140–157. (28) Fukushi, K.; Aoyama, K.; Yang, C.; Kitadai, N.; Nakashima, S. Surface Complexation Modeling for Sulfate Adsorption on Ferrihydrite Consistent with in Situ Infrared Spectroscopic Observations. Appl. Geochem. 2013, 36, 92–103. (29) Rossberg, A.; Ulrich, K.-U.; Weiss, S.; Tsushima, S.; Hiemstra, T.; Scheinost, A. C. Identification of Uranyl Surface Complexes on Ferrihydrite: Advanced EXAFS Data Analysis and CD-MUSIC Modeling. Environ. Sci. Technol. 2009, 43 (5), 1400–1406. (30) Villalobos, M.; Leckie, J. O. Surface Complexation Modeling and FTIR Study of Carbonate Adsorption to Goethite. J. Colloid Interface Sci. 2001, 235 (1), 15–32. (31) Ponthieu, M.; Juillot, F.; Hiemstra, T.; van Riemsdijk, W. H.; Benedetti, M. F. Metal Ion Binding to Iron Oxides. Geochim. Cosmochim. Acta 2006, 70 (11), 2679–2698. (32) Gu, C.; Wang, Z.; Kubicki, J. D.; Wang, X.; Zhu, M. X-Ray Absorption Spectroscopic Quantification and Speciation Modeling of Sulfate Adsorption on Ferrihydrite Surfaces. Environ. Sci. Technol. 2016, 50 (15), 8067–8076. (33) Wang, X.; Wang, Z.; Peak, D.; Tang, Y.; Feng, X.; Zhu, M. Quantification of Coexisting Inner- and Outer-Sphere Complexation of Sulfate on Hematite Surfaces - ACS Earth and Space Chemistry (ACS Publications). ACS Earth Space Chem. 2018, 2 (4), 387–398. (34) Ingri, N.; Andersson, I.; Pettersson, L.; Yagasaki, A.; Andersson, L.; Holmström, K.; Undheim, K.; Rosendahl, C. N.; Haugg, M.; Trabesinger-Rüf, N.; et al. LAKE -- A Program System for Equilibrium Analytical Treatment of Multimethod Data, Especially Combined Potentiometric and Nuclear Magnetic Resonance Data. Acta Chem. Scand. 1996, 50, 717–734.


Page 30 of 42

Page 31 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(35) Stumpf, S.; Stumpf, T.; Lützenkirchen, J.; Walther, C.; Fanghänel, T. Immobilization of Trivalent Actinides by Sorption onto Quartz and Incorporation into Siliceous Bulk: Investigations by TRLFS. J. Colloid Interface Sci. 2008, 318 (1), 5–14. (36) Hiemstra, T.; Venema, P.; Riemsdijk, W. H. V. Intrinsic Proton Affinity of Reactive Surface Groups of Metal (Hydr)Oxides: The Bond Valence Principle. J. Colloid Interface Sci. 1996, 184 (2), 680–692. (37) Oregon State University., J. C.; Westall. FITEQL : A Computer Program for Determination of Chemical Equilibrium Constants from Experimental Data; Dept. of Chemistry, Oregon State University,: Corvallis, Or. :, 1982. (38) Sahai, N.; Sverjensky, D. A. GEOSURF: A Computer Program for Modeling Adsorption on Mineral Surfaces from Aqueous Solution. Comput. Geosci. 1998, 24 (9), 853–873. (39) Keizer, M.; an Riemsdijk, W. H. ECOSAT (Ver. 4.7) Equilibrium Calculation of Speciation and Transport. Wagening. Agric. Univ. Wagening. 1999. (40) Kinniburgh, D. G. FIT Non-Linear Optimization Algorithm and User Manual. Br. Geol. Surv. Nottm. UK 1993. (41) Xie, X.; Giammar, D. E.; Wang, Z. MINFIT: A Spreadsheet-Based Tool for Parameter Estimation in an Equilibrium Speciation Software Program. Environ. Sci. Technol. 2016, 50 (20), 11112–11120. (42) Doherty, J. PEST: A Unique Computer Program for Model-Independent Parameter Optimisation. Water 94 GroundwaterSurface Hydrol. Common Interest Pap. Prepr. Pap. 1994, 551. (43) Janot, N.; Pinheiro, J. P.; Botero, W. G.; Meeussen, J. C. L.; Groenenberg, J. E. PESTORCHESTRA, a Tool for Optimising Advanced Ion-Binding Model Parameters:



Derivation of NICA-Donnan Model Parameters for Humic Substances Reactivity. Environ. Chem. 2016, 14 (1), 31–38. (44) Hwang, Y. S.; Lenhart, J. J. The Dependence of Hematite Site-Occupancy Standard State Triple-Layer Model Parameters on Inner-Layer Capacitance. J. Colloid Interface Sci. 2008, 319 (1), 206–213. (45) Katz, L. E.; Hayes, K. F. Surface Complexation Modeling. J. Colloid Interface Sci. 1995, 170 (2), 477–490. (46) Hayes, K. F.; Redden, G.; Ela, W.; Leckie, J. O. Surface Complexation Models: An Evaluation of Model Parameter Estimation Using FITEQL and Oxide Mineral Titration Data. J. Colloid Interface Sci. 1991, 142 (2), 448–469. (47) Christl, I.; Kretzschmar, R. Competitive Sorption of Copper and Lead at the Oxide-Water Interface: Implications for Surface Site Density. Geochim. Cosmochim. Acta 1999, 63 (19– 20), 2929–2938. (48) Goldberg, S. Sensitivity of Surface Complexation Modeling to the Surface Site Density Parameter. J. Colloid Interface Sci. 1991, 145 (1), 1–9. (49) Sverjensky, D. A. Standard States for the Activities of Mineral Surface Sites and Species. Geochim. Cosmochim. Acta 2003, 67 (1), 17–28. (50) Sverjensky, D. A. Prediction of Surface Charge on Oxides in Salt Solutions: Revisions for 1:1 (M+L−) Electrolytes. Geochim. Cosmochim. Acta 2005, 69 (2), 225–257. (51) Villegas-Jiménez, A.; Mucci, A. Estimating Intrinsic Formation Constants of Mineral Surface Species Using a Genetic Algorithm. Math. Geosci. 2010, 42 (1), 101.


Page 32 of 42

Page 33 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(52) Johnston, C. P.; Chrysochoou, M. Mechanisms of Chromate, Selenate, and Sulfate Adsorption on Al-Substituted Ferrihydrite: Implications for Ferrihydrite Surface Structure and Reactivity. Environ. Sci. Technol. 2016, 50 (7), 3589–3596. (53) Sabur, M. A.; Goldberg, S.; Gale, A.; Kabengi, N.; Al-Abadleh, H. A. TemperatureDependent Infrared and Calorimetric Studies on Arsenicals Adsorption from Solution to Hematite Nanoparticles. Langmuir 2015, 31 (9), 2749–2760. (54) Mitchell, W.; Goldberg, S.; Al-Abadleh, H. A. In Situ ATR–FTIR and Surface Complexation Modeling Studies on the Adsorption of Dimethylarsinic Acid and p-Arsanilic Acid on Iron-(Oxyhydr)Oxides. J. Colloid Interface Sci. 2011, 358 (2), 534–540. (55) Kabengi, N. J.; Chrysochoou, M.; Bompoti, N.; Kubicki, J. D. An Integrated Flow Microcalorimetry, Infrared Spectroscopy and Density Functional Theory Approach to the Study of Chromate Complexation on Hematite and Ferrihdyrite. Chem. Geol. 2017, 464, 23–33. (56) Hiemstra, T.; Van Riemsdijk, W. H. A Surface Structural Model for Ferrihydrite I: Sites Related to Primary Charge, Molar Mass, and Mass Density. Geochim. Cosmochim. Acta 2009, 73 (15), 4423–4436. (57) Bompoti, N.; Chrysochoou, M.; Machesky, M. Surface Structure of Ferrihydrite: Insights from Modeling Surface Charge. Chem. Geol. 2017, 464, 34–45. (58) Hiemstra, T. Surface and Mineral Structure of Ferrihydrite. Geochim. Cosmochim. Acta 2013, 105, 316–325. (59) Rietra, R. P. J. J.; Hiemstra, T.; van Riemsdijk, W. H. The Relationship between Molecular Structure and Ion Adsorption on Variable Charge Minerals. Geochim. Cosmochim. Acta 1999, 63 (19), 3009–3015.



(60) Machesky, M. L.; Wesolowski, D. J.; Palmer, D. A.; Ichiro-Hayashi, K. Potentiometric Titrations of Rutile Suspensions to 250°C. J. Colloid Interface Sci. 1998, 200 (2), 298–309. (61) Hawkins, T.; Allen, N.; Machesky, M. L.; Wesolowski, D. J.; Kabengi, N. Ion Exchange Thermodynamics at the Rutile–Water Interface: Flow Microcalorimetric Measurements and Surface Complexation Modeling of Na–K–Rb–Cl–NO3 Adsorption. Langmuir 2017, 33 (20), 4934–4941. (62) Machesky, M. L.; Předota, M.; Ridley, M. K.; Wesolowski, D. J. Constrained Surface Complexation Modeling: Rutile in RbCl, NaCl, and NaCF3SO3 Media to 250 °C. J. Phys. Chem. C 2015, 119 (27), 15204–15215. (63) Ridley, M. K.; Machesky, M. L.; Wesolowski, D. J.; Palmer, D. A. Surface Complexation of Neodymium at the Rutile-Water Interface: A Potentiometric and Modeling Study in NaCl Media to 250°C. Geochim. Cosmochim. Acta 2005, 69 (1), 63–81. (64) Wesolowski, D. J.; Machesky, M. L.; Palmer, D. A.; Anovitz, L. M. Magnetite Surface Charge Studies to 290°C from in Situ PH Titrations. Chem. Geol. 2000, 167 (1–2), 193– 229. (65) Machesky, M. L.; Andrade, W. O.; Rose, A. W. Adsorption of Gold(III)-Chloride and Gold(I)-Thiosulfate Anions by Goethite. Geochim. Cosmochim. Acta 1991, 55 (3), 769– 776. (66) Machesky, M. L.; Wesolowski, D. J.; Palmer, D. A.; Ridley, M. K. On the Temperature Dependence of Intrinsic Surface Protonation Equilibrium Constants: An Extension of the Revised MUSIC Model. J. Colloid Interface Sci. 2001, 239 (2), 314–327. (67) Nelder, J. A.; Mead, R. A Simplex Method for Function Minimization. Comput. J. 1965, 7 (4), 308–313.


Page 34 of 42

Page 35 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(68) Olsson, D. M.; Nelson, L. S. The Nelder-Mead Simplex Procedure for Function Minimization. Technometrics 1975, 17 (1), 45–51. (69) Kocsis, L.; György, A. Efficient Multi-Start Strategies for Local Search Algorithms. In Machine Learning and Knowledge Discovery in Databases; Springer, Berlin, Heidelberg, 2009; pp 705–720. (70) Huang, M.; Aine, C. J.; Supek, S.; Best, E.; Ranken, D.; Flynn, E. R. Multi-Start Downhill Simplex Method for Spatio-Temporal Source Localization in Magnetoencephalography. Electroencephalogr. Clin. Neurophysiol. Potentials Sect. 1998, 108 (1), 32–44. (71) Martí, R.; Resende, M. G. C.; Ribeiro, C. C. Multi-Start Methods for Combinatorial Optimization. Eur. J. Oper. Res. 2013, 226 (1), 1–8. (72) Hiemstra, T. Ferrihydrite Interaction with Silicate and Competing Oxyanions: Geometry and Hydrogen Bonding of Surface Species. Geochim. Cosmochim. Acta 2018, 238, 453– 476. (73) Hiemstra, T. Formation, Stability, and Solubility of Metal Oxide Nanoparticles: Surface Entropy, Enthalpy, and Free Energy of Ferrihydrite. Geochim. Cosmochim. Acta 2015, 158, 179–198. (74) Boily, J.-F.; Lützenkirchen, J.; Balmès, O.; Beattie, J.; Sjöberg, S. Modeling Proton Binding at the Goethite (α-FeOOH)–Water Interface. Colloids Surf. Physicochem. Eng. Asp. 2001, 179 (1), 11–27. (75) Kulik, D. A. Gibbs Energy Minimization Approach to Modeling Sorption Equilibria at the Mineral-Water Interface: Thermodynamic Relations for Multi-Site-Surface Complexation. Am. J. Sci. 2002, 302 (3), 227–279.



(76) Wang, Z.; Giammar, D. E. Mass Action Expressions for Bidentate Adsorption in Surface Complexation Modeling: Theory and Practice. Environ. Sci. Technol. 2013, 47 (9), 3982– 3996. (77) Benjamin, M. M. Modeling the Mass-Action Expression for Bidentate Adsorption. Environ. Sci. Technol. 2002, 36 (3), 307–313. (78) Lützenkirchen, J.; Marsac, R.; Kulik, D. A.; Payne, T. E.; Xue, Z.; Orsetti, S.; Haderlein, S. B. Treatment of Multi-Dentate Surface Complexes and Diffuse Layer Implementation in Various Speciation Codes. Appl. Geochem. 2015, 55, 128–137.


Page 36 of 42

Page 37 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Abstract graphic



Figure 1. Flow chart of MUSE algorithm. 178x80mm (300 x 300 DPI)


Page 38 of 42

Page 39 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. MUSE results for the FH spectroscopy dataset using two different SSAs: a) MSE as a function of the number of starting points in the multi – start optimization algorithm, and b) extracted parameters (Capacitance (triangles), monodentate (circles) and bidentate (squares) log K values as a function of the number of starting points. 50x31mm (300 x 300 DPI)



Figure 2. MUSE results for the FH spectroscopy dataset using two different SSAs: a) MSE as a function of the number of starting points in the multi – start optimization algorithm, and b) extracted parameters (Capacitance (triangles), monodentate (circles) and bidentate (squares) log K values as a function of the number of starting points. 50x31mm (300 x 300 DPI)


Page 40 of 42

Page 41 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3. Adsorption profiles produced from ATR spectroscopy and MCR-ALS (a). Adsorption data (points) and model fits for two SSA values (dashed lines) are shown for total chromate adsorbed and individual surface complexes; corresponding fitted log Ks (b), MSE (c), and capacitance values (d) as a function of SSA. Trend lines represent the relationships for log Ks (c) and Capacitance (d) with SSA: log K MD: -0.664 ln (SSA) + 14.85, log K MD: -0.717 ln (SSA) + 22.215, Capacitance: 108.54 (SSA)-0.736. 154x95mm (300 x 300 DPI)



Figure 4. Simulated macroscopic adsorption data for FH at ionic strength 0.01 M and SSA 600 m2/g. Data (points) and model simulation (dashed lines) are shown for adsorbed chromate. 81x50mm (300 x 300 DPI)


Page 42 of 42

Assessment of Modeling Uncertainties Using a Multistart Optimization

Recommend Documents