MINFIT: A Spreadsheet-Based Tool for Parameter Estimation in an

Subscriber access provided by Heriot-Watt | University Library

Article

MINFIT: A Spreadsheet-based Tool for Parameter Estimation in an Equilibrium Speciation Software Program Xiongfei Xie, Daniel E. Giammar, and Zimeng Wang Environ. Sci. Technol., Just Accepted Manuscript • DOI: 10.1021/acs.est.6b03399 • Publication Date (Web): 23 Sep 2016 Downloaded from http://pubs.acs.org on September 24, 2016

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 free 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 accessible to all readers and 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.

Environmental Science & Technology 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 33

Environmental Science & Technology

MINFIT: A Spreadsheet-based Tool for Parameter Estimation in an Equilibrium Speciation Software Program Xiongfei Xie†, Daniel E. Giammar‡ and Zimeng Wang§* † City of Lakeland Water Utilities Department, Lakeland, Florida, United States ‡ Department of Energy, Environmental and Chemical Engineering, Washington University in St. Louis, St. Louis, Missouri, United States § Department of Civil and Environmental Engineering, Louisiana State University, Baton Rouge, Louisiana, United States

*Corresponding author [email protected] 3316H Patrick Taylor Hall, Louisiana State University, Baton Rouge, LA 70803 Phone: (225) 578-1591 / Fax: (225) 578-4945 MINFIT Website: http://minfit.strikingly.com

Revised Manuscript Submitted to Environmental Science & Technology September 2016

ACS Paragon Plus Environment


1

Page 2 of 33

Abstract

2

Determination of equilibrium constants describing chemical reactions in the aqueous phase

3

and at solid-water interface relies on inverse modeling and parameter estimation. Although there

4

are existing tools available, the steep learning curve prevents the wider community of

5

environmental engineers and chemists to adopt those tools. Stemming from classical chemical

6

equilibrium codes, MINEQL+ has been one of the most widely used chemical equilibrium

7

software programs. We developed a spreadsheet-based tool, which we are calling MINFIT, that

8

interacts with MINEQL+ to perform parameter estimations that optimize model fits to

9

experimental datasets. MINFIT enables automatic and convenient screening of a large number of

10

parameter sets towards the optimal solutions by calling MINEQL+ to perform iterative forward

11

calculations following either exhaustive equidistant grid search or randomized search algorithms.

12

The combined use of the two algorithms can securely guide the searches for the global optima. We

13

developed interactive interfaces so that the optimization processes are transparent. Benchmark

14

examples including both aqueous and surface complexation problems illustrate the parameter

15

estimation and associated sensitivity analysis. MINFIT is accessible at http://minfit.strikingly.com.

16

1


Page 3 of 33

17


Introduction

18

Chemical equilibrium models calculate speciation of elements by simultaneously solving mole

19

balance and mass action equations. The applications of the models include both forward

20

calculations with known equilibrium constants and inverse calculations to determine the unknown

21

constants using experimental data. Determination of intrinsic equilibrium constants describing

22

chemical reactions in the aqueous phase and at solid-water interfaces is essential for developing

23

speciation-based models for element mobility in aquatic systems.1 In the last few decades, a

24

number of computer programs have been developed and popularized to varied extent and

25

geographical coverage in the aquatic chemistry community2, including MINEQL+3, MINSORB,4

26

FITEQL5, MINTEQ6, PHREEQC7, ECOSAT8, ORCHESTRA9, GEOSURF10, and Geochemist’s

27

Workbench.11 While there are various modeling software programs, the fundamental numerical

28

algorithm (the tableau approach) is similar for a vast majority of them.

29

As an advanced category of chemical equilibrium models and an aspect of aquatic chemistry

30

in which parameter estimation is very common, surface complexation models (SCM) simulate

31

adsorption by considering electrostatic effects and can automatically modify the equilibrium

32

constants according to various relationships between surface potential and surface charge.12,

33

Compared with the conventional distribution coefficient (Kd) and empirical isotherm models (e.g.,

34

Langmuir and Freundlich), SCM can predict the impact of solution chemistry on the binding of 2


13


Page 4 of 33

35

aqueous solutes to surfaces with a single set of parameters over a broad range of pH and other

36

conditions.14, 15

37

Forward modeling with known equilibrium constants and total concentrations is enabled by

38

default in those programs. In contrast, inverse modeling (calibration of models from experimental

39

data) has been less straightforward. FITEQL, a derivative-based nonlinear least squares

40

optimization program initially developed in the early 1980s and occasionally updated until the late

41

1990s, is a widely used inverse modeling software program for this purpose. The advantage of

42

FITEQL is that the user interface allows full customizability and transparency. However, FITEQL

43

does not have a built-in thermodynamic database or activity coefficient model. The requirement of

44

providing initial guesses for the fitting parameters in FITEQL and other similar programs may

45

lead to non-unique solutions or solutions representing local rather than global minima. Additional

46

pain points for users are (1) the troubleshooting for non-convergence that often occurs for poorly

47

defined systems with multiple fitting parameters (usually when two or more are to be determined16,

48

17

49

that FITEQL convergence could not be achieved due to various reasons.17-21 In some cases,

50

repetitive manual trial-and-error fitting using multiple forward calculation had to be conducted.22-

51

27

), and (2) the re-calculation of ionic strength dependent equilibrium constants. It was common

The most recent operating system that supports FITEQL (Windows XP) is being phased out.

3


Page 5 of 33


52

We are aware that in addition to FITEQL, other programs such as ECOSAT8 (in combination with

53

a program named FIT28), ORCHESTRA9 and GEOSURF10 have various capabilities for parameter

54

estimation. Although those codes are powered with advanced functions, their dissemination has

55

only extended to limited geographic regions and academic groups. ProtoFit, an easier-to-use

56

alternative of FITEQL, is powerful for optimization problems for protonation reactions on sorbent

57

surfaces from acid-base titration,29 but it is not ready for full models describing ion adsorption.30

58

In comparison, those programs with interactive user interfaces and widely distributed application

59

examples, represented by MINEQL+, Visual MINTEQ and PHREEQC, are popular in a wide user

60

community, which is reflected by their thousands of citations compared with a few dozens for the

61

other mentioned programs, as of 2016.

62

Connecting external parameter estimation programs (such as PEST31, UCODE32 and

63

homemade codes in Matlab or Excel) with geochemical codes can be performed to estimate any

64

parameters.33-37 However, implementing such a connection requires fairly advanced knowledge

65

about data I/O (input/output) and the structure of those programs. There are also successful

66

developments of advanced algorithms for inverse modeling (e.g., genetic algorithms30, 38, Gibbs

67

Energy Minimization39, 40). Such applications are mostly employed by those who have programing

68

experience. The most recent version of Visual MINTEQ (3.1) has a module that generates files to

4



Page 6 of 33

69

be read by PEST. The module solves the optimization problem with iterative calls of Visual

70

MINTEQ by PEST until convergence. The Gauss-Marquardt-Levenberg algorithm adopted by

71

PEST also relies on good initial guesses. PHREEQC allows “tagged” parameters to change when

72

called by external programs, providing a means to optimize selected parameters through the fitting

73

routine. A recent program PhreePlot41 used in conjunction with PHREEQC enables the fitting

74

functionality similarly as UCODE does and offers several algorithms for optimization.42 However,

75

its powerful functionality also means higher learning curve for users who are not familiar with

76

PHREEQC and command line interfaces.

77

From our interactions with colleagues and our own experiences, we believe that there is a need

78

for easy-to-use complementary tools that perform parameter estimation for the most widely used

79

chemical equilibrium modeling programs without requiring advanced knowledge of programming.

80

Among the existing software programs, MINEQL+ has basic functionality and is very popular

81

with relatively easy introduction to a wide range of users.43 Unfortunately MINEQL+, including

82

its most recent version 5.0 (November 2015), does not contain a parameter estimation capability.

83

MINEQL+ is hard-coded and does not have an Application Programming Interface (API) package

84

to allow developers to connect it with external programs.43 Instead of modifying its source code,

5


Page 7 of 33


85

we were motivated to develop an external tool that can fit experimental data without altering the

86

widely-accepted user experience of MINEQL+.

87

The objective of this study was to develop a software program, which we are calling MINFIT,

88

that specifically interacts with MINEQL+ to optimize model fits to experimental datasets by

89

minimizing the residual sum of squares. Two algorithms, exhaustive equidistant grid search and

90

randomized parameter optimization, both of which do not rely on good initial guesses, were

91

implemented to ensure that the fit is at global rather than local minimum and is unique. We also

92

designed functions in MINFIT to enable straightforward sensitivity analysis of the individual

93

fitting parameter. A fundamental consideration of the program design was to maintain the widely

94

adopted user experience of MINEQL+ when integrating it with a transparent and fully user-

95

supervised optimization process.

96 97

98 99

Methods General Approach. The tableau algorithm of MINEQL+ to solve chemical equilibrium problems was comprehensively documented,44,

45

and its robustness and popularity have been

100

testified by thousands of applications in the literature. While MINEQL+ does not allow for

101

dynamic monitoring and control of the calculation routes, it fortunately has a multi-run calculation

102

option named “Field Data” analysis where the user can supply an external data file for processing. 6



Page 8 of 33

103

The “Field data” analyses were originally designed for chemical equilibrium simulations of field

104

samples with distinct chemistry conditions.3 Upon definition in MINEQL+, the external data file

105

can contain any experimental conditions to be simulated with variable equilibrium constants. This

106

option opens up the models to be truly data-driven, but has never been considered for parameter

107

optimization. MINFIT takes advantage of this function and tells MINEQL+ to perform a large

108

number of forward calculations within pre-defined ranges of parameter sets.

109

analyzes the calculation results by comparing them with the experimental results (Figure 1).

110

Although MINEQL+ requires inputs of all equilibrium constants, total concentration of each

111

component (or fixed concentration of certain species), temperature, solid concentration and

112

specific surface area, in its “Multi-run” option those input values can be varied for individual

113

calculations. Upon definition in MINEQL+, each line of the text, corresponding to a calculation,

114

provides the values of the necessary input parameters.

MINFIT then

115

Given the relatively light computation demand of most problems addressed by MINEQL+, the

116

efficiency of the optimization algorithm is not a limiting factor. Therefore, MINFIT adopted the

117

most straightforward trial-and-error approach and placed most of its effort in providing user-ease

118

and robustness. MINFIT was developed based on Microsoft Excel Visual Basic for Applications

119

(VBA), only requiring basic knowledge of spreadsheet calculations.

7


Page 9 of 33


120

Exhaustive Equidistant Grid Search. Solution of parameter estimation problem using

121

MINFIT starts with an exhaustive equidistant grid search that covers the widest ranges of

122

parameter sets. The number of fitting parameters (f), number of experiments (e), the number of

123

input parameters (i), and the number of datasets to fit in each experiment (d) are entered when

124

MINFIT initializes the spreadsheet. The detailed and illustrative explanations of those parameters

125

are presented in Figure S1 and Table S1. The experimental conditions (varies within Field Data

126

series as input parameters), the results to fit, and the corresponding weighting factor for each result

127

are tabulated in the “Problem” tab. With the upper (Un), lower (Ln) bounds, and the step lengths of

128

grid (Sn) of the nth fitting parameter, MINFIT automatically calculates (formula in nomenclature)

129

the number of grids for each dimension (Gn), the total number of parameter sets (p) and the total

130

number of MINEQL+ runs (r). MINFIT then generates the fitting parameter sets together with the

131

corresponding input conditions in a tabulated text file (Figure S1). A total number of d blank text

132

files are created to save the results that will be extracted from MINEQL+ after computation.

133

Randomized Parameter Search. MINFIT also offers a randomized search optimization

134

algorithm that generates random sets of parameter values to be screened. In this approach, only Un,

135

Ln and p are used to generate the “Field Data”, although Sn is still entered for each fitting

136

parameter to calculate Gn and subsequently p. The random parameter sets follow a normal

8



Page 10 of 33

137

distribution for each fitting parameter with a mean (µ) of (Un + Ln)/2 and a standard deviation (σ)

138

of (Un − Ln)/2. The selected standard deviation enables 31.8% of the randomized generated

139

parameters to fall larger than Un or smaller than Ln, allowing the search to dabble beyond the

140

“bounds” with a minor probability. When the search is narrowed down to smaller ranges, the

141

randomized search may or may not have a better chance to capture the minimum than the

142

exhaustive equidistant grid searches. We recommend randomized search as a verification step to

143

confirm that the RSS cannot be further minimized. Nevertheless, users have the flexibility to apply

144

hybrid strategies for their specific problems.

145

Data Extraction and Processing. MINEQL+ loads the “Field Data" file and runs the problem

146

under the “Field Data” option, where the meanings of each column of the text file are defined.

147

MINFIT is programed to require MINEQL+ to output the results for each species in a certain

148

format and save them as a two-column txt file (name and values, details in the Supporting

149

Information S2). The results are saved to the MINFIT-created blank text files in the order that is

150

consistent with the datasets. Subsequently, MINFIT reads the files and performed computation in

151

the “Calibration” tab. MINFIT allows users to freely define the objective function by manually

152

editing the formula that calculates RSS in the “Calibration” tab, if needed, in forms not readily

9


Page 11 of 33


153

reportable by MINEQL+ (e.g., the ratio of two concentrations). Screenshots of illustrative

154

examples are provided in the tutorial (T26).

155

MINFIT computes the residual sum of squares (RSS) from calculations for each parameter set

156

and compiles the optimal solutions in the “Summary” tab. A “Review Residual” function allows

157

users to selectively examine a certain range of simulation results and to visualize the quality of the

158

fitting (Supporting Information S2 and Tutorial T11). The searches can be repeated until the

159

parameter sets are narrowed down to a smaller range where a minimum RSS value may be located

160

(Figure 1).

161

Illustrative Examples. We use three examples to illustrate the use of MINFIT. The examples

162

cover a simple aqueous complexation problem and then metal adsorption using both the

163

generalized double layer and the more complex triple layer model. The objective functions include

164

concentration, fraction of adsorbed ion, and fraction of individual surface species. The values of f,

165

e and d also differ (Table S1). A step-by-step tutorial provides more details about those models

166

(Supporting Information). None of the examples involves multi-dentate surface complexation

167

reactions, whose mass action equations computed by MINEQL+ may not be valid if sorbent

168

concentration varies15, 46-48 (Supporting Information S1).

169 170

Results 10



Page 12 of 33

171

Cu complexation with pyromellitic acid. Giammar and Dzombak determined the formation

172

constants of copper complexes with mellitic acids using potentiometric titrations to generate

173

experimental data and FITEQL to identify equilibrium constants that provided the optimal fit.49.

174

The simple 1:1 complexes of CuL and CuHL (L denotes the fully deprotonated acid) could

175

describe all the data. We selected the data set of Cu complexation with pyromellitic acid to

176

illustrate MINFIT’s capability to reproduce their fitting (Table S2). All reactions were in the

177

aqueous phase and the logarithm concentrations were considered in the fitting. Figures 2a-2e plot

178

the experimental data and the optimized fit, after two exhaustive grid searches (Figures 2f and 2g)

179

and one randomized parameter search (Figure 2h) in a process of gradually narrowing down the

180

search ranges. There were noticeable differences between the values obtained by MINFIT and

181

FITEQL (Table S2). The optimal fits obtained by MINFIT (RSS = 1.11) were better than those by

182

FITEQL (RSS = 5.12). It was probably because the goodness of fit by FITEQL during the

183

iteration already reached the convergence criterion, so that FITEQL did not pursue further

184

refinements.

185

Chromate adsorption to goethite. Mathur and Dzombak compiled a generalized double layer

186

model dataset for the adsorption of various cations and anions to goethite.50 We selected chromate

187

adsorption to goethite as an example with three equilibrium constants as the fitting parameters.

11


Page 13 of 33


188

We reproduced Mathur and Dzombak’s model fit to the dataset of Mesuere and Fish51 with a total

189

chromate concentration of 0.2 mM and 1.8 g/L goethite at 0.1 M ionic strength (Table S3). We

190

benchmarked the ability of MINFIT to fit the data and inversely sought for the optimal set of three

191

equilibrium constants (Figure 3). MINFIT started with an exhaustive grid search with

192

incrementally narrowing down of the search ranges (Figures 3b-3e). Randomized parameter

193

optimization confirmed that the preceding grid search had already pinpointed the optimal solution

194

(Figure 3f). The obtained equilibrium constants and SOS/DF (sum of squares divided by degrees

195

of freedom) turned out to be very close to the original values obtained by Mathur (Table S3).52

196

Sulfate surface complexations on ferrihydrite. Using X-ray absorption near edge

197

spectroscopy (XANES), Gu et al. quantified the relative contribution and inner sphere and outer

198

sphere surface complexes of sulfate on ferrihydrite.53 They also developed an extended triple layer

199

model (ETLM) to fit experimental data. Compared with most inverse modeling problems in the

200

literature which only fit the macroscopic adsorption fractions, this study fit two datasets

201

simultaneously (d = 2). These datasets were the respective fractions of inner and outer sphere

202

surface complexes out of the total sulfate. MINFIT solved the optimization problem for the two

203

equilibrium constants (outer and inner sphere surface complexation reactions) by fitting the

204

fraction based surface speciation data as a function of pH at 0.02 M ionic strength and 1 mM total

12



Page 14 of 33

205

sulfate loading (Figure 4). Figures 4b-4d illustrate the process of narrowing down the ranges of

206

fitting parameters to the optimal values through two grid searches and one randomized search.

207

Sensitivity analysis of the fit to individual fitting parameter. The randomized search

208

screens a large number of fitting parameter sets surrounding the optimal solution, and the “Review

209

Residual” function of MINFIT sorts those results with residual from low to high. Users can

210

straightforwardly process and plot the data as in Figure 5, which indicates the sensitivity of the fit

211

to each fitting parameter. The range of log K is calculated from the maximum minus minimum of

212

the log K values as their corresponding residual increased from low to high (details in the tutorial

213

T16). Theoretically the curves should be smooth and continuous, but they are plotted from

214

randomly generated discretized data. For example, for the Cu-pyromellitic complexation model

215

(Figure 5a), departing from the optimal solution (residual 1.1108), those less optimal solutions (up

216

to residual 1.6) could allow both the log K for CuL2− and the log K for CuHL− to vary in a range

217

of 0.4, indicating that the overall goodness of fit was equally sensitive to the two equilibrium

218

constants. For the chromate-goethite model (Figure 5b), departing from the optimal solution

219

(residual 12.4), those less optimal solutions (up to residual equal to 14.2) could allow the log K for

220

≡FeCrO4− to vary in a range of 0.05, that for ≡FeHCrO4 and ≡FeOHCrO42− to vary by more than

221

0.1. It indicated that the goodness of fit was most sensitive to the log K for ≡FeCrO4−. This was

13


Page 15 of 33


222

intuitive as the ≡FeCrO4− species contributed predominantly to the most dramatic range of the

223

adsorption edge, and the other two species merely refined the fits where the edge was flat and not

224

as sensitive to pH (Figure 3a). For the sulfate-ferrihydrite model (Figure 5c), the fit was slightly

225

more sensitive to the log K for the outer sphere reaction. Equivalently, it could be interpreted that

226

varying the log K of the outer sphere reaction departing from the optimal values resulted in a more

227

drastic increase of residual. Alternative visualization methods for sensitivity analysis are described

228

in the tutorial (T26).

229 230

Discussions

231

Features and Advantages of MINFIT. (1) Free of Non-Convergence Problem. Compared

232

with other inverse modeling software programs, MINFIT is free of non-convergence problems. As

233

long as the forward model can be run on MINEQL+, then MINFIT can perform any parameter

234

estimation tasks within that model. It avoids the pitfalls related to problematic initial guesses, and

235

MINFIT is unlikely to return local minimum solutions as the equidistant grid search is exhaustive.

236

(2) Flexibility and Transparency of Data Organization. MINFIT can generate “Field Data”

237

files according to the user-defined experimental results in any format. Formats include adsorption

238

isotherm, potentiometric titration or even randomly organized datasets. The transparent structure

14



Page 16 of 33

239

of MINFIT also allows for the user to define any objective function that may be calculated from

240

metrics that MINEQL+ outputs. The users can modify the formula in calculating the individual

241

residual on the “Calibration” tab.

242

(3) Convenient Interface for Graphics. The “Review Residual” button in the “Calibration” tab

243

of MINFIT sorts those results and transfers selected results into the “Summary” tab using

244

hyperlink. This feature allows users to customize graphic visualization of the goodness of the fit as

245

well as the progress of optimization using built-in graphing functions of Excel that can be

246

automatically updated upon each click of the hyperlinks.

247

(4) Weighting Factor Table to Account for Experimental Error. The information generated by

248

MINFIT can calculate several useful statistical metrics that are reported by other software

249

programs. MINFIT allows the use of a weighing factor table that is assigned to each data point. A

250

well-known metric for the goodness of fit in FITEQL is WSOS/DF (weighted sum of

251

squares/degrees of freedom), which can be expressed as WSOS ∑ (measured − calcuated )2 = （1） DF ×−

252

where W is the weighing factor for the specific data point, Rmeasured and Mcalcuated are the

253

experimental and simulated data, respectively. FITEQL uses an estimate of the experimental error

254

associated with each experimental data point to calculate the weighting factors. If users choose to

15


Page 17 of 33


255

follow FITEQL’s WSOS/DF convention, then the value of W in MINFIT can be entered as equal

256

to !1/#, where s is the estimate of the experimental error. The weighting factors in MINFIT can

257

be supplied individually for each data point. Additional discussions about how to utilize W and

258

interpret WSOS/DF are presented in the Supporting Information (S4).

259

(5) Sensitivity Analysis. MINFIT enables convenient evaluation of sensitivity of the objective

260

function in response to the variation of the individual fitting parameter. When the search narrows

261

down to the ranges where the optimal solution resides, the “Review Residual” function can

262

compile and sort parameter sets that return residuals smaller than the user-provided threshold

263

value. Examining the ranges of those fitting parameters as the residual increases from the lowest

264

to highest gives information about how sensitive the fit is to each fitting parameter (Figure 5).

265

Although MINFIT does not compute confidence intervals or standard deviations, the sensitivity

266

analysis can provide useful information indicating the relative uncertainty of the optimal

267

parameters.

268

Exhaustive Grid Search vs Randomized Parameter Optimization. Exhaustive equidistant

269

grid search is recommended to conduct the initial screening. The known and fixed grid length

270

allow users to refine the next grid search. Randomized parameter search is recommended to verify

271

the optimal solution obtained from equidistant grid search, because it can generate a large number

16



Page 18 of 33

272

of random parameter sets surrounding the values to be verified. The minimum residual sum of

273

squares would be easily visualized in the built-in plot in the “Calibration” tab (i.e., the “wells” in

274

exhaustive grid search and the “bound” in randomized search). Nevertheless, before concluding

275

the optimization, it is always advisable to use the graphic functions of MINFIT to visually

276

examine the goodness of the fits.54

277

Limitations of MINFIT. An obvious limitation of MINFIT is that it is specifically developed

278

to interface with MINEQL+, although the potential user community is already quite substantial. In

279

addition, we summarize several other limitations related to the algorithms and the MINEQL+

280

working platform, and discuss possible alternatives to circumvent those limitations.

281

(1) Computation Demands. Our major objective is to develop an easy-to-use tool to enable

282

parameter estimation. Therefore, the novelty or efficiency of the algorithm is not our focus as long

283

as the learning curve to adopt MINFIT is reduced. In fact, a personal laptop (8GB RAM,

284

CPU1.8GHz) can finish a computation with 8,000 runs of a triple layer surface complexation

285

model in about three minutes and a non-electrostatic aqueous speciation model in about one

286

minute. Data extraction from MINEQL+ and the data processing in MINFIT can be performed

287

efficiently.

17


Page 19 of 33


288

(2) Column Number Limitation. As MINEQL+ limits the number of columns of the “Field

289

Data” below 25, MINFIT requires the sum of i and f to be smaller than 25. Nevertheless, most of

290

the parameter estimation problems in the related literature or applications would not exceed this

291

limit. For input parameters that are not varied throughout the rows of the “Field Data”, they can be

292

entered in the “Calculation Wizard” of MINEQL+, reducing the number of columns taken in the

293

text file.

294

(3) Output Manager Limitation. We identified a bug in MINEQL+ 4.6 that when the total line

295

number of the “Field Data” is larger than a number around 8000, the “Output Manager” crashes

296

when displaying results even though all the calculations can be finished. This becomes a limitation

297

for MINFIT when the number of experiments (e) is large, which means that users should carefully

298

allocate the density of the grids on each dimension. From our experiences, each round of

299

parameter search with 8000 runs would satisfy most of our needs. In the most recently released

300

version of MINEQL+ (5.0), this bug is fixed and we tested that MINEQL+ 5.0’s Output Manager

301

can handle up to 100,000 runs in each “Field Data” processing, although the computing time and

302

output file size also increase dramatically. As the problem complexity varies, MINFIT would let

303

the users decide on the tradeoff between the total number of searches needed and the computing

304

time needed for each search.

18



Page 20 of 33

305

(4) Fitting Parameter Limitation. MINFIT allows up to three independent fitting parameters to

306

be optimized simultaneously. A larger number of fitting parameters usually would give further

307

refinement of the fits, but they also run the risk of over-fitting and getting non-unique solutions. In

308

most cases, a robust model would seldom have too many parameters to be determined at a time

309

unless justified by specific physical/chemical constraints. Instead, a bottom-up approach is

310

recommended to build models from simpler to more complex situations, which essentially reduces

311

the degree of freedoms by dividing the problem into several simpler sub-problems.12,

312

Therefore, MINFIT would be sufficient to address most needs of the user community. If more

313

fitting parameters must be optimized simultaneously, then other advanced software programs (that

314

cannot be integrated with MINEQL+) such as PEST31 and UCODE32 can still be used. The types

315

of inputs allowed in “Field Data” are sufficient for most chemical equilibrium problems, but the

316

most recent MINEQL+ 5.0 still does not include capacitances. More detailed discussions about the

317

fitting parameter limitations are presented in the Supporting Information (S5).

24, 55-60

318

Environmental Implications The availability of a user-friendly parameter estimation tool for

319

use in connection with MINEQL+ can facilitate the performance of parameter estimation by more

320

environmental science and engineering users. We anticipate that the applications of FITEQL,

321

ECOSAT-FIT, ProtoFit, PEST, UCODES will persist in various research communities, but

19


Page 21 of 33


322

MINFIT offers a valuable tool with a relatively easy learning curve and remarkable robustness.

323

MINFIT enables broader applications of a systematic approach for a range of parameter

324

optimization problem types. Parameter estimation can be used in a broad range of studies

325

requiring chemical equilibrium calculations for systems with incompletely known equilibrium

326

constants. In addition to the complexation and adsorption examples presented here, MINFIT could

327

be used for problems involving metal precipitate solubility and acid-base problems.

328

A shallower barrier to parameter estimation also means that more published datasets that could

329

be used to yield equilibrium constants will actually yield such values. Presenting not just data but

330

reaction parameters that can simulate the data makes the results of a study easier to be compared

331

with or used in other research. For example, chemical speciation based reactive transport models

332

should be built on robust datasets of the thermodynamic constants for the critical reactions, and

333

there are increasing amounts of interesting datasets from which reaction parameters can be

334

estimated for use in those reactive transport models.

335

Acknowledgements

336

The work was supported by a new faculty start-up award from Louisiana State University and a

337

Chinese Government Award for Outstanding Self-Financed Students Abroad. We acknowledge

338

the technical support and encouragement from MINEQL+ developer Dr. Bill Schecher

20



Page 22 of 33

339

(Environmental Research Software). Prof. Dave Dzombak provided the original data source of

340

chromate adsorption that they used to build the goethite surface complexation database.

341

Comments and suggestions of Associate Editor Timm Starthmann and three anonymous reviewers

342

improved the clarify of the manuscript.

343

Supporting Information

344

The supporting information (extended discussion, additional information about the three

345

illustrative examples, the links to download the MINEQL+ input files and a step-by-step user

346

tutorial) is available free of charge on the ACS Publications website.

347 348

Nomenclature

349

Entered by users

350

f

number of fitting parameters, f = 1, 2, 3

351

i

number of input parameters, i + f ≤ 25

352

d

number of datasets to fit in each experiment, i + f ≤ 25

353

e

number of experiments

354

Ln

the lower bound of the nth fitting parameter

355

Un

the upper bound of the nth fitting parameter

356

Sn

Step length of each grid

357

W

weighing factor for the specific data point to fit defined in MINFIT, W = 1 by default

358

Rmeasured experimental data

359

Mcalculated simulated data

360

Calculated by MINFIT

361

Gn

number of grids, int((Un− Ln)/Sn+1)

362

p

total number of parameter sets, ∏&() %&

363

r

total number of MINEQL+ runs, × ∏&() %& . e ≤ 8,000 for V4.6, e ≤ 100,000 for V5.0.

'

'

21


Page 23 of 33


364

µ

mean of the normal distributed fitting parameter in randomized parameter search, (Un + Ln)/2

365

σ

standard deviation the normal distributed fitting parameter in randomized parameter search, (Un − Ln)/2

366

RSS

residual sum of squares

367

WSOS weighted sum of squares, ∑( (*+,-./+0 − 1,23.,4+0 )5 )

368

SOS

369

Others

370

DF

degree of freedom, × −

371

s

experimental error

372

References

sum of squares, ∑((*+,-./+0 − 1,23.,4+0 )5 )

373 374

1.

Miller, A. W.; Rodriguez, D. R.; Honeyman, B. D., Upscaling sorption/desorption processes in

375

reactive transport models to describe metal/radionuclide transport: A critical review. Environ. Sci. Technol.

376

2010, 44, (21), 7996-8007.

377

2.

378

Technol. 2012, 46, (11), 5650-5657.

379

3.

380

Environmental Research Software: Hallowell, ME, 2007.

381

4.

382

montmorillonite Part II: modelling. J. Contam. Hydrol. 1997, 27, (3), 223-248.

383

5.

384

equilibrium constants from experimental data, Oregon State University: Corvallis, OR, 1999.

385

6.

Gustafsson, J. P. Visual MINTEQ 3.0 user guide; Royal Institute of Technology: 2011.

386

7.

Parkhurst, D. L.; Appelo, C. A. J., Description of input and examples for PHREEQC version 3—A

387

computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical

388

calculations. In U.S. Geological Survey Techniques and Methods, Book 6, Modeling Techniques, Chapter

389

43, 2013.

390

8.

391

calculation of speciation and transport in soil-water systems, Department of Environmental Sciences, Soil

392

Quality, Wageninggen University: Wageningen, The Netherlands, 2002.

393

9.

394

equilibrium models. Environ. Sci. Technol. 2003, 37, (6), 1175-1182.

395

10.

396

surfaces from aqueous solution. Comput. Geosci. 1998, 24, (9), 853-873.

397

11.

398

Hydrogeology Program, University of Illinois: 2010.

Brezonik, P. L.; Arnold, W. A., Water chemistry: Fifty years of change and progress. Environ. Sci. Schecher, W. D.; McAvoy, D. C. MINEQL+: A chemical equilibrium modeling system, version 4.6, Bradbury, M. H.; Baeyens, B., A mechanistic description of Ni and Zn sorption on NaHerbelin, A. L.; Westall, J. C. FITEQL 4.0: A computer program for the determination of chemical

Keizer, M. G.; Van Riemsdijk, W. H. ECOSAT version 4.7. a computer program for the

Meeussen, J. C. L., ORCHESTRA: An object-oriented framework for implementing chemical Sahai, N.; Sverjensky, D. A., GEOSURF: A computer program for modeling adsorption on mineral Bethke, C. M.; Yeakel, S. The Geochemist's Workbench, Release 8.0 GWB Essentials Guide, 22



Page 24 of 33

399

12.

Payne, T. E.; Brendler, V.; Ochs, M.; Baeyens, B.; Brown, P. L.; Davis, J. A.; Ekberg, C.; Kulik, D.

400

A.; Lutzenkirchen, J.; Missana, T.; Tachi, Y.; Van Loon, L. R.; Altmann, S., Guidelines for thermodynamic

401

sorption modelling in the context of radioactive waste disposal. Environ. Modell. Softw. 2013, 42, 143-156.

402

13.

403

Geochem. 1990, 23, (1), 177-260.

404

14.

405

geochemist's perspective. J. Hydrol. 2000, 230, (3), 127-171.

406

15.

407

complexation modeling: Theory and practice. Environ. Sci. Technol. 2013, 47, (9), 3982-3996.

408

16.

409

model parameter estimation using FITEQL and oxide mineral titration data. J. Colloid Interface Sci. 1991,

410

142, (2), 448-469.

411

17.

412

phosphate adsorption on γ-alumina and kaolinite: triple-layer model. Soil Sci. Soc. Am. J. 1997, 61, (3),

413

784-793.

414

18.

415

Sons: 2011.

416

19.

417

E., Nanoscale size effects on uranium(VI) adsorption to hematite. Environ. Sci. Technol. 2009, 43, (5),

418

1373-1378.

419

20.

420

a comparison of some humic – metal complexation models. Chem. Speciation Bioavailability 2003, 15, (4),

421

93-100.

422

21.

423

complexation modeling of proton and metal sorption onto graphene oxide. Colloids Surf., A 2015, 466, 28-

424

39.

425

22.

426

montmorillonites. II. Surface complexation modeling. Environ. Sci. Technol. 2014, 48, (15), 8698-8705.

427

23.

428

and head groups of cationic surfactants on retention of cesium by organoclays. Environ. Sci. Technol. 2010,

429

44, (13), 5142-5147.

430

24.

431

D. E., Adsorption of uranium(VI) to manganese oxides: X-ray absorption spectroscopy and surface

432

complexation modeling. Environ. Sci. Technol. 2013, 47, (2), 850-858.

Davis, J. A.; Kent, D. B., Surface complexation modeling in aqueous geochemistry. Rev. Mineral. Koretsky, C., The significance of surface complexation reactions in hydrologic systems: a Wang, Z.; Giammar, D. E., Mass action expressions for bidentate adsorption in surface Hayes, K. F.; Redden, G.; Ela, W.; Leckie, J. O., Surface complexation models: An evaluation of

He, L.; Zelazny, L.; Martens, D.; Baligar, V.; Ritchey, K., Ionic strength effects on sulfate and

Karamalidis, A. K.; Dzombak, D. A., Surface complexation modeling: gibbsite. John Wiley & Zeng, H.; Singh, A.; Basak, S.; Ulrich, K.-U.; Sahu, M.; Biswas, P.; Catalano, J. G.; Giammar, D.

Evans, L. J.; Sengdy, B.; Lumsdon, D. G.; Stanbury, D. A., Cadmium adsorption by an organic soil:

Duster, T. A.; Szymanowski, J. E. S.; Na, C.; Showalter, A. R.; Bunker, B. A.; Fein, J. B., Surface

Soltermann, D.; Baeyens, B.; Bradbury, M. H.; Fernandes, M. M., Fe (II) uptake on natural Wang, T.-H.; Hsieh, C.-J.; Lin, S.-M.; Wu, D.-C.; Li, M.-H.; Teng, S.-P., Effect of alkyl properties

Wang, Z.; Lee, S.-W.; Catalano, J. G.; Lezama-Pacheco, J. S.; Bargar, J. R.; Tebo, B. M.; Giammar,

23


Page 25 of 33


433

25.

Bradbury, M. H.; Baeyens, B., Modelling the sorption of Mn (II), Co (II), Ni (II), Zn (II), Cd (II),

434

Eu (III), Am (III), Sn (IV), Th (IV), Np (V) and U (VI) on montmorillonite: Linear free energy

435

relationships and estimates of surface binding constants for some selected heavy metals and actinides.

436

Geochim. Cosmochim. Acta 2005, 69, (4), 875-892.

437

26.

438

to bind heavy metal ions: Based on formation environment, adsorption reversibility and ageing. Appl.

439

Geochem. 2014, 45, 114-119.

440

27.

441

reductive dissolution: Role of lead(II) adsorption and surface speciation. J. Colloid Interface Sci. 2013, 389,

442

(1), 236-243.

443

28.

444

Keyworth: 1993.

445

29.

446

titration data. Comput. Geosci. 2006, 32, (9), 1344-1356.

447

30.

448

species using a genetic algorithm. Math. Geosci. 2010, 42, (1), 101-127.

449

31.

450

Computing, Corinda, Australia 1994, 122.

451

32.

452

computer codes for universal sensitivity analysis, calibration, and uncertainty evaluation constructed using

453

the JUPITER API; 2328-7055; 2006.

454

33.

455

PEST with PHREEQC. Freiberg online geology (FOG), Freiberg 2011, 28.

456

34.

457

modelling study. Chem. Geol. 2015, 417, 279-288.

458

35.

459

and its application to a column experiment. Geochim. Cosmochim. Acta 1999, 63, (19-20), 3039-3048.

460

36.

461

N.; Hatcher, P. G.; Santschi, P. H., Plutonium immobilization and remobilization by soil mineral and

462

organic matter in the far-field of the Savannah River Site, US. Environ. Sci. Technol. 2014, 48, (6), 3186-

463

3195.

464

37.

465

spectroscopic identification of neptunium(V) inner-sphere complexes on the hematite–water interface.

466

Environ. Sci. Technol. 2015, 49, (4), 2560-2567.

Meng, S.; Wang, H.; Liu, H.; Yang, C.; Wei, Y.; Hou, D., Evaluation of the ability of ferrihydrite

Wang, Y.; Wu, J.; Wang, Z.; Terenyi, A.; Giammar, D. E., Kinetics of lead(IV) oxide (PbO2)

Kinniburgh, D. FIT User Guide; BGS Technical Report WD/93/23, British Geological Survey, Turner, B. F.; Fein, J. B., Protofit: a program for determining surface protonation constants from Villegas-Jiménez, A.; Mucci, A., Estimating intrinsic formation constants of mineral surface Doherty, J.; Brebber, L.; Whyte, P., PEST: Model-independent parameter estimation. Watermark Poeter, E. E.; Hill, M. C.; Banta, E. R.; Mehl, S.; Christensen, S. UCODE_2005 and six other

Bachmaf, S.; Merkel, B. J., Estimating water chemistry parameters from experimental data using Gustafsson, J. P.; Tiberg, C., Molybdenum binding to soil constituents in acid soils: An XAS and Appelo, C. A. J.; Postma, D., A consistent model for surface complexation on birnessite (δ-MnO2) Xu, C.; Athon, M.; Ho, Y.-F.; Chang, H.-S.; Zhang, S.; Kaplan, D. I.; Schwehr, K. A.; DiDonato,

Müller, K.; Gröschel, A.; Rossberg, A.; Bok, F.; Franzen, C.; Brendler, V.; Foerstendorf, H., In situ

24



Page 26 of 33

467

38.

Fakhraei, H.; Driscoll, C. T., Proton and aluminum binding properties of organic acids in surface

468

waters of the northeastern US. Environ. Sci. Technol. 2015, 49, (5), 2939-2947.

469

39.

470

optimization of geochemical model parameters and inverse modeling. Appl. Geochem. 2015, 55, 28-45.

471

40.

472

tool for geochemical activity models. Computat. Geosci. 2014, 18, (2), 227-242.

473

41.

Kinniburgh, D.; Cooper, D., PhreePlot: Creating graphical output with PHREEQC. 2011.

474

42.

Mahoney, J., Using PhreePlot to calibrate mining related geochemical models: A user’s perspective.

475

In 10th International Conference on Acid Rock Drainage & IMWA Annual Conference,

476

https://www.imwa.info/docs/imwa_2015/IMWA2015_Mahoney_193.pdf, 2015.

477

43.

Which Model is For You? http://www.mineql.com/BigView6.html (2016/02/05),

478

44.

Benjamin, M. M., Water Chemistry. 1st ed.; McGraw-Hill: New York, 2002.

479

45.

Brezonik, P. L.; Arnold, W. A., Water Chemistry: An Introduction to the Chemistry of Natural and

480

Engineered Aquatic Systems. Oxford University Press: New York, 2011.

481

46.

482

Treatment of multi-dentate surface complexes and diffuse layer implementation in various speciation codes.

483

Appl. Geochem. 2015, 55, 128-137.

484

47.

485

phosphate concentration in soils: Results of theoretical modelling” by Marek Duputel, Nicolas Devau,

486

Michel Brossard, Benoît Jaillard, Davey L. Jones, Philippe Hinsinger and Frédéric Gérard (2013). Appl.

487

Geochem. 2014, 46, (0), 85-89.

488

48.

489

Powell, B. A., A self-consistent model describing the thermodynamics of Eu(III) adsorption onto hematite.

490

Geochim. Cosmochim. Acta 2013, 122, 430-447.

491

49.

492

Chem. 1998, 27, (1), 89-105.

493

50.

494

Science and Technology, Lützenkirchen, J., Ed. Elsevier: 2006; Vol. Volume 11, pp 443-468.

495

51.

496

complexation models. Environ. Sci. Technol. 1992, 26, (12), 2357-2364.

497

52.

498

modeling. Carnegie Mellon University, Pittsburgh, 1995.

499

53.

500

quantification and speciation modeling of sulfate adsorption on ferrihydrite surfaces. Environ. Sci. Technol.

501

2016, 50, (15), 8067-8076.

Miron, G. D.; Kulik, D. A.; Dmytrieva, S. V.; Wagner, T., GEMSFITS: Code package for Hingerl, F. F.; Kosakowski, G.; Wagner, T.; Kulik, D. A.; Driesner, T., GEMSFIT: a generic fitting

Lützenkirchen, J.; Marsac, R.; Kulik, D. A.; Payne, T. E.; Xue, Z.; Orsetti, S.; Haderlein, S. B.,

Gustafsson, J. P.; Lumsdon, D. G., Comment on “Citrate adsorption can decrease soluble

Estes, S. L.; Arai, Y.; Becker, U.; Fernando, S.; Yuan, K.; Ewing, R. C.; Zhang, J.; Shibata, T.;

Giammar, D. E.; Dzombak, D. A., Copper complexation with the mellitic acid series. J. Solution Mathur, S. S.; Dzombak, D. A., Chapter 16 Surface complexation modeling: goethite. In Interface Mesuere, K.; Fish, W., Chromate and oxalate adsorption on goethite. 1. Calibration of surface Mathur, S. S. Development of a database for ion sorption on goethite using surface complexation Gu, C.; Wang, Z.; Kubicki, J. D.; Wang, X.; Zhu, M., X-ray absorption spectroscopic

25


Page 27 of 33


502

54.

Lützenkirchen, J., Comparison of 1-pK and 2-pK versions of surface complexation theory by the

503

goodness of fit in describing surface charge data of (hydr)oxides. Environ. Sci. Technol. 1998, 32, (20),

504

3149-3154.

505

55.

506

to metal oxide minerals. Environ. Sci. Technol. Lett. 2015, 2, (8), 227-232.

507

56.

508

quartz-sand dominated sediments. Environ. Sci. Technol. 2014, 48, (12), 6569-6577.

509

57.

510

capabilities of diffuse double layer model. Geochim. Cosmochim. Acta 2005, 69, (11), 2725-2734.

511

58.

512

and sulfate on iron oxide surfaces. Geochim. Cosmochim. Acta 2015, 158, 130-146.

513

59.

514

Comparison of approaches to calibrate a surface complexation model for U(VI) sorption to weathered

515

saprolite. Transport Porous Media 2009, 78, (2), 185-197.

516

60.

517

Effect of phosphate on U (VI) sorption to montmorillonite: Ternary complexation and precipitation barriers.

518

Geochim. Cosmochim. Acta 2016, 175, 86-99.

Wang, Z.; Ulrich, K.-U.; Pan, C.; Giammar, D. E., Measurement and modeling of U(IV) adsorption Dong, W.; Wan, J., Additive surface complexation modeling of uranium(VI) adsorption onto Richter, A.; Brendler, V.; Nebelung, C., Blind prediction of Cu(II) sorption onto goethite: Current Hinkle, M. A. G.; Wang, Z.; Giammar, D. E.; Catalano, J. G., Interaction of Fe(II) with phosphate Zhang, F.; Parker, J. C.; Brooks, S. C.; Kim, Y.-J.; Tang, G.; Jardine, P. M.; Watson, D. B.,

Troyer, L. D.; Maillot, F.; Wang, Z.; Wang, Z.; Mehta, V. S.; Giammar, D. E.; Catalano, J. G.,

519

26



520 521

Figure 1. Work flow of MINFIT in support of MINEQL+ to perform parameter estimation in chemical

522

equilibrium models.

27


Page 28 of 33

Page 29 of 33


523 524

Figure 2. Illustration of using MINFIT to determine the equilibrium constants for Cu complexation with

525

pyromellitic acid at a total concentration of 10 mM. Ionic strength is 0.1 M. Symbols and lines indicate

526

data and simulation respectively (a-e). With gradually refining search ranges, Panels f-h were generated by

527

MINFIT that plot the RSS values for each parameter set that was screened. The number of parameter set is

528

a unique MINFIT-assigned identifier for the screened log K or a set of log Ks (details in Figure S1 and

529

tutorial). The inset text presents the search parameters of exhaustive grid search (f and g) and randomized

530

search (h) and the values of the log K and RSS that give the optimal fit in the specific round of search

531

(details in Table S2). Note that Panels f-h are directly extracted from the MINFIT interface for illustrative

532

purposes. The bounds and units of both axes can be edited manually.

28



Page 30 of 33

533 534

Figure 3. Illustration of using MINFIT to determine the equilibrium constants for chromate surface

535

complexation on goethite within the generalized double layer model. Ionic strength is 0.1 M. The data were

536

from Mesuere and Fish51, and the basic model parameters (surface (de)protonation, goethite site

537

density/concentration) were from Mathur and Dzombak50. Symbols and lines indicate data and simulation

538

respectively (a). The inset text presents the search parameters of exhaustive grid search (b to e) and

539

randomized search (f). The values of the log K and RSS that give the optimal fit in the specific round of

540

search, as well as how those results are compared with those obtained by FITEQL, are presented in Table

541

S3.

29


Page 31 of 33


542 543

Figure 4. Illustration of using MINFIT to determine the equilibrium constants for sulfate surface

544

complexation on ferrihydrite within the extended triple layer model. Symbols and lines indicate data and

545

simulation respectively. Ionic strength is 0.02 M. The fractions of inner and outer sphere surface complexes

546

of sulfate out of the total sulfate, as determined by XANES, were considered in the model fitting. The total

547

adsorbed sulfate data and simulation were included for visual references, but were not considered in the

548

optimization calculation. Symbols and lines indicate data and simulation respectively (a). The inset text

549

presents the search parameters of exhaustive grid search (b and c) and randomized search (d) and the values

550

of the log K and RSS that give the optimal fit in the specific round of search. The detailed model

551

formulation was reported in Gu et al.53 and summarized in Table S4.

30



Page 32 of 33

552 16 p1CuL4− p2CuHL3−

15

1.5

≡FeOH2+--HSO4− p1

1600

1.3 13

(a) 0

0.1 0.2 0.3 Range of Log K

0.4

(≡FeO)2SO2 p2

1700

14

1.1

553

1800 P1 ≡FeCrO4− P2 ≡FeHCrO4 ≡FeOHCrO42− P3

Residual

Residual Sum of Squares

1.7

1500

(b)

12 0

0.05 0.1 Range of Log K

0.15

(c)

1400 0

0.2 0.4 0.6 0.8 Range of Log K

1

554

Figure 5. Sensitivity analysis for the problems of (a) Cu complexation with pyromellitic acid, (b) double

555

layer model of chromate adsorption to goethite and (c) triple layer model of sulfate adsorption to

556

ferrihydrite. The plots were prepared with the simulation results from the randomized parameter search (i.e.,

557

the data of Figures 2h, 3f and 4d) surrounding the optimal solutions obtained from the precedent grid

558

search screenings. The range of log K indicates the deviation away from the optimal log K. The range of

559

log K was calculated from the maximum minus minimum of the log K values as their corresponding

560

residual increased from low to high. The formula to calculate them in the spreadsheet is presented in the

561

tutorial (T16).

562

31


Page 33 of 33

563


Table of Content Art

564

32


MINFIT: A Spreadsheet-Based Tool for Parameter Estimation in an

Recommend Documents