Geochemical Rate Models
Geochemical Rate Models
Geochemical Rate Models
J. Donald Rimstidt
Published in the United States of America by Cambridge University Press, New York
www.cambridge.org
Information on this title: www.cambridge.org/9781107029972
© J. Donald Rimstidt 2014
First published 2014
Printed and bound in the United Kingdom by TJ International Ltd., Padstow, Cornwall
A catalogue record for this publication is available from the British Library
Preface page ix
1 Geochemical models 1
2 Modeling tools 9
3 Rate equations 36
4 Chemical reactors 56
5 Molecular kinetics 79
6 Surface kinetics 102
7 Diffusion and advection 128
8 Quasi-kinetics 156
9 Accretion and transformation kinetics 182
10 Pattern formation 205
References 210
Index 228
vii
Preface
ix
x Preface
with actual outcomes to test the validity of the underlying knowledge. Every
scientist should understand and appreciate the importance of this model-
building and -testing scenario.
This book does not delve into computer models of rate processes because
there are already many good sources of information about this approach.
Instead it focuses on fundamental and relatively simple models that are the
foundations of the models incorporated into the various computer codes.
I have enjoyed the process of writing this book and learned a great deal
from teaching this material to students. I hope that the readers will find it
equally interesting and enlightening.
Chapter 1
Geochemical models
Modern geochemistry studies the distribution and amounts of the chemical elem-
ents in minerals, ores, rocks, soils, waters, and the atmosphere, and the circulation
of the elements in nature, on the basis of the properties of their atoms and ions.
(Goldschmidt, 1958)
The distribution and circulation of the chemical elements in and on the Earth
is influenced by a myriad of chemical and physical factors, many of which
have changed over geological time. Understanding the role of these factors
in geological processes requires us to condense information about elemental
abundances and distributions into models. This book is about geochemical
models for situations where time plays a key role. Geoscientists have always
appreciated the importance of time in fashioning the Earth. Many geological
processes require time spans that are far too long for human observation,
but we can use models to extrapolate rates based on short-term observations
to predict geochemistry in deep time. Equally important are models that
forecast the future behavior of geochemical systems because those models
are needed for environmental management and resources recovery projects.
Some of the models described in this book were developed by geochemists
but many others come from applied sciences and engineering. Because of
this diverse provenance, the models in their original form used a confusing
mix of units, terminology, and notation. This book attempts to remedy that
problem by recasting the models using internally consistent notation, units,
and terminology familiar to geochemists. Furthermore, whenever possible
the models are developed from fundamental theory showing a sufficient
number of intermediate steps to allow the reader to follow the derivations.
Thermodynamic models have been a mainstay of geochemistry since the
early twentieth century. They are especially effective for deep earth condi-
tions where local equilibrium conditions prevail. However, at and near the
Earth’s surface extensive amounts of mass transport and low temperatures
keep many reactions from reaching equilibrium. Kinetics models are needed
to properly describe these situations. This makes the models of kinetic and
dynamic processes described in this book complementary to thermodynamic
1
2 Geochemical models
Model construction
Models organize our knowledge about the world. Doing science is the pro-
cess of making observations, using those observations to develop a model,
and then verifying the model’s effectiveness by comparing its predictions
with additional observations. The methods and scope of model-making
vary from discipline to discipline, but the goal of creating reliable predictive
tools is always the same. Developing a useful model requires a combination
of creative thinking and disciplined use of modeling tools. There are many
good discussions of modeling strategies (Aris, 1994; Bender, 1978; Bunge,
1997; Hall and Day, 1977; Hall et al., 1977; Harte, 1988; Overton, 1977;
Rescigno and Thakur, 1988; Shoemaker, 1977) but developing an effective
model is in many ways a creative process. Good models are elegant and
powerful. Elegance means that the model is expressed in the simplest, easi-
est to understand terms. Power means that the model explains the widest
Model construction 3
Model reliability
No model is perfect (Oreskes et al., 1994). Even the best models fail under
some conditions. In addition, making mistakes is a natural part of the mod-
eling process. Systematic methods should be used to find, analyze, and cor-
rect mistakes and to define the valid range of a model. Determining the cause
of a model’s failure and repairing the model is a key task in model building.
Models are used to understand how complex behaviors arise from the
interaction of simple processes. Ideally models are built upon reliable prin-
ciples such as the conservation of matter, energy, and charge or the prin-
ciple of detailed balancing. Often less-reliable relationships such as empir-
ical rate equations must also be used. All quantitative models require input
data and relationships that come from measurements, which always contain
some error. This error can propagate through a model in unexpected ways,
especially if the model simulates nonlinear interactions. This challenges the
modeler who must decide whether an unexpected result is a legitimate pre-
diction or simply an artifact arising from an unfortunate combination of
errors. Three classes of errors occur in models. Formal errors are incorrect
assumptions and/or formulations. They include errors in the conceptual
foundation of the model as well as errors in the input data. Structural errors
are errors in mathematical manipulations such as programming errors or
algebraic errors. They include software bugs. Computational errors are
errors in numbers caused by incorrect rounding or by addition or multi-
plication errors. Because there are so many ways that a model might fail,
models and their predictions must be verified and validated to delineate the
bounds of their reliability.
Verification tests whether the model is internally consistent, incorpo-
rates the correct relationships in the correct ways, and uses correct data.
It is a good idea to develop a set of standards and practices that can be
used in geochemistry model-making to insure the validity of the resulting
model. Table 1.1 is an example of a checklist that might be used. Models
should use a consistent system of notation and units to avoid structural
errors. Equations should be tested using dimensional analysis. The input
data should be reviewed to insure its correctness and internal consistency.
Verification should determine the expected precision of the model’s predic-
tions based on the propagation of errors through the model. Error ana-
lysis is probably the most undervalued part of model development. For
simple models, error propagation can be done using simple algebraic meth-
ods described in Chapter 2. For complex models, simple error propagation
is often not practical. In such cases sensitivity analysis is a good strategy.
Sensitivity analysis uses various schemes to systematically vary, within
the expected range of error, the values of important modeling parameters
to see how those variations affect the model’s predictions. The utility of
Model reliability 5
Verification Questions
Validation Questions
Are the predictions geologically reasonable? Are they consistent with reasonable
estimates?
Are the predictions consistent with other scientific observations and knowledge?
Are the quantitative predictions sufficiently close to the behavior of one or more
natural analogs?
into the water table, the concentration of dissolved species increases to perhaps
500 mg/L. The total mass of rock and soil dissolved away each year is 5 × 104 mg
or ~50 g. If the density of the rock and soil minerals is ~3 g/cm3, then 17 cm3
(= 1.7 × 10–5 m3) of material is carried away annually. If that volume is removed
from under the 1 m2 of land surface, the surface elevation will be lowered by
1.7 × 10–5 m (~0.02 mm). So the land surface elevation is reduced by
0.02 mm/yr or ~2 mm/century.
These estimates show that chemical weathering denudation rates should be
on the order of a few millimeters per century. They might be as high as a few
centimeters per century in humid tropical settings or as low as a few hundred
microns per century in arid cold environments. Because these estimates con-
strain model predictions to fall within this range of values, they can be used as a
validation test of a chemical weathering denudation model.
Harte (1988) and Weinstein and Adam (2008) provide detailed explanations
about how to make meaningful estimations and they give a large number of
practice examples.
Interpretation of results
The results of simple deterministic models are generally easy to interpret
and apply. For example, a model of radioactive decay will predict that the
radioactivity of a substance will be 1/8 of the original value after three
half-lives have passed. As models become more complex they become more
difficult to interpret. This suggests that model building should begin with
highly idealized and simplified cases that are easily understood and tested.
The next stage of model building involves testing the effect of potentially
important variables on the model output to determine which variables are
important enough to include and which can be neglected. This process of
building models in a stepwise fashion not only eliminates a large number of
unneeded independent variables but it helps the modeler develop a concep-
tual understanding of the situation being modeled.
The most important task of model building is to communicate expert
knowledge developed by the model builder to others who can use that
knowledge to solve problems. The modeling endpoint should be a report
that contains thorough documentation of the model and an explanation
of how to use its predictions. This report should contain a clear descrip-
tion of the model’s conceptual basis in a well-written narrative section with
accompanying illustrations and explanations of the mathematical methods.
The report should illustrate relationships between the model’s predictions
and input parameters using response maps. These graphs of predictions ver-
sus parameters are much easier to understand than equations and tables.
8 Geochemical models
Many of the figures in this book are response maps. The report should also
describe the implementation of the modeling algorithm, including descrip-
tions of special computational methods. The rationale for selecting the
input data, including a discussion of their uncertainty, should be described.
The range and domain of the model should be specified. Finally, the report
should explain how to interpret the model’s output.
Chapter 2
Modeling tools
Before any model is ready for use it must be verified. The verification step
is greatly simplified if the model is constructed using conventional com-
putational methods, notation, and units. This chapter reviews some pro-
cedures and conventions that are recommended for geochemical model
construction.
9
10 Modeling tools
1. Write the chemical formulae for the reactant and product species on opposite
sides of an equal sign.
2. Find the element that is least represented in these species and insert integer (or
sometimes small fraction) coefficients into the equation to make equivalent
amounts of this element on each side of the equation.
3. Repeat this step with the second least represented element and so on until all the
elements are balanced. For reactions involving aqueous solutions, it is usually
necessary to add H2O, O2, H2, or H+ to balance the hydrogen and oxygen.
4. Adjust the coefficients of the charged species to produce an equivalent amount
of charges on each side of the equation.
5. Finally, divide the entire equation through by integers to reduce these coefficients
to small integer values (although sometimes small fractions are acceptable).
1. We know that during the formation of acid mine drainage, pyrite is destroyed
and goethite is formed so the first step is to put pyrite on the left side of the
equation and goethite and sulfate on the right.
FeS2 = FeOOH + SO42−
3. However, there must be two sulfates on the right-hand side to balance the two
sulfur atoms in the pyrite:
FeS2 = FeOOH + 2 SO42−
4. Next, there must be 0.5 H2O on the left-hand side to balance the H in the
goethite on the right-hand side.
FeS2 + 0.5 H2O = FeOOH + 2 SO42−
5. Because the sulfur in the pyrite is reduced, with a formal charge of −1, and the
sulfur in the sulfate is oxidized, with a formal charge of +6, we know that an
oxidant, in this case O2, must occur on the left-hand side of the reaction. There
are already 10 oxygen atoms on the right-hand side of the reaction and only
0.5 on the left so 9.5 oxygen atoms (4.75 O2) are needed on the left-hand side.
FeS2 + 4.75 O2 + 0.5 H2O = FeOOH + 2 SO42−
6. Finally, four hydrogen ions are needed on the right-hand side to charge bal-
ance the sulfate. These can come from the addition of two more water mol-
ecules to the left-hand side.
FeS2 + 4.75 O2 + 2.5 H2O = FeOOH + 4 H+ + 2 SO42−
However, this adds too much oxygen to the left-hand side, so one O2 must be
removed.
FeS2 + 3.75 O2 + 2.5 H2O = FeOOH + 4 H+ + 2 SO42−
7. This reaction is often written with decimal fractions as shown above but the
coefficients can be written as simple fractions.
FeS2 + 15/4 O2 + 5/2 H2O = FeOOH + 4 H+ + 2 SO42−
Then the relationships between the coefficients that are required by the conser-
vation of mass and charge are expressed by simple equations.
Fe: a = d
S: 2a = f
O: 2b + c = 2d +4f
H: 2c = d + e
chg: e = 2 f
12 Modeling tools
These equations are solved in a stepwise fashion to find the numerical value for
each coefficient.
Let a = 1.
From Fe: d = a = 1.
From S: f = 2a = 2.
From chg: e = 2f = 4.
From H: 2c = 1 + 4 so c = 5/2
From O: 2b = 2d +4f – c = 2 + 8 – 5/2 = 15/2 so b = 15/4
In all cases, the value for one of the coefficients must be arbitrarily chosen. The
other coefficients are then scaled to that value.
Notation
The notation used in this book follows geochemistry conventions as closely
as possible but has been modified where necessary to avoid redundancy or
confusion. Geochemical modeling requires an easy to recognize and intern-
ally consistent set of notation. Inconsistent notation complicates model
verification and leads to errors. The notation used in geochemistry is often
different from that used in engineering, chemistry, and biology, and some-
times models require the introduction of new notation. In all cases, notation
should be clearly defined in the model’s documentation and it is good prac-
tice to include a table of notation, with units, in every report.
Thermodynamic functions, especially free energy, enthalpy, and entropy,
are often used in kinetics models. Because thermodynamics is so widely used
in science, many different kinds of notation have developed. The convention
chosen for this book is illustrated using Gibbs free energy as an example.
For the free energy change related to a reaction, the delta in front of ∆Gr°
means that this variable reflects the change in free energy in going from the
reactants to the products. The subscript, r , designates this variable as the
free energy of reaction. The superscript, °, indicates that both the reactants
and the products are in the standard state, 298.15 K (25°C) and either 1 atm
or 100 kPa (1 bar) total pressure. For the free energy of formation of a sub-
stance from the elements, the subscript, f, associated with ∆Gf° means that
this is the free energy related to the formation from the elements (this is
essentially the same as ∆Gr° except that the free energy of formation of the
elements in their most stable state at 298.15 K (25°C) and either 1 atm or
1 bar total pressure is defined as zero. It is often useful to use an appended
Concentration 13
notation for ∆Gf° and ∆Gr° to indicate what reaction or substance the vari-
able refers to. For example, the free energy of reaction for a reaction listed
in equation (1) would be ∆Gr°(1) and the free energy of formation of quartz
would be ∆Gf°(qz).
Concentration and activity notation can be very confusing. Chemists have
long used square brackets, [i], to indicate the activity of species i and paren-
theses, (i), to indicate the concentration of species i, but geochemists have
occasionally reversed this convention. To make things worse, curly brackets
are used in some engineering disciplines to indicate activity. This confusion
is best avoided by abandoning the bracket-parenthesis notation in favor of
using ai to denote activity and mi to indicate molal concentration. Although
molar concentration units are best avoided in geochemistry, when used they
can be represented as Mi. Concentrations given other units such as mol/m3,
ppm, or mg/L units are best represented as ci.
pX notation indicates the negative logarithm (base 10) of quantity X.
Most people are familiar with pH, which is the negative logarithm of the
hydrogen ion activity, i.e. pH = –log aH . By analogy pK = –log K and
+
pCl = –log a Cl .
−
Units
The International System of Units was devised to facilitate communication
between scientific and engineering disciplines. Using SI units is highly rec-
ommended because it minimizes errors that arise from complicated units
conversions. The International System of Units is based on seven base units
(Table 2.1). Decimal multiples of these units and derived units are named
using the prefixes listed in Table 2.2.
Concentration
Concentration units listed in Table 2.3 are frequently used by geochemists.
Thermodynamic calculations are typically based on molal concentration
and concentration should be expressed in molal units whenever possible.
Some disciplines use molar or millimolar concentration units supposedly
because of the convenience of mixing of solutions using volumetric flasks.
Unfortunately molar concentrations change whenever the solution density
changes because of changing temperature, pressure, or composition. This
14 Modeling tools
SI base units
length meter m
mass gram g
time second s
electric current ampere A
temperature (thermodynamic) kelvin K
amount of substance mole mol
luminous intensity candela cd
Derived units
volume liter 1 L = 10–3 m3
force newton 1 N = 1 kg m–1
sec–2
energy joule 1J=1Nm
calorie1 1 cal = 4.184 J
pressure pascal 1 Pa = 1 N m–2
bar 1 bar = 105 Pa
atmosphere 1 atm = 101325 Pa
power watt W
electrical charge coulomb C
electrical potential volt V
1
Widely used but not recommended.
Table 2.3. The most commonly used concentration units for aqueous solutions
Unit Definition
This means that both CO2(aq) and H2CO3(aq) are present in the solution.
Because there is no simple way to determine the concentration of these
individual species, the concentration of unionized, dissolved CO2 is usu-
ally reported as the sum of the concentrations of CO2(aq) and H2CO3(aq).
This convention is quite useful for thermodynamic calculations and the
reported equilibrium constants for reactions involving H2CO3(aq) as well
as the ∆Gf°(H2CO3(aq)) are based on the convention of lumping these two
species together. However, this convention may cause a problem in kinetics
models because the rate of reaction of CO2(aq) is likely to be different from
the rate of reaction of H2CO3(aq).
16 Modeling tools
Dimensional analysis
A dimension is a quantity such as time, mass, temperature, or distance,
that can be expressed numerically in terms of one or more standard units.
Dimensions have homogeneous linear scales so that a 2-meter interval mea-
sured at the beginning of a 1-kilometer traverse is identical in length to a
2-meter interval measured near the end of the traverse. The magnitude of
a dimension is expressed in terms of a pure number and a unit of measure.
Units are standards to which dimensional quantities are compared to obtain
a numerical ratio, which is the pure number. For example, comparison of
18 Modeling tools
80 mg 1g 1 mol mol
80 ppm ≈ × × = 2 × 10 −3 = 2 × 10 −3 m (2.4)
1 kg 1000 mg 40 g kg
The universal gas constant can be converted from volume and pressure units
to energy units.
The surface area to volume ratio of a sphere is used in several models in this
book. Dimensional analysis can be used to find the relationship between the sur-
face area (A) with dimensions of [L2] and the internal volume (V) with dimen-
sions of [L3] for any convex regular solid. The surface area is related to the vol-
ume by a dimensionless constant, b.
[ L2 ] = b[ L3 ]2 3 (2.6)
A = bV 2 / 3 (2.7)
The first step toward finding the value of b, substitutes the definitions of A and
V for a sphere into Eq. (2.7).
2 /3
4
4π r 2 = b π r 3 (2.8)
3
4π 4π 1/ 3
b= 2 /3
= 2 /3
= 4.84 (2.9)
4 4
π 2 /3
3 3
Logarithms
Logarithmic transformations are quite common in geochemical models. The
variables used in geochemical calculations often range over many orders
of magnitude so that it is more convenient to work with the logarithms of
these numbers rather than the numbers themselves. In addition, many rela-
tionships in thermodynamics and kinetics are linearized using logarithmic
transformations.
A logarithm is the power (p) to which a base (a) must be raised to produce
a number (N).
ap = N (2.10)
Note that the largest errors associated with this method are
for the logarithms of 3 and 6, both of which are incorrect by
~0.02 log units.
Errors
There is some uncertainty in all data, and model building must take this
error into account. The first step in error management is error detection,
error reduction, and error quantification. There are three types of error:
systematic error, random error, and blunders. Improved experimental pro-
tocol can reduce all these, but designing progressively better experiments
eventually leads to diminishing returns so that at some point it is necessary
to use some kind of error analysis to manage the uncertainty in the variable
being quantified.
Systematic errors affect the accuracy but not the precision of the result.
They are usually errors in calibration or observation where the same incor-
rect protocol is applied to all measurements. They displace all measurements
from the true value by the same amount so they cannot be detected by a
statistical analysis of only one data set. However, systematic error can be
detected and reduced by comparing data sets from several different sources
using meta-analysis and systematic review (Rimstidt et al., 2012).
Random error displaces individual measurements from the true value, but
as the number of determinations increases their average becomes closer to
the true value. In a series of independent measurements that contain only
random errors, the most reliable estimate of the true value is the mean of all
the measurements. So that for a series of N measurements, x1, x2, x3, … xN,
the best estimate of the true value of x is x .
∑x i
x= i =1
(2.11)
N
22 Modeling tools
∑ (x − x)
2
i
σx = (2.12)
N −1
σx
σx = (2.13)
N
xs − x
ts = (2.14)
σs
3. Use a probability table to determine the probability (Ps) that a normally distrib-
uted measurement will differ from the mean by ts standard deviations (Ps = 1
− Ptable).
Significant digits 23
4. Determine the number of measurements in the data set that are expected to be as
bad as xs by multiplying the total number of data (N) by the probability (Ps).
N s = NPs (2.15)
5. If Ns < 0.5, then xs fails Chauvenet’s criterion and is rejected along with all data
that differ from the mean by more than xs.
6. After the data are rejected, calculate a new mean and uncertainty using the
remaining data.
7. Chauvenet’s criterion can be applied to a data set only one time.
Significant digits
The significant digits in a number give an indication of the uncertainty asso-
ciated with that value. This means that the number has been rounded to the
same order of magnitude as its uncertainty. The number of significant digits
is determined by the following rules:
1. The leftmost nonzero digit is the most significant one.
2. In numbers with no decimal point, the rightmost nonzero digit is the least
significant one.
3. In numbers with a decimal point, the rightmost digit is the least significant one,
even if it is a 0.
4. All digits between the least and most significant are counted as significant digits.
For example, 65000 has two significant digits; 0.1010 has four significant
digits; and 65000.1010 has nine significant digits. When numbers are trun-
cated to the proper number of significant digits, the least significant digit
must be rounded up or down according to the following rules:
1. If the digit following the one to be retained is 6, 7, 8, or 9, round the retained digit
up by one.
2. If the digit following the one to be retained is 0, 1, 2, 3, or 4, simply drop it with-
out changing the retained digit.
3. If the digit is 5 and the digit to be retained is even, do not change it. If it is odd,
round it up by one.
reported uncertainty, the least significant digit in the stated value should be
of the same order of magnitude (in the same decimal position) as the most
significant digit of the uncertainty. It is customary to perform the calcula-
tion using numbers with additional significant digits, to avoid introducing
rounding errors and then rounding the final answer to the correct number of
significant digits. For a number without a reported uncertainty, one should
assume that all reported digits are significant. Geochemical measurements
typically have more than 1% error so that, if no uncertainty is stated, three
significant digits are likely to be sufficient to express the reliable part of the
number.
Propagation of errors
Precision is lost as variables are manipulated through formulae. This loss of
precision is taken into account by propagating the uncertainty associated
with the model’s variables through the equations. There are several standard
references that provide detailed discussions of error propagation and man-
agement (Bevington, 1969; Deming, 1943; Mandel, 1964; Taylor, 1982).
For the propagation of uncertainties in equations with addition or sub-
traction, the uncertainty in the calculated value is the square root of the sum
of the squares of the absolute uncertainties.
q = x + y + za − b − c (2.16)
xyz
q= (2.18)
abc
2 2 2 2 2 2
δq δx δ y δz δa δb δc
= + + + + + (2.19)
q x y z a b c
q = Bx (2.20)
Propagation of errors 25
δq = B δx (2.21)
q = q (x) (2.22)
dq
δq = δx (2.23)
dx
log e
δq = δx (2.26)
x
0.434 δx
δ (log x ) = δ x = 0.434 (2.27)
x x
q = xn (2.28)
δq δx
= n (2.29)
q x
q = f ( x, y,... z ) (2.30)
26 Modeling tools
2 2 2
∂q ∂q ∂q
δq = δx + δ y + ... + δ z (2.31)
∂x ∂y ∂z
δq δx δ y ... δz
≤ + + + (2.32)
q x y z
δ q = 2 2 + 2 2 + 2 2 = 3.46 (2.33)
2 2
δq 1.22 × 10 −6 0.00131
= −4
+ = 0.0118 (2.34)
q 2 .31 × 10 0.124
The fractional error in the rate is 1.18% so the absolute error is (0.0118)(1.53 ×
10−5) = 1.81 × 10−7. This means that the rate is 1.53(0.0181) × 10−5 mol/sec.
Multiplication by a constant with no error: A common laboratory activity is to
prepare a solution by weighing a mass of solid and then dissolving it in water.
The concentration is usually reported in terms of the number of moles of solid
rather than the mass. If the mass of NaCl (WM = 58.44) is 40.00(0.01) g, the
number of moles is 40.00/58.44 = (40.00)(0.0171) = 0.6845 moles. We can assume
that there is no random error associated with the molecular weight of NaCl, so
the error for this value is found using Eq. (2.21).
This conversion shows that the amount of NaCl added to the solution is
0.6845(0.00017) mol.
Regression models 27
Regression models
Fitting data to equations is an important modeling activity that is supported
by an abundance of computer software. The mathematical basis for this
software is described in many statistics books and will not be elaborated
here. Methodologies that are especially suited for kinetics data are reviewed
in Chapter 7 of van Boekel (2009). This section explains some strategies that
can guide fitting kinetic data to equations.
One reason to fit data to a function is to summarize the data and report it
in a way that allows easy interpolation or other mathematical manipulation.
If this is the main objective, the data can be fit to an arbitrary function so
long as that function makes a smooth and parsimonious summary of the
data and is compatible with desired mathematical manipulations. Typically
a low-order polynomial is suitable for this purpose. It is important to avoid
over-fitting the data. The goal is to filter out random error without losing
the information contained in the data. Increasing the number of regression
variables tends to fit the data better as reflected by smaller residuals and a
higher correlation coefficient. This is because the additional variables make
the fitting function more flexible, with more minima and maxima, so it can
pass closer to the data points. However, this additional flexibility can cause
the function to flex wildly and to predict values that are beyond the range
of the data. This unreasonable result can be avoided by beginning the fitting
process using the simplest feasible equation and then adding variables one
at a time until the residuals of the fit are about the same size as the errors in
the data. Sometimes data can be converted to a linear or other simple form
using a logarithmic or power transformation. The fitted equation should not
be extrapolated beyond the domain of the original data.
28 Modeling tools
Table 2.4 lists Fe3+ concentration versus contact time with pyrite from a plug
flow reactor experiment (PFR10) (Rimstidt and Newcomb, 1993). The rate of
ferric iron consumption is the time derivative of these data. Finding this deriva-
tive for t = 0 is especially useful because the initial conditions of the experiment
are accurately known so their effect on the rate can be clearly established. One
strategy for finding the initial rate is to fit the first few concentration versus time
0 1.80 × 10−3
4.7 1.16 × 10−3
6.3 1.12 × 10−3
9.1 1.04 × 10−3
20.5 8.24 × 10−4
0.002
0.0015
• •
mFe3+, molal
0.001 •
0.0005
0
0 5 10 15 20 25
t, sec
Figure 2.1. Fit of the PFR10 data to a second-order polynomial (solid line, R 2 =
0.96) and a fourth-order polynomial (dashed line, R 2 = 1.00).
Regression models 29
dmFe3+
= b + 2 ct + ... (2.39)
dt
When t = 0, the rate equals b. Figure 2.1 shows fits of a second and a fourth-
order polynomial to the data.
The second-order fit predicts an initial rate of −1.23 × 10−4 m/sec and the fourth
fit predicts a much faster rate of −4.85 × 10−4 m/sec. Both functions fit the data
well and the fourth-order polynomial fits it so well that the function passes
through all the points so there are no residuals. However, the rate based on the
fourth-order polynomial appears to be too fast because that function does not
smooth out the random errors in the data. Both functions incorrectly predict
that the concentration of Fe3+ increases after ~18 seconds of reaction because
they are simply approximations for a short interval of the true concentration
versus time function.
A second reason to fit data to a function is to test whether the data are con-
sistent with a model or to use a theoretical function to extrapolate experimen-
tal results to conditions otherwise not attainable. The goodness of fit of the
data to a theoretical function can be gauged by the coefficient of determin-
ation (R2), which is the correlation coefficient squared. R2 is often interpreted
as the fraction of the variability of the response variable that is explained
by its functional relationship to the independent variable. For example, if
R2 = 0.8 then 80% of the variation in the response variable is explained by
the model and 20% of the variation is the result of factors, including ran-
dom error, that are not part of the model. Some geochemical processes occur
under conditions or over time intervals that are difficult to simulate by exper-
iments, making it necessary to use a theoretical model to predict their behav-
ior. This approach uses data from experiments performed under easily attain-
able conditions to quantify parameters that are used to calibrate a theoretical
model. These parameters are then used to predict the value of the variable at
30 Modeling tools
0.8
0.6
T, °C k, sec–1
–6
–8
•
••
–10 •
–12
log k
–14
–18
100°C
–20
0.0015 0.002 0.0025 0.003 0.0035
1/T, K–1
Figure 2.3. Arrhenius plot used to extrapolate rate constants for the kaolinite to
illite transformation reaction from high temperature measurements to find the
rate constant at 100°C.
k = Ae − Ea RT
(2.42)
− Ea 1
log k = log A + (2.43)
R T
32 Modeling tools
Fitting the data in Table 2.5 to this equation gives the line shown in Figure 2.3.
This fit explains 97% of the variation in log k (R2 = 0.97).
−8039(992 )
log k = + 5.45(1.8) (2.44)
T
So log k = −16.1(2.7). Even though the data fit the equation quite well, the
uncertainty in log k is nearly 3 log units. Most of this error is due to the uncer-
tainty associated with the slope of the line. This example stands as a warning to
model builders who extrapolate data without considering error propagation.
Numerical differentiation
Kinetic processes cause a change in a quantity over time. If this quantity is
measured at evenly spaced time intervals, it is relatively easy to find the rate
by computing a numerical derivative (Pollard, 1977). The simplest method
of numerical differentiation is based on the observation that the chord of a
graph has a slope that closely approximates the tangent at an intermediate
point.
This leads to a simple scheme for calculating the numerical derivative
from a data set (Pollard, 1977). The derivative at the first point in the data
set ( f−′1) gives more weight to the first two points.
1
f −′1 = ( −3 f−1 + 4 f0 − f1 ) (2.46)
2h
The derivative at each intermediate point ( f0′ ) is the slope of the chord con-
necting the surrounding points (Figure 2.4).
1
f0′ = ( − f−1 + f1 ) (2.47)
2h
The derivative at the last point in the data set ( f+′1) gives more weight to the
last two points.
1
f+′1 = ( f−1 − 4 f0 + 3 f1 ) (2.48)
2h
Numerical differentiation 33
f–1 f0 f+1
h
Figure 2.4. Schematic
f(x) f 0 = 1 ( f 1+ f1)
2h illustration of how
the slope of a chord
approximates the
tangent of a graph at an
x intermediate point.
Note that the equations for the endpoints of the data set are less dependable
than the equation for the intermediate points and they tend to give unrea-
sonable results if the slope is steep. Because numerical differentiation takes
the difference between values of the dependent variable it magnifies the scat-
ter in those data. Some of this scatter can be eliminated by various smooth-
ing schemes, but they tend to bias the data and are best avoided.
Ferric iron reacts with pyrite to produce ferrous iron and sulfuric acid.
14 Fe3+ + FeS2(py) + 4 H2O = 15 Fe2+ + 2 SO42− + 8 H+ (2.49)
1
f −′1 = ( −3 (9.79 × 10 ) + 4 (9.61 × 10 ) − 9.45 × 10 ) = −6.33 × 10
−4 −4 −4 −8
( 2 )(300 )
(2.50)
The rate for the next point (t = 300) is found using Eq. (2.47).
1
f0′ = ( −9.79 × 10 + 9.45 × 10 ) = −5.67 × 10
−4 −4 −8
(2.51)
(2 )(300 )
The subsequent points are treated the same way. The rate for the last point (t =
3000) is found using Eq. (2.48).
1
f+′1 = (9.04 × 10 − 4 (8.98 × 10 ) + 3 (8.91 × 10 )) = −2.50 × 10
−4 −4 −4 −8
( 2 )(300 )
(2.52)
34 Modeling tools
9.5× 10–4
concentration, molal
9.0× 10–4
8.5× 10–4
0 1000 2000 3000
time, sec
(b) 0.0
–2.0× 10–8
rate, molal/sec
–4.0× 10–8
–6.0× 10–8
Figure 2.5. (a) Concentration of Fe3+ versus time from experiment BR5 (Rimstidt
and Newcomb, 1993). (b) Numerical derivatives of the concentration of Fe3+
versus time data. Early in the experiment the rates are fast but they slow to a
nearly constant rate after ~300 sec.
Numerical differentiation 35
dn
r= = k ∏ aini (3.1)
dt
no = ∑ ni (3.2)
36
Rate equations 37
dm r 1
R= = = k ∏ aini (3.3)
dt M M
If the reaction takes place on an interface that separates two phases and the
rate equation expresses how fast a component is transferred to or from that
interface, this quantity is actually a flux (J, mol/m2sec).
r 1
J= = k ∏ aini (3.4)
A A
From Eq. (3.3) we see that r = RM and from Eq. (3.4) we see that r = JA. If
we combine these expressions we find how the rate of change of concentra-
tion in a solution is related to the dissolution flux.
dm A A
R= = J= k ∏ aini (3.5)
dt M M
Equations (3.3) and (3.5) contain both activity and concentration variables.
In order to integrate the equation and to make the equations compatible
with transition-state theory (Chapter 5), the concentration in the derivative
term of these equations can be converted to activity.
γ m
a= (3.6)
γ ° m°
38 Rate equations
da 1 γ γ k
= r = ∏ aini = k ′∏ aini (3.7)
dt m°M γ ° M
da A γ A
= J = γ k ∏ aini = k ′∏ aini (3.8)
dt m°M γ ° M
More often, the activity terms on the right-hand side are converted to con-
centration terms by rearranging Eq. (3.6). This allows the equation to be
integrated to give the concentration as a function of time, but can result
in complicated and often meaningless units for the apparent rate constant.
For example, if one or more of the reaction orders is a fraction, as often
is the case for mineral dissolution rates, the apparent rate constant will
have units of molal raised to a fractional power, which has no physical
meaning.
dm k k ∏ γ ini
R= = ∏ γ ini mini = ∏ mi = k ′∏ mi
ni ni
(3.9)
dt M M
dm A A
R= == k ∏ γ ini mini = k ∏ γ ini ∏ mini = k ′∏ mini (3.10)
dt M M
( ni )t = ( ni )t = 0 + νiξ (3.11)
The reaction rate is then expressed as the time derivative of the extent of
reaction.
d ξ 1 dni
= (3.12)
dt ν i dt
Unopposed reactions (n ≠ 1) 39
The oxidation of pyrite by ferric iron produces ferrous sulfate and sulfuric acid.
14 Fe3+ + FeS2(py) + 8 H2O = 15 Fe2+ + 2 SO42− + 16 H+ (3.13)
The reaction rate can be written in terms of any of the species in Eq. (3.13)
where the negative sign indicates that the species is consumed.
A R
+
→ products (3.15)
40 Rate equations
In the simplest case, only one aqueous species affects the reaction rate and
the rate of consumption of species A is proportional to the rate constant (k)
multiplied by its concentration (m, mol/kg) raised to an arbitrary power (n).
Values of n and k are determined empirically.
dm
= − km n (3.16)
dt
m t
dm
∫m n
= − k ∫ dt (3.17)
mo 0
(− n +1)
m
= 1−
( − n + 1) kt (3.19)
m mo(− n +1)
o
m
p= (3.20)
mo
1
( − n + 1) kt (− n +1)
p = 1 − (3.21)
m(− n +1) o
( − n + 1)kt
( 1− n +1)
α = 1 − 1 − (3.22)
m( − n +1) o
When air-saturated surface waters infiltrate pyritic mine wastes, the dissolved
oxygen (DO) is consumed by reaction with pyrite.
FeS2 + 7/2 O2 + H2O = Fe2+ + 2 H+ + 2 SO42− (3.23)
Equation (3.21) can be used to find how long it will take for all the DO to be
consumed if 1 kg of air-saturated solution (mDO = 2.7 × 10−4 = 10−3.57 mol/kg)
contacts 1 m2 of pyrite surface.
Unopposed reactions (n ≠ 1) 41
0.8
0.6
p
0.4
0.2
0
0 10 20 30 40 50 60 70
time, hr
Williamson and Rimstidt (1994) report the following rate equation for reaction
(3.23).
dξ 2
= − J py = − J DO (3.25)
dt 7
Equations (3.5), (3.24), and (3.25) can be combined to produce a rate equa-
tion for DO consumption for a system with 1 kg of water and 1 m2 of pyrite
surface area.
t=
mo(1− n ) (10 −3.57 )0.5
dm
= − km (3.28)
dt
m t
dm
∫ m
= −k ∫ t (3.29)
m0 0
m
ln = − kt (3.30)
mo
It is often useful to express this result in terms of the fraction reacted (α) or
the fraction remaining (p).
m
α = 1− p = = e − kt (3.31)
mo
Characteristic time
Reaction rates are often compared using a characteristic time. For example,
the characteristic time might be defined as the time needed for the destruc-
tion of 50% of the reactants (α = p = 0.5). This time is called the reaction’s
half-life, t½. When t = t½, α = p = 0.5. The half-life for an unopposed reac-
tion where n ≠ 1 is found by setting m = 0.5 in Eq. (3.18).
(− n +1)
mo(− n +1) − (0.5 mo )
t1 ⁄ 2 = (3.32)
( − n + 1) k
Characteristic time 43
ln 0.5 0.693
t1 ⁄ 2 = = (3.33)
−k k
After 5 tc, p = 99.3 and for all practical purposes the reaction can be said to
be complete.
The concentration of dissolved oxygen in equilibrium with air is 2.7 × 10−4 molal.
This can be combined with the expression for the rate constant above to find an
overall rate constant, k′.
44 Rate equations
200
180
160
140
120
t , hr
100
80
60
40
20 H2S HS–
0
4 5 6 7 8 9 10
pH
Combining Eqs (3.35) and (3.37) gives an equation that expresses the rate as a
function of hydrogen sulfide concentration.
m
ln = ln p = − k ′t (3.39)
mo
p = e − (k ′t) (3.40)
Substituting α = 0.5 into this equation gives a relationship for the half-life of
H2S (t½, hr) as a function of hydrogen ion concentration.
Opposed reactions: principle of detailed balance 45
A←
→B
R+
(3.43)
R−
For the very simple reaction described by Eq. (3.43), the rate of the for-
ward reaction is the product of the forward rate constant (k+) and the
activity of A.
dmB
R− = = k + aA (3.44)
dt
The rate of the backward reaction is the product of the backward rate con-
stant (k−) and the activity of B.
dmA
R+ = = k − aB (3.45)
dt
The overall reaction rate is the difference between the forward and back-
ward rates.
R = k+ aA − k− aB (3.46)
1.2
0.8
R/R+
0.6
0.4
Figure 3.3. Reaction rate
normalized to the far-from-
equilibrium rate versus 0.2
chemical potential driving
the reaction. The overall
rate approaches zero as 0
the driving force for the –10 –8 –6 –4 –2 0
reaction approaches zero. ∆µr /RT
k+ aB , eq
= =K (3.47)
k− aA, eq
This relationship appears to be valid even when the reaction is not at equi-
librium so long as the reaction mechanism (i.e. transition state) remains the
same. This means that the rate equation can also be written to make the rate
proportional to the activity ratio, Q = aA/aB.
ka Q
R = − k+ aA 1 − − B = − k+ aA 1 − (3.48)
k+ aA K
Q
∆µ r = − RT ln (3.49)
K
− ∆µ r
Q
= e RT (3.50)
K
Combining Eqs (3.48) and (3.50) gives an equation that relates the reac-
tion rate to the chemical potential driving the reaction. This relationship
can also be expressed in terms of affinity where A = −∆µr. Affinity is a key
concept in non-equilibrium thermodynamics, which is discussed briefly in
Chapter 10.
First-order forward and first-order backward 47
− ∆µ r
R = k+ aA 1 − e RT (3.51)
When the reaction is far from equilibrium (∆µr < ~−3 RT) the reverse reac-
tion rate is effectively zero, so the far-from-equilibrium reaction rate (R+)
≈ k+aA. This means that Eq. (3.51) can be written to show how the overall
rate declines as the reaction approaches equilibrium (Figure 3.3).
− ∆µ r
R R
= = 1 − e RT (3.52)
R+ k+ aA
This relationship is valid only if the reaction mechanism remains the same
over the entire range of conditions. If the reaction mechanism changes, the
relationship between R and µ becomes complicated.
RA = k+ mA
A← →B
(3.53)
RB = k− mB
The net rate of appearance of B is the difference between the forward and
reverse rates.
dmB
= RB − RA = k+ mA − k− mB (3.54)
dt
Aris (1989) presents three different schemes to integrate this differential rate
equation for an initial condition where the concentration of B is zero, i.e.
mBo = 0 when t = 0. The following derivation is a more general derivation
where mBo ≥ 0. The trick to this derivation is to identify relationships, which
allows Eq. (3.54) to be recast into a form that is easily integrated. The most
useful of these relationships are the conservation of mass and the principle
of detailed balance. The conservation of mass requires that the concentra-
tion of A and B must sum to a constant value so that the sum of concentra-
tions of A and B at any time must equal the sum of their concentrations at
equilibrium (mAe + mBe).
The principle of detailed balance requires that the quotient of the equilib-
rium concentrations be equal to the equilibrium constant. This applies for
dilute solutions where the activity coefficients are very near one.
k+ mBe
K= = (3.56)
k− mAe
The conservation of mass equation (3.55) can be used to write the mA term
in Eq. (3.54) in terms of mB.
dmB
= k+ (( mAe + mBe ) − mB ) − k− mB (3.57)
dt
dmB
= k+ ( mAe + mBe ) − ( k+ + k− ) mB (3.58)
dt
mB t
dmB
∫ k+ ( mAe + mBe ) − ( k+ + k− ) mB ∫0
= dt (3.59)
mB 0
k ( m + mBe ) − ( k+ + k− ) mB
ln + Ae = − ( k+ + k − ) t (3.60)
k+ ( mAe + mBe ) − ( k+ + k− ) mBo
+ ( K Be )
k mBe + m − k + k m
( + −) B
= − ( k+ + k − ) t
( )
ln (3.61)
k+ mBe + mBe − ( k+ + k− ) mBo
K
( k + k− ) mBe − ( k+ + k− ) mB mBe − mB
ln + = ln = − ( k+ + k − ) t (3.62)
( k + + k − ) mBe − ( k + + k − ) mBo mBe − mBo
mB
mBe − mB 1 − mBe
ln = ln = ( k+ + k− ) t (3.63)
mBe − mBo 1 − mBo
mBe
First-order forward and first-order backward 49
This equation can be transformed and rearranged to give the ratio of the
concentration of B to its equilibrium concentration.
mB mBo − (k+ + k− )t
m = 1 − 1 − m e (3.64)
Be Be
mB − (k+ + k− )t
m = 1 − e (3.65)
Be
0.5 =
mB
mBe
( − ( k + k )t
= 1 − e + − 12 ) (3.67)
CO2 dissolves in water as CO2(aq), which then reacts with water to form car-
bonic acid, H2CO3.
CO2 (aq ) + H2 O ←
R+
→ H2CO3 (3.69)
R−
dmH2 CO3
R+ = = k+ aCO2 aH2 O = k+ aCO2 (3.70)
dt
The rate of the reverse reaction is the reverse rate constant times the concentra-
tion of H2CO3.
50 Rate equations
0.0275 31
0.0257 25.9
0.026 20.6
0.043 25.5
0.0358 27
0.0375 15.1
0.044 17.5
0.036 13.7
0.037 18
16.7
Average 0.0347(0.002) 21.1(1.9)
dmCO2
R− = = k− aH2 CO3 (3.71)
dt
The reacting species have zero charge, so we can set their concentrations equal to
their activity. The overall rate of appearance of H2CO3 is the difference between its
rate of formation and the rate of destruction as defined by Eqs (3.70) and (3.71).
dmH2 CO3
R+ − R− = = k+ mCO2 − k− mH2 CO3 (3.72)
dt
This relatively small value means that at equilibrium only 0.16% of the dissolved
CO2 is in the form of H2CO3. This is due to the small size of the forward rate
First-order forward and first-order backward 51
1.0
0.8
0.6
mB/mBe
0.4
0.2 50 % 99 %
0
0 0.1 0.2 0.3
time, sec
Figure 3.4. Graph of mB/mBe versus time showing the time needed for the H2CO3
concentration reaction to reach 50% and 99% of the equilibrium value.
constant relative to the backward rate constant, i.e. H2CO3 forms slowly and
dissociates quickly.
Equation (3.66) can be rearranged to find the time needed for H2CO3 to reach
a fraction of the equilibrium concentration.
1 m
t=− ln 1 − B (3.74)
( k+ + k − ) mBe
1
t=− ln (1 − 0.5) = − (0.0473) ( −0.301) = 0.0327 sec (3.75)
(21.13)
If mB/mBe = 0.99,
1
t=− ln (1 − 0.99) = − (0.0473) ( −4.61) = 0.218 sec (3.76)
( .13)
21
Although these times may seem short compared to the characteristic times for
many geochemical processes, they are actually very slow compared to the rates
of interaction of most aqueous ions with water molecules.
52 Rate equations
daB
RB − RA = = k+′ − k−′ aB (3.79)
dt
k ′ − k−′ aB
ln + = − k−′ t (3.81)
k+′ − k−′ aBo
k−′ aB aB
1 − k+′ 1− K
ln = ln = − k−′ t (3.82)
1 − k−′ aBo a
1 − Bo
k+′ K
The equilibrium constant for Eq. (3.77) is equal to the equilibrium activity
of species B (K = aBe) so K in Eq. (3.82) can be replaced by aBe.
aB aBo − k ′ t
a = 1 − 1 − a e (3.83)
−
Be Be
aB
a =1− e
− k−′ t
(3.84)
Be
This equation can be rearranged to find a characteristic time (t½, sec), which
is defined as the time needed for aB /aBe to change from 0 to 0.5.
aB
0.5 = = (1 − e − k−′ t1 ⁄ 2 ) (3.85)
aBe
Quartz makes up about 15% of the Earth’s crust, so it is not surprising that it
often controls the dissolved silica concentration in adjacent pore waters, espe-
cially at high temperatures. However, at low temperatures the ability of quartz
to buffer the dissolved silica concentration is limited by its slow dissolution and
precipitation rates. We can visualize the time needed for a solution to equilibrate
with quartz as a function of temperature using the rate constants from Rimstidt
and Barnes (1980) to construct a graph showing the concentration (activity) of
dissolved silica versus time for solutions that initially contain no dissolved silica
at temperatures of 25°, 50°, and 75°C. This model considers a system consisting
of 1 m2 of quartz in contact with 1 kg of water.
Quartz dissolves to release silicic acid into solution.
SiO2(qz) + 2 H2O(l) = H4SiO4(aq) (3.87)
K = aH4 SiO4 e (3.88)
daH4 SiO4
= k+′ − k−′ aH4 SiO4 (3.89)
dt
T, °C K k−′
25 1.11 × 10−4 3.80 × 10−10
50 2.50 × 10−4 1.79 × 10−9
75 4.94 × 10−4 6.77 × 10−9
0.0005
0.0004 75°C
silica concentration, molal
0.0003
50°C
0.0002
0.0001
25°C
0
0 20 40 60 80 100
time, yr
Because H4SiO4 is an uncharged species, we can assume that its activity coef-
ficient is near 1, and because the activity of quartz and water are 1 we can set
aH4 SiO4 e = K .
This allows Eq. (3.91) to be rewritten in terms of the equilibrium constant and
the precipitation rate constant.
The values of the precipitation rate constant and the equilibrium constant at
25°, 50°, and 75°C (Rimstidt and Barnes, 1980) given in Table 3.2 were used to
calculate the concentration versus time graphs shown in Figure 3.5.
Zero-order forward and first-order backward 55
100
80
60
t1⁄2, yr
40
20
0
0 20 40 60 80 100
temperature, °C
Figure 3.5 shows that even though the equilibrium concentration of dis-
solved silica increases considerably with increasing temperature, the time needed
to reach 50% of that equilibrium concentration decreases significantly from
58 years at 25°C to 3.2 years at 75°C.
Another way to visualize this is to create a graph that shows t½ as a function
of temperature using the relationship t1 ⁄ 2 = 0.693 k−′ where log k− = −0.707–
2589/T (K) (Rimstidt and Barnes, 1980). This graph is shown in Figure 3.6. At
temperatures above 100°C the equilibration rates become relatively fast.
Chapter 4
Chemical reactors
Ideal chemical reactors can be linked in many different ways so that the
output of one reactor becomes the input of the next so that the reactions
occur in stages. Staged reactor models can be adapted to simulate most, if
not all, real-world scenarios. This means that the behavior of complex inter-
acting systems can be simplified to combinations of simple reactors, each
of which is relatively easy to model. Ideal chemical reactor models can be
applied as a first approximation to most natural situations and the concept
easily leads to mathematical descriptions of those situations. For example,
56
Chemical reactors 57
a shallow pond or lake that is well stirred by wind action can be modeled
as a batch reactor. A pool in a low gradient stream, which is stirred by the
stream flow, can be modeled as a mixed flow reactor; and a parcel of fluid
flowing in a rock fracture can be modeled as a plug flow reactor. In addition
to using chemical reactors models to understand natural situations, they are
the basis for designing experiments to measure rates.
Chemical reactors can be embellished to provide effective mixing and
to establish well-defined physical relationships among reacting phases.
Figure 4.2 illustrates some examples of methods devised to bring solids into
contact with solutions under controlled hydrodynamic conditions. Mixing by
an overhead stirrer can suspend small particles to produce a homogeneous
slurry. Magnetic stirrers can also be used for mixing, but the motion of the
magnet on the reactor bottom can grind the particles, which changes the
surface area of the solid to the mass of the solution in an unpredictable way.
Shaking, rocking, or rotating a reactor that has a free space (gas phase) can
also produce effective mixing. This scheme is especially useful for batch reac-
tors that must be sealed to avoid the gain or loss of a vapor phase. Reactors
containing a slurry, packed bed, or monolith can be mixed using an external
recycle loop containing a pump that provides a mixing flow. If the loop is
closed, the system is a batch reactor. If solution is added and removed from
the loop, the system is a mixed flow reactor. Larger grains can be held in a
packed bed through which the solution is forced by pumping. If the bed is
long compared to its lateral dimensions, the reactor will perform as a plug
58 Chemical reactors
overhead externally
stirrer recycled
packed bed
shaken
Figure 4.2. Examples of
schemes bring solids
into contact with fluids in magnetic rotating disk
chemical reactors. stirrer
flow reactor. If the bed is thin relative to its lateral dimensions, it will approxi-
mate the behavior of a mixed flow reactor. Monolith samples, which can be
solid blocks of material or particles imbedded in an inert solid matrix, are
desirable because they can be cleaved or ground and polished to present well-
defined surfaces to the solution and these surfaces can be observed during and
after their reaction with the solution. Monoliths can be presented to the solu-
tion in several ways ranging from simply suspending them in a stirred solu-
tion to flowing the solution past them using various schemes. A widely used
scheme is to embed the monolith into the face of a disk, which is attached to
a rotating shaft. This rotating disk exposes the solid to a well-defined hydro-
dynamic situation, which makes this apparatus useful for studying the rate
of transport-limited dissolution (Alkattan et al., 1998; Levich, 1962; Sjöberg
and Rickard, 1983). Many other schemes are possible, but regardless of the
scheme, the reactor should be designed to operate as closely as possible to an
ideal chemical reactor to simplify modeling its behavior.
Tracer kinetics 59
Tracer kinetics
Tracer kinetics interpreted using the continuity equation is an effective
way to characterize the mixing and flow behavior of chemical reactors.
Tracers are inert substances that are added to reactors, and the tracer con-
centration in the reactor or in the effluent stream is used to infer some-
thing about the reactor performance. Ideal tracers are inert so they do
not react with anything in the reactor. This means that the rgen term in the
continuity equation (4.1) equals zero. Tracers are frequently used to deter-
mine reactor dimensions or flow rates in field settings where direct phys-
ical measurements are impractical or impossible. In addition to displaying
conservative behavior, field tracers must be environmentally benign, inex-
pensive, and easy to quantify at low concentrations. Bromide and lithium
ions meet these criteria and are frequently used. The environmental and
cost requirements for tracers used in laboratory, pilot scale, and indus-
trial settings are less stringent because smaller amounts are used and the
resulting solutions can be captured and treated before disposal. Dyes or
even radioactive tracers might be used under these more controlled con-
ditions. Natural situations are often complicated so the models in this
section should be viewed as approximations of real systems. The models
described by Rubinow (1975) illustrate how tracer models are applied to
more challenging cases.
Tracer dosing follows one of two end member scenarios; the tracer is either
added as an instantaneous spike or it is added continuously beginning at a
known time. Adding spikes is the simplest method but if the reactor is large
a single spike may not produce a detectable signal. On the other hand, con-
tinuous addition requires a metering pump and a larger amount of tracer.
Serendipitous natural tracers can be used, but they are often introduced in
unpredictable and irregular amounts. This makes them useful for tracing
flow paths even if they cannot be used to develop quantitative models of
reactor volumes or flow rates.
Tracer behavior is different in each of the ideal chemical reactors, so the
first step in developing a tracer kinetics model is to decide which reactor
type best simulates the real situation and which kind of tracer dosing has
occurred. If the amount and timing of tracer introduction is known, a for-
ward model can be developed. If the amount and timing of tracer detected
is known, an inverse model can be produced. The equations derived in the
following section use concentration units of mg/L because these units are
typical of field studies. For laboratory experiments the models would use
molal concentration units. The mass (M) variables would be replaced by
mole quantities (n). The volume variables (V) would be replaced by mass of
water (M) and the flow rate (Q) would have units of kg/sec.
60 Chemical reactors
The concentration in the reactor after the spike is added (C, mg/L) is the
sum of the initial mass of tracer (Mo, mg) and the mass of tracer added in
the spike (Ms, mg) divided by the total volume of solution (Vs, L).
M o + M s CoVo + C sVs
C = = (4.3)
Vo + Vs Vo + Vs
The initial amount of the tracer in the reactor is the original concentration
of tracer (Co, mg/L) times the original volume of solution (Vo, L).
M o = CoVo (4.4)
The mass of the tracer in the spike is its concentration (Cs, mg/L) times the
volume of the spike solution.
M s = C sVs (4.5)
Hot spring pools are usually well stirred by thermal convection and rising gas
bubbles, so they can be modeled as ideal batch reactors and the volume of water
in the pool can be determined using a tracer spike.
A typical spike might consist of 1 L of a 4.6 × 105 mg/L bromide solution. Once
the spike is added to the hot spring, it takes a short time for it to mix throughout
the spring pool. If the completely mixed pool water contains 4.6 mg/L, Eq. (4.3)
can be rearranged to calculate the volume of the pool (Vo,, L). This relationship
can be further simplified by realizing that Cs >> C so C can be neglected in the
numerator and that C >> Co so that Co can be neglected in the denominator.
C −C C
V = Vs s ≈ Vs s (4.6)
C − Co C
mg
4.6 × 105
Cs L = 1 × 105 L
V ≈ Vs = 1 L (4.7)
C mg
4.6
L
The rate of removal of the tracer by the effluent solution is directly pro-
portional to its flow rate (Q, L/sec) divided by the reactor volume (V, L).
This quotient (kQ, sec−1) behaves like the rate constant in a first-order rate
equation.
dC Q
= − C = − kQC (4.9)
dt V
Equation (4.9) can be rearranged and integrated from the initial tracer
concentration (Co, mg/L) to the tracer concentration at some later time
(C, mg/L).
C t
dC
∫ C = − kQ ∫ dt (4.10)
Co 0
C
ln = − kQt (4.11)
Co
Once the tracer concentration has homogenized in the hot spring, as described
in Example 4.1, discharge from the spring carries away some of the tracer so the
concentration in the spring pool declines as shown in Figure 4.3.
Equation (4.11) can be used along with the data from this graph to find the
discharge of the hot spring. After 5 hours (18,000 sec) have passed, the concen-
tration is 0.76 mg/L so C/Co = 0.165 and ln(C/Co) = −1.80. This means that kQ =
−1.8/−1.8 × 104 = 1.0 × 10−4 sec−1. The volume of water in the spring pool is 1.0
× 105 L so the discharge Q = kQV = (1.0 × 10−4 sec−1)(1.0 × 105 L) = 10 L/sec.
62 Chemical reactors
3
C, mg/L
0
0 2 4 6 8 10
time, hr
Qt
V= (4.13)
φ
The mass of tracer in the reactor (M, mg) is the initial mass (Mo, mg) plus the
mass added over the interval of tracer injection (ΔM, mg) and that amount
equals the product of concentration in the tracer feed (Cin, mg/L) and the
feed flow rate (Qin, L/sec) times the time interval (t, sec).
M = M o + ΔM = M o + CinQint (4.15)
The volume of solution in the reactor at the end of the tracer injection (V,
L) is the initial volume (Vo, L) plus the volume added by the injection flow.
V = Vo + Qint (4.16)
The concentration of tracer (C, mg/L) in the reactor at time = t (sec) is the
mass of tracer from Eq. (4.15) divided by the volume of solution from Eq.
(4.16).
M M o + CinQt
C= = (4.17)
V Vo + Qt
dC Q Q
= Cin − C = kQ (Cin − C ) (4.19)
dt V V
C t
dC
∫ Cin − C = kQ ∫ dt (4.20)
Co 0
64 Chemical reactors
C − Co
ln in = kqt (4.21)
Cin − C
Vφ
t= (4.23)
Q
All subsequent slugs contain the same tracer concentration as the first.
data to integrated rate equations, either using linearized forms or using non-
linear regression, will produce dubious results unless the data span a very
large extent of reaction (Rimstidt and Newcomb, 1993).
A more common approach to analyzing batch reactor data is to deter-
mine rates by taking the derivative of the concentration versus time data.
If the concentration is determined at equally spaced times, the numerical
derivative of the concentration versus time data can be used to find the rate
of generation at each sampling time. A rate equation can then be developed
by correlating the rate at each sample time with the conditions that existed
at that time. Single experiments seldom span a large enough extent of reac-
tion to produce a reliable correlation, so rates from multiple experiments
that span a large range of concentrations must be pooled.
An alternative strategy is to perform numerous short-term experiments
that span a wide range of initial conditions and then use the derivative of
a polynomial fit of the concentration versus time data for each experiment
to find the rate at the beginning of each experiment. An advantage of find-
ing this initial rate is that all the rate-controlling parameters are fully con-
strained by the experimental design at the beginning of the experiment. This
method finds the derivative, at t = 0, of an arbitrary function that fits the
concentration versus time data. The simplest fitting function is a straight line
and the slope of the line equals the rate. This works for cases where the rate
is nearly constant over the duration of the experiment. If the rate changes
significantly during the experiment, the concentration versus time graph will
be curved and a straight line underestimates the rate. Fitting the data to a
second (or higher) order polynomial will accommodate some curvature.
m = a + bt + ct 2 (4.24)
dm
= b + 2ct (4.25)
dt
the initial rate from the slopes, (mt − mo)/∆t, of chords drawn between each
concentration (mt) point and the starting concentration (mo). The slopes of
these chords are graphed versus time and this function is extrapolated to t =
0. This gives the slope of the tangent to the concentration versus time graph
at t = 0, which is the initial rate. Graphs of the chords’ slopes versus time
are nearly linear for small extents of reaction but do show some curvature,
which can be fit by a second-order polynomial. The initial rate is found by
solving this polynomial for the slope at t = 0.
t m ∆m/∆t
sec molal molal/sec
0 1.89 × 10 −6
The heavy line in Figure 4.4 is a second-order polynomial fit of all the concen-
tration versus time data in Table 4.1.
The “b” term for this fit is 6.89 × 10−10. A comparison of this fit to the first few
data points shows that the slope of the line at t = 0 is too small. This problem
can be remedied by fitting only the first three data points to produce the light
line shown in Figure 4.4.
Batch reactor experiments 67
8×10–6
•
•
6×10–6 •
•
concentration, molal
4×10–6 •
2×10–6 •
0
0 4000 8000 12000 16000
time, sec
The “b” term for this fit is 1.31 × 10−9, which is nearly twice the value obtained
by fitting all the data points. This example illustrates that a relatively large
amount of error can result from fitting data for too large an extent of
reaction.
The chord method is an improvement over the polynomial fit method but it
requires slightly more manipulation of the data. First the slopes of the chords
shown in Figure 4.5a are calculated and tabulated. Then the chord slopes are
graphed versus time as shown in in Figure 4.5b and these points are fit to a
second-order polynomial.
When t = 0, the slope of the chord, which approximates the tangent, is 1.31 ×
10−9 molal/sec, which is the initial rate.
68 Chemical reactors
(a) 8× 10–6
6× 10–6
concentration, molal
4× 10–6
2× 10–6
0
0 4000 8000 12000 16000
time, sec
1.0× 10–9
5.0× 10–10
0
0 4000 8000 12000 16000
time, sec
mol
rgen = ( min − mout )Q = (10.0 × 10 −4 − 1.86 × 10 −4 ) (2.70 × 10 −5 ) = 2.20 × 10 −8
sec
The flux of Fe3+ to the surface of the pyrite was determined by dividing this
rate by the total surface area of the pyrite.
mol
2.20 × 10 −8
J Fe3+
r 3+
= Fe = sec = 2.34 × 10 −7 mol (4.30)
A m2 m 2sec
(2 g ) 0.047
g
70 Chemical reactors
Vφ
Δt = (4.31)
Q
Δm
= 6.25 × 10 −8 t 2 − 2.23 × 10 −6 t + 2.42 × 10 −5 (4.32)
Δt
Plug flow reactor experiments 71
Table 4.2. Data for Fe3+ reaction with pyrite plug flow experiment.
Q mout ∆t ∆m ∆m/∆t
L/sec molal sec molal molal/sec
2.5×10–5
2.0×10–5
1.5×10–5
∆m/∆t
1.0×10–5
5.0×10–6
0
0 5 10 15 20 25
time, sec
Figure 4.6. Slope of the chords versus time for the Fe3+ reaction with pyrite in a
plug flow reactor experiment. The intercept of this graph is the initial rate.
r = k ∏ mini (4.33)
This equation is nonlinear; direct fitting requires the use of nonlinear regres-
sion, which can be a daunting task. This problem can be avoided by log-
transforming Eq. (4.33) to a linear form.
Fitting the data to this equation using linear regression is a simple way to
find the rate constant and the apparent reaction orders. This equation can
be further generalized by transforming the log k term using the Arrhenius
equation (see Chapter 5). The activation energy (Ea, kJ/mol) term in this
equation expresses how the rate changes with changing temperature (T, K).
E 1
− a
k = Ae R T
(4.35)
A is a pre-exponential factor with the same units as the rate constant and R
is the gas constant (8.314 J/mol K). Equation (4.35) can be log-transformed
to express log k in terms of 1/T.
Ea 1
log k = log A − (4.36)
2.303R T
Combining Eqs (4.34) and (4.36) produces a function that relates the
observed rate to temperature and solution composition.
Ea 1
2.303R T ∑
log r = log A − + n log mi (4.37)
Multiple linear regressions must be used to fit data to Eq. (4.37). See for
example Rimstidt et al. (2012).
Non-ideal chemical reactors 73
Chaïrat et al. (2007) measured the rate of calcium release from dissolving fluora-
patite over the pH range of 3 to 7. The predominant reaction produces calcium
and fluoride ions and undissociated phosphoric acid.
Ca5(PO4)3F + 9 H+ = 5 Ca2+ + 3 H3PO4 + F− (4.38)
Figure 4.7 shows that the logarithm of the Ca2+ release flux versus pH is linear
over this pH range so these data can be fit to Eq. (4.34).
This equation can be transformed to give the Ca2+ flux as a function of the activ-
ity of H+.
(0.020)
J = 10 −9.91(0.09) aH0.88
+ (4.40)
–7
–8
log J
–9
–10
–11
3 4 5 6 7
pH
sample n
n= ∑ ( Δm) M (4.41)
sample 1
Dividing n at each sample time by the initial mass of solution in the reactor
corrects the concentration to what it would have been if no solution had
been withdrawn by previous samples. This corrects the concentration versus
time data to match that expected for an ideal batch reactor.
Alternatively, using the chord or polynomial fit method to find the initial
rate approximately corrects for the changing A/M ratio in a sampled batch
reactor provided the sample size is small relative to the total amount of solu-
tion in the reactor.
Mixing in reactors
Modeling incomplete mixing in chemical reactors can be quite challenging.
Most models use a distribution function of the residence time for fluid elem-
ents in the reactor (Danckwerts, 1953; Denbigh and Turner, 1984). Nauman
(2008) provides a concise history and a useful review of this theory.
Many reaction-transport models treat parcels of fluid as if they are slugs
flowing in a plug flow reactor. In real systems, hydrodynamic mixing and
diffusion (Aris, 1956; Taylor, 1953) cause dispersion of species among the
slugs. Gelhar et al. (1992) review the factors that affect dispersion in aquifer
systems.
Non-steady state mixed flow reactors 75
The (rout − rin) term is the species concentration in the feed solution (min, mol/
kg) minus its concentration in the effluent (mout, mol/kg) multiplied by the
flow rate of the effluent stream (Q, kg/sec).
The rate of accumulation of the species in the reactor is the time rate of
change of its concentration multiplied by the mass of solution in the reactor
(M, kg).
dm
racc = out M (4.44)
dt
Equations (4.42), (4.43), and (4.44) are combined to find the rate of reaction
for non-steady state conditions.
dm
rgen = ( mout − min )Q + out M (4.45)
dt
The first term in this equation is simple to evaluate using the concentration
and flow rate data that are typically collected for MFR experiments. The
second term is evaluated by fitting the effluent concentration to a power law
function and then taking its derivative.
mout = at b (4.46)
dmout
= abt
b −1
(4.47)
dt
The a and b parameters are then used to find the rate of change in the efflu-
ent concentration. The rate of generation at each sample time is found by
combining Eqs (4.45) and (4.47).
Equation (4.44) shows that the larger the mass of solution in the reactor
the larger the rate of accumulation relative to the rate of generation and the
longer it takes for the reactor to reach steady state. This means that a reactor
designed to contain a very small amount of solution will reach steady state
relatively quickly.
Weissbart and Rimstidt (2000) measured the concentration of Ca2+ in the efflu-
ent of a mixed flow reactor in which wollastonite was dissolving. In their experi-
ments, the Ca2+ concentration in the effluent decreased continuously over the
Table 4.3. Calculation of non-steady state rates of Ca2+ release from dissolving wollastonite. The
mass of solution in the reactor 0.047 kg. Data fromWeissbart and Rimstidt (2000).
U1 4.25 × 10−5 1632 1.43 × 10−3 6.08 × 10−8 −8.65 × 10−8 −2.57 × 10−8
U2 2.47 × 10−5 15960 8.13 × 10−4 2.01 × 10−8 −1.02 × 10−9 1.91 × 10−8
U3 1.70 × 10−5 83220 2.16 × 10−4 3.67 × 10−9 −4.08 × 10−11 3.63 × 10−9
U4 1.17 × 10−5 97200 1.99 × 10−4 2.32 × 10−9 −3.02 × 10−11 2.30 × 10−9
U5 4.02 × 10−5 99600 3.40 × 10−5 1.37 × 10−9 −2.88 × 10−11 1.34 × 10−9
U6 3.10 × 10−5 159900 5.79 × 10−5 1.79 × 10−9 −1.14 × 10−11 1.78 × 10−9
U7 9.67 × 10−5 161700 2.20 × 10−5 2.13 × 10−9 −1.12 × 10−11 2.12 × 10−9
U8 1.05 × 10−4 177060 1.85 × 10−5 1.95 × 10−9 −9.38 × 10−12 1.93 × 10−9
U9 8.43 × 10−5 346020 1.94 × 10−5 1.64 × 10−9 −2.54 × 10−12 1.63 × 10−9
U10 7.95 × 10−5 509700 1.94 × 10−5 1.54 × 10−9 −1.20 × 10−12 1.54 × 10−9
U11 2.25 × 10−4 600000 4.80 × 10−6 1.08 × 10−9 −8.70 × 10−13 1.08 × 10−9
course of the experiment and appears to not reach a steady state. The Ca2+ flux
from the dissolving wollastonite can be found using Eq. (4.48) and the data in
Table 4.3.
The (Q)(mCa) column of Table 4.3 is the rate that the effluent removes Ca from
the reactor. It is the first term of Eq. (4.48). The second term in Eq. (4.48) is
found by fitting the concentration of Ca2+ in the effluent versus time to a power
law (Figure 4.8a).
(a) 10×10–4
8×10–4
6×10–4
mCa, mol/kg
4×10–4
2×10–4
0
0 200000 400000 600000
time, sec
1.5 ×10–8
rCa, mol/sec
1.0 ×10–8
5.0 ×10–9
0
0 200000 400000 600000
time, sec
The time rate of change of Ca concentration in the reactor solution is the deriva-
tive of Eq. (4.49).
dmCa
= (3.52 ) ( −0.948) t = −3.34t −1.948
−1.948
(4.50)
dt
Equation (4.50) is evaluated for each sample time and the (dmCa/dt)(M) column
reports the dissolved Ca that is retained in the reactor causing the Ca concentra-
tion to increase. Adding the (Q)(mCa) and (dmCa/dt)(M) columns gives the rate of
calcium release from the wollastonite at each sampling time (Figure 4.8b).
Minerals usually dissolve incongruently during the initial stages of the dis-
solution process. In this case, Ca is released from wollastonite faster than Si. In
addition, the rate of both Ca and Si declines with increasing extent of reaction.
Eventually the difference between release rates, as well as the change in the rates,
from one sample time to the next becomes smaller than the resolution of the
rate measurements. When both of these conditions exist the rates are assumed
to represent the dissolution rate that is appropriate for long-term geochemical
processes.
Chapter 5
Molecular kinetics
Rate equations are quantitative models of the time course of chemical reac-
tions. Although rate equations are based on macroscopic observations, they
reflect processes that occur at the molecular scale. This chapter reviews some
of the important models that link these two scales. These models are espe-
cially useful because they constrain the mathematical form of rate equa-
tions and they provide a conceptual basis for thinking about the reactions.
Because water is so important in geochemical systems, this chapter focuses
on models for reaction rates in the aqueous phase.
Transition-state theory
Chemical reactions break and reform bonds within and between molecules
so that the bonding arrangement of the reactants gives way to the bond-
ing arrangement of the products. For a reaction to occur: (1) the molecules
must “collide” with each other to form a cluster; (2) the atoms in that cluster
must be configured in the approximate geometry of the products; and (3)
the cluster must contain sufficient energy to allow the rearrangement of the
electrons from the breaking bonds to the developing bonds.
The idea that molecules must collide to form a cluster correctly predicts
that the reaction rate increases with concentration (or pressure) because
increasing the number of molecules per unit volume will result in more col-
lisions. This is the reason that rate equations stipulate that rates increase
with increasing concentration, pressure, or activity.
The Franck–Condon principle states that the rearrangement of electrons
to form new bonds occurs so quickly that the nuclei of the reacting atoms do
not change position during the bond breaking–bond forming process. This
means that for a reaction to occur, the reacting species must collide in such a
way that they form a cluster with the atoms arranged with the geometry of
the product molecules. This requirement means that most collisions do not
lead to a reaction. A reactive cluster that meets these geometric constraints
is called an activated complex.
79
80 Molecular kinetics
cool
Quasi-equilibrium model
The quasi-equilibrium model considers a simple reaction where atom A
reacts with molecule BC and displaces atom C to produce molecule AB.
A + BC = AB + C (5.1)
For this reaction to occur, molecule A must approach molecule BC from the
correct direction with enough translational energy to form a transition-state
cluster.
A + BC = ABC* (5.2)
aABC*
K* = (5.3)
aA aBC
The activity of the transition-state cluster, ABC*, is found from the product
of the equilibrium constant and the activity of each reactant.
1 *
mABC* = K aA aBC (5.5)
γ ABC*
k T k T K*
R = B mABC * = B a a (5.6)
h h γ ABC* A BC
82 Molecular kinetics
(a)
A B- C Reactants
BC distance
A
AB C distance
(b)
A-B-C
Ea
Energy
The Boltzmann constant (kB) has a value of 1.381 × 10−23 J/K and the Plank
constant (h) is 6.626 × 10−34 J sec. This equation is further modified by intro-
ducing a transmission coefficient (κ, no units) to account for the fact that a
few of the transition-state clusters fail to form products and instead revert
Quasi-equilibrium model 83
back to being reactants. The transmission coefficient often has a value near
one and is usually ignored.
κ kBT *
R= K aA aBC (5.7)
γ ABC* h
Unit analysis of this equation shows that it has units of molal/sec. The first
three terms are combined into a rate constant.
κ kBT *
k= K (5.8)
γ ABC* h
R = −100 mFeS
2
+
2 O3
(5.11)
R 1 × 10 −4 mol
m (FeS2O3 )22+ * = = = 1.61 × 10 −19 (5.12)
( B )
k T h 6 .21 × 10 12
kgg
This concentration is well below detection limits for aqueous species and any
[(FeS2O3)22+]* cluster that forms has a lifetime of ~1.6 × 10−13 sec, which makes
this exercise entirely hypothetical but hopefully instructive.
84 Molecular kinetics
Because ∆Gr = ∆Hr − T∆Sr, Eq. (5.13) can be expanded to relate the enthalpy
and entropy of activation to K*.
∆H r* − T ∆Sr* = − RT ln K * (5.14)
Combining Eqs (5.8) and (5.14) gives a relationship between the enthalpy
and entropy of activation and the rate constant (k).
− ∆H r* 1 ∆Sr* κ kBT
ln k = + + ln (5.15)
R T R γ ABC* h
For an ideal gas, the activation energy (Ea kJ/mol) is related to ∆H*.
∆H * = ∆ ( E + PV ) = ∆E a + ∆nRT (5.16)
∆S * k T
ln A = r + ln B (5.17)
R h
(a) 30
25
20
∆T, °C
10 °
15
0°
75 0°
5 °
25
10 Q10
0
20 40 60 80 100 120
Ea, kJ/mol
(b) 5
− Ea 1
ln k = + ln A (5.18)
R T
− Ea 1
log k = + log A (5.19)
2.303R T
86 Molecular kinetics
k − Ea 1 1
log 2 = − (5.20)
k1 2.303R T2 T1
This equation is the basis for the Q10 rule of thumb, which states that for
many reactions the rate doubles for every 10°C increase in temperature. This
approximation seems to work well for biochemical reactions near room
temperature but is somewhat problematic for geochemical reactions that
occur over a wider range of temperatures and have a wider range of activa-
tion energies. The effect of temperature and activation energy on doubling
is shown in Figure 5.3a. Figure 5.3b shows how many times the rate will
increase for every 10°C increase in temperature. These figures show that the
Q10 rule applies for reactions that have an activation energy near 50 kJ/
mol and occur near 25°C. Note that for a reaction with a low activation
energy, the temperature must increase by a large amount to double the rate
(Figure 5.3a); and for a reaction with a high activation energy only a small
temperature increase will double the rate.
Williamson and Rimstidt (1993) determined the rate of the ferric thiosulfate
reaction at temperatures ranging from 12° to 34°C. When those data are plotted
on a graph of log k versus 1/T (Figure 5.4) the equation for the best-fit straight
line is
log k = 23.9 − 6426/T
The activation energy for this reaction can be calculated from the slope of
this line.
The relatively large value of this activation energy is explained by the fact that
the two positively charged ions must have enough energy to overcome their elec-
trostatic repulsion. A postulated geometry for the activated complex is shown in
Figure 5.5. Note that the reactants approach each other in a way that brings the
negatively charged thiosulfate adjacent to the positively charged ferric iron ion
of the other reactant molecule and vice versa. This geometry keeps the highly
charged ferric iron ions separated as far as possible.
The very large activation energy of the ferric thiosulfate reaction means that
a 5°C temperature increase, from 25° to ~30°C, will cause the rate to double
(Figure 5.3).
Compensation effect 87
• •
• •
2.5
•
log k
2
•
1.5
•
1
0.0032 0.0033 0.0034 0.0035 0.0036
1/T, K–1
Figure 5.4. Arrhenius plot for the ferric thiosulfate reaction. The slope of the log k
versus 1/T graph is 6426 so the activation energy for the reaction is 123 kJ/mol.
Compensation effect
Several studies have reported that for a series of related reactions the pre-
exponential in the Arrhenius equation correlates with the activation energy.
This correlation has a physical significance for homogeneous, gas phase
reactions but there is no theoretical justification for reactions involving con-
densed phases (Garn, 1975). Compensation effect models should be viewed
with caution because errors in real data are known to cause a correlation
between the slope and intercept of regression models (Liu and Guo, 2001;
88 Molecular kinetics
Rimstidt et al., 2012). Comparing the log k values at 298 K to the activation
energy can minimize the correlation due to these errors (Ohlin et al., 2010).
At infinite dilution (I = 0) where the activity coefficient for the activated
complex equals one (and assuming that κ = 1) the rate constant for this
reaction is related to the equilibrium constant for the formation of the acti-
vated complex.
k T
ko = B K * (5.22)
h
γ γ k T γ γ
k = A B B K * = A B ko (5.23)
γ AB* h γ AB*
The Debye–Hückel limiting law relates the activity coefficients of all the spe-
cies, including the activated complex, to the ionic strength.
− log γ i = AZ 2 I (5.24)
Zi is the charge on the ion; I is the ionic strength of the solution; and A =
0.5092 (at 25°C). Equations (5.22) and (5.23) can be combined to create a
relationship between log k and ionic strength.
( 2
)
log k = log ko + ZA2 + ZB2 − (ZA + ZB ) A I (5.25)
This equation predicts that a graph of log k versus √I will be a straight line
with a slope of 2AZAZB. At 25°C, 2A = 1.02 ≈ 1, so the slope of the line is
Effect of ionic strength on rates 89
approximately equal to the product of the charges on the A and B ions. This
relationship is not only an effective way to predict the effect of ionic strength
on rates but it can also be used to confirm the stoichiometry of postulated
transition states.
In the previous two examples, the activated complex for the ferric thiosulfate
reaction was postulated to be [(FeS2O3)22+]*. This molecule is a combination of
two FeS2O3+ ions, each with a +1 charge, so the product ZA ZB is +1. If this pos-
tulated stoichiometry for the activated complex is correct, the slope of a graph
of log k versus √I (Figure 5.6) should be +1. The Williamson and Rimstidt
(1993) data for the rate as a function of ionic strength produces a straight line
with a slope of 0.928.
Multiplying the slope of the line by 1.02 gives a value of 0.945 for the product of
the charges on the reacting ions, which is very close to the expected value of +1.
Note that this equation can be used to predict that log ko = 2.17. This test gives
further support to the postulated composition of the activated complex.
2.6
•
• •
2.5 27°C •
•
2.4
log k
2.3
2.2
2.1
2
0 0.1 0.2 0.3 0.4 0.5
I
Figure 5.6. The graph of log k versus √I for the ferric thiosulfate reaction. This
graph has a slope of 0.928 (≈ 1), which is consistent with each of the interacting
molecules having a charge of +1.
90 Molecular kinetics
∆V * = VAB* − VA − VB (5.27)
van Dldik et al. (1989) have compiled a large amount of data pertaining to
the effect of pressure on reactions in aqueous media. For the reaction, A +
B = AB* → products, the rate constant change with pressure is proportional
to the volume of activation.
∂ ln k ∆V * ∂ log k ∆V *
= − or
∂P T =− (5.28)
∂P T RT 2.303RT
Two kinds of volume of change occur when dissolved species react to form
an activated complex. The volume occupied by the species themselves
changes and the associated water molecules undergo rearrangements that
change their volume.
When Eq. (5.28) is integrated between 1 bar where the rate constant is
k1 and P where the rate constant is k2, log (k2/k1) is predicted to be a linear
function of pressure.
k ∆V *
log 2 = − P (5.29)
k1 2.303RT
Actual log k versus P relationships usually display some curvature over large
ranges of pressure so the log k values are often fitted to a parabolic function
of pressure.
log k = a + bP + cP 2 (5.30)
The activation volume is calculated from the slope of this fit at P ≈ 0 (usually
at 0.1 MPa)
∂ log k
=b (5.31)
∂P T
A warning is appropriate here. If the rates are expressed using volume units
such as molarity, the compressibility of the solvent must be considered when
Effect of pressure on rates 91
Huss Jr. et al. (1982) measured the rate of Mn(II)-catalyzed oxidation of SO2 by
dissolved oxygen at pressures up to about 140 MPa (1400 bars). A fit of log R
versus P for the data shown in Figure 5.7 gives the equation
d log k
= 2 (7.76 × 10 −7 ) P − 1.95 × 10 −3 (5.33)
dP
m 3MPa 1
∆V * = − (2.303)(298.15K ) 8.314 × 10 −6 −1.95 × 10
−3
K mol MPa (5.34)
m 3
∆V * = 1.11 × 10 −5
mol
–5.4
•
–5.5
•
•
•
log R
–5.6
•
•
–5.7
•
–5.8
0 25 50 75 100 125 150
P, MPa
For this reaction to occur, the water molecules coordinating the ferrous ion
at an average distance of 2.21Å must move closer to the ion because the aver-
age distance of water molecules coordinating a ferric iron ion is 2.05Å. The
electrostatic force that keeps each water molecule near the ion is constantly
perturbed by thermal vibrations that cause each water molecule to move
toward and away from the ion as if it was attached by a spring. According
to Hook’s law, the potential energy associated with each bond is a parabolic
function of the vibrational displacement, x.
1 2
Em = kx (5.36)
2
Figure 5.8 shows the parabolic relationships between the energy (E) of the
donor and acceptor species and the average displacement (x) of the solvent
molecules from their most stable configuration. The intersection of these
two parabolas is the energy barrier to the electron transfer step (i.e. free
energy of activation, ∆G*). The free energy of activation is a function of the
Photochemical reactions 93
∆G * =
( E m + ∆G °)2 (5.37)
4E m
The rate constant for the electron exchange reaction is the product of a fre-
quency term and ∆G*.
kT ( E m + ∆G °)2
k= exp − (5.38)
h 4 E m kT
Photochemical reactions
Solar energy delivers enormous amounts of energy to the Earth’s sur-
face, making photochemical reactions very important in geochemical and
94 Molecular kinetics
biological processes. The global average solar energy incident at the top of
the atmosphere is 342 W/m2. Some 102 W/m2 is reflected back into space
and between 65 and 98 W/m2 is absorbed by the atmosphere (Arking, 1996).
That leaves between 142 and 175 W/m2 as the average solar energy flux at
the Earth’s surface, so that over a year sunlight delivers between 4.5 and
5.5 million kJ/m2 of energy to the Earth’s surface. Much of that radiant
energy is converted to heat but a significant fraction drives photochemical
reactions, the most important of which is photosynthesis. Many other pho-
tochemical reactions are important to geochemistry. For example, photo-
oxidation of Fe2+ (Southworth, 1995) and Mn2+ (Nico et al., 2002) exert an
important control on the redox chemistry of stream water. Most of those
photochemical reactions occur in surface waters and in films of water on
mineral surfaces.
Planck’s law states that the amount of energy that each photon carries
depends upon its frequency (ν, sec−1), multiplied by Planck’s constant (6.63
× 10−34 J/sec) and the frequency depends on the speed of light (c, 2.99 × 108
m/sec) and the wavelength (λ, m).
hc
E ( J/photon ) = hν = (5.39)
λ
This means that light with shorter wavelengths delivers a larger amount of
energy per photon. Blue green light (λ = 480 to 500 nm), the most intense
sunlight at the Earth’s surface, delivers 240 to 250 kJ per Einstein. An
Einstein is 1 mole of photons (NA = 6.02 × 1023) regardless of wavelength.
N Ahc
E ( J/Einstein ) = (5.40)
λ
Figure 5.9 shows the relative intensity of sunlight and the amount of
energy delivered by one Einstein of light as a function of wavelength.
When a molecule or ion absorbs a photon, that photon’s energy can be
dissipated in several different ways, but one way is for that energy to cause
a chemical reaction to occur. The first law of photochemistry is that a com-
pound must absorb light for a photochemical reaction to occur (Grotthuss–
Draper law). The second law of photochemistry is that each photon that
is absorbed activates only one molecule for a subsequent reaction (Stark–
Einstein law). The quantum yield (φ) is defined as the number of molecules
that react divided by the number of photons absorbed. It can also be defined
in terms of moles.
ultraviolet
infrared
orange
yellow
green
violet
blue
red
400
(b)
energy, kJ/E
300
200
Figure 5.9. (a) Graph of the
relative intensity of solar
radiation at the Earth’s
surface. (b) Graph of the
100 energy contained in 1 mole
800 700 600 500 400 300 of photons as a function of
wavelength, nm wavelength of light.
Radiolysis
Radiolysis reactions are brought about by the absorption of ionizing radi-
ation. The radiation may belong to the light group, i.e. X-rays, γ-rays, elec-
trons–positrons, muons; or the heavy group, i.e. protons, α-particles, fission
fragments. When these particles pass through a medium, some of their
energy is transferred to the molecules of that medium raising their electronic
states. This can lead to neutral dissociation or to ionization if the energy
transfer is sufficiently large.
Radiolysis of water produces hydroxyl free radicals (•OH) (Le Caër, 2011)
along with hydrated electrons and •H atoms. Subsequent reactions produce
H2O2, H2, O2, and, depending upon the solution composition, various other
species. Because the H2 can rapidly diffuse away from the site of genera-
tion, radiolysis reactions often alter the local redox conditions. For example,
96 Molecular kinetics
Free radicals
Free radicals are chemical species that have one or more unpaired electrons.
They can have a positive, negative, or zero charge. Although a few radi-
cals are stable (e.g. O2 and NO) or highly persistent, most are highly reac-
tive. Thermally activated, electron transfer, photochemical, and radiolysis
processes can all produce free radicals. Free radicals are common reaction
intermediates in organic and biochemical reactions and they are no doubt
important in many organic geochemical processes (Dominé et al., 2002).
Free radicals play important roles in many inorganic geochemical proc-
esses. For example, hydroxyl free radicals are produced during pyrite oxida-
tion by O2. Presumably they form as the result of a modified Fenton reaction
(Rimstidt and Vaughan, 2003). Free radicals are further involved as various
sulfoxy species are oxidized to sulfate (Druschel et al., 2003). Oxidation of Fe2+
is another example of a reaction that involves free radical intermediates.
Breaking Si-O bonds by fracturing quartz produces O3–Si• and Si–O• free
radicals that can subsequently react with water to produce •OH radicals
(Narayanasamy and Kubicki, 2005). The •OH species is quite toxic to cells
and may be responsible for the lung scarring that is symptomatic of silicosis.
•OH radicals can also alter DNA and may be responsible for the develop-
ment of lung cancer (Fubini, 1998).
Free radicals are studied by generating them using radiolysis. These stud-
ies have produced extensive tabulations of rate constants for free radical
reaction rates (Bielski et al., 1985; Buxton et al., 1988; Buxton et al., 1995;
Neta et al., 1988).
Chemical bonding and reaction mechanisms 97
Reaction mechanisms
A reaction mechanism is a model that explains, elementary step by elemen-
tary step, how reactants become products. Reaction mechanism models also
explain how the rate of each elementary step is influenced by the structure
of its transition state. Because several different reaction paths may be possi-
ble for a particular reaction, there is no way to a priori determine a reaction
mechanism. Transition states cannot be observed directly because they are
present at very low concentrations and have a fleeting existence. This means
that reaction mechanisms must be inferred using a combination of forward
models based on fundamental physical principles and reverse models cal-
ibrated by the outcomes of experimental tests. Detailed studies of many
kinds of reactions have generated some guiding principles for working out
reaction mechanisms (Ašperger, 2003; Moore and Pearson, 1981). Casey
(2001) and Casey and Swaddle (2003) give some guidelines that are specific
to the dissolution of oxide and hydroxide compounds, which are common
in geochemical environments.
In this reaction, ferrous iron, a Lewis acid, reacts with the lone pair electrons
on water molecules, which are Lewis bases, to form a hydrated cation, a
Lewis adduct. Although the sulfate ion also interacts with water molecules,
this interaction is much weaker than the Lewis acid–base interaction of the
ferrous iron. Many, perhaps all, chemical reactions can be viewed as Lewis
acid–base interactions, making this theory very useful for developing mod-
els of reaction mechanisms.
Organic chemists refer to Lewis acids as electrophiles because they are
attracted to electron-rich sites on other molecules, and they refer to Lewis
bases as nucleophiles because they are attracted to electron-deficient molecu-
lar sites. Ingold (1969) developed a classification of reaction mechanisms
based on this idea and this classification is the foundation for modeling the
mechanisms of organic reactions (Bruckner, 2002; Grossman, 1999). Casey
(2001) and Casey and Swaddle (2003) adapted some of these principles to
apply to the dissolution of oxides. The electron-rich, and therefore nucleo-
philic, sites on molecules are nicely visualized using the electron localization
function (Gibbs et al., 2005).
Donor–acceptor theory (Gutmann, 1978) further extends the Lewis bond-
ing theory to explain how bonds in reacting molecules lengthen (become
weaker) or shorten (become stronger) as a result of the Lewis base (electron
donor) and Lewis acid (electron acceptor) interactions. Gutmann proposed
three rules for predicting bond length variations during donor–acceptor
interactions:
1. As donor and acceptor atoms approach each other, the induced rearrangement
of electrons in both the donor and the acceptor causes the adjacent bonds in each
to lengthen.
2. Donor–acceptor interactions cause electron density shifts throughout the entire
molecule and those shifts produce changes in bond lengths.
3. As the coordination number of an atom in an adduct increases so do the bond
lengths of any bonds originating from the coordination center.
Chemical bonding and reaction mechanisms 99
The donor–acceptor rules are consistent with quantum mechanical models that
show how electrons in a bond respond as a donor molecule approaches the elec-
trophilic site on an acceptor molecule. A well-known quantum mechanics model
by Lasaga and Gibbs (1990) illustrates the attack of a H2O molecule on a Si–O
bond in a H6Si2O5 molecule. Their model shows that the electron pair on the
water molecule’s oxygen atom mounts a nucleophilic attack on one of the Si
atoms leading to a 5-coordinated Si atom. At the same time, one of the water
molecule’s hydrogen atoms mounts an electrophilic attack on the bridging oxy-
gen in the H6Si2O5 molecule. The simultaneous alignment of the oxygen end of
the water molecule with the Si atom and the electrophilic attack of the hydrogen
atom on the bridging oxygen atom creates the geometry of the transition state
(Figure 5.10). Once that transition state has formed, the curly arrows show the
shifts in electron density from preexisting bonds to the developing bonds. This
shift divides the H6Si2O5 molecule into two H4SiO4 molecules.
(a) (b)
(c)
Figure 5.10. Schematic illustration showing how nucleophilic and electrophilic interactions cause (a) a water
molecule to align with a Si–O bond to form (b) a transition state that results in electron rearrangement,
which disassociates the Si–O bond and forms (c) two Si–OH moieties.
100 Molecular kinetics
11. Elementary reactions must be feasible from a bond energy and geometry stand-
point. Extensive bond and atom rearrangement does not occur in a single
elementary step.
12. The charge on the transition state of the rate-determining step can be deduced
from the effect of ionic strength on the rate.
The Cu+ is rapidly re-oxidized by a very fast reaction with ferric iron, which
converts it back to Cu2+.
Cu+ + Fe3+ = Cu2+ + Fe2+ (5.43)
This meets the definition of a catalyst, which is a species or site that lowers the
activation energy of the reaction and is not consumed by the reaction.
The ferric thiosulfate reaction can be inhibited by the addition of Ag+, which
forms a very strong complex with thiosulfate and that complex is unreactive
because Ag+ cannot be reduced to Ag0 by reaction with thiosulfate or oxidized
by reaction with Fe3+.
Chapter 6
Surface kinetics
dG
σ = (6.1)
dA T , P ,etc.
102
Geometric surface area models 103
The surface free energy term is always positive; this means that the equi-
librium surface should acquire a geometry that minimizes its surface free
energy. An atomically smooth sphere has the lowest overall surface free
energy for solids with isotropic bonding. For example, amorphous silica
tends to develop a spherical form. For crystalline solids, the faces with the
highest atomic packing densities and shortest bond lengths have the lowest
σ values. These “low index” faces tend to occur more frequently in crystal
forms, although the growing crystals are not at equilibrium with the sur-
rounding solution so most crystals actually display growth forms rather
than true equilibrium forms. A crystal that is at perfect equilibrium with
the surrounding solution would have atomically smooth faces, but actual
crystals usually display growth hummocks and other surface irregularities
related to the growth process.
The value of σ for a crystal face relative to the other crystal faces can be
found using Wulff’s theorem, which states that there exists an interior point
such that the surface free energy of a crystal face is proportional to the perpen-
dicular d
istance from that point to the crystal face (Wolff and Gualtieri, 1962).
Dmax − Dmin
De = (6.2)
Dmax
ln
Dmin
104 Surface kinetics
The effective diameter given by this model is very close to the arithmetic
average of the maximum and minimum grain diameter for a narrow grain
size range, but for a wider size range it is lower because the smaller particles
contribute more surface area.
The simplest model that relates the specific surface area of grains to their
diameter assumes that the grains are spherical. This model expresses the
specific geometric surface area (Ageo, m2/g) as a function of the molar vol-
ume (Vm, m3/mol) of the substance, the diameter of the grains (De, m), and
the molecular weight of the substance (Wm, g/mol).
6Vm m 2
Ageo = (6.3)
DWm g
Equation (6.3) can be transformed to give the specific surface area as a func-
tion of the density (ρ, g/cm3) and diameter (D, m) of the grains.
6 × 10 −6 m 2
Ageo = (6.4)
Dρ g
If a sample of quartz sand is sieved to recover the 1 to 2 mm diameter size frac-
tion, the effective diameter of the grains in the sample can be found using Eq.
(6.2).
2 −1
De = = 1.44 mm = 1.44 × 10 −3 m
2
ln
1
The density of quartz is 2.66 g/cm3, so the geometric specific surface area of the
sample is found using Eq. (6.4).
6 × 10 −6 m2
Ageo = = 0.00157
(1.44 × 10 )(2.66)
−3
g
2
V 3
Afu = m (6.5)
NA
The area occupied by 1 mole of formula units (Amol, m2/mol) is the area per
formula unit multiplied by Avogadro’s number.
1
Dfu = (6.7)
Amol
1
t fu = (6.8)
JAmol
Forsterite (Mg2SiO4) has a molar volume of 4.365 × 10−5 m3/mol. The average
area occupied by one formula unit of forsterite is found using Eq. (6.5).
2
m3 3
4.365 × 10 −5
mol m2
Afu = = 1.738 × 10 −19 (6.9)
fu fu
6.023 × 10 23
mol
The surface area occupied by 1 mole of forsterite is found using Eq. (6.6).
m2 fu m2
Amol = 1.738 × 10 −19 6.023 × 10 23 = 1.047 × 105 (6.10)
fu mol mol
1 mol
Dfu = = 9.553 × 10 −6 2 (6.11)
1.047 × 105 m
1
t fu = = 480 sec ≈ 8 min (6.12)
mol m2
1.99 × 10 −8
1.047 × 10 5
m 2sec mol
The porosity of packed beds ranges from about 50% for loosely packed grains
to about 25% for close packed grains. The A/M for a packed bed of 1 mm diam-
eter grains with 30% porosity filled with a solution with a density of 1 g/cm3 is
about 14 m2/kg.
A 1− φ 6 0.7 6 m2
M = φ 1000 ρD = 0.3 (1000 )(0.001)
= 14
kg
Figure 6.1 shows how the A/M ratios for packed beds filled with a solution with a
density near 1 g/cm3 varies with grain diameter and porosity. Real grains are not
perfect spheres so these values should be increased by 5 to 10 times to account
for their surface roughness.
Reactions at surfaces
Surface reactions comprise four steps.
1. Adsorption of a fluid phase species (A) to a surface site (S): A(aq or v) +
S = A…S
2. Surface diffusion of the reacting species to a reactive site (S#): A…S → A…S#
3. Reaction between the reactive site and the adsorbed species to produce a product
species (B): A…S# = A–S#→ B–S
Adsorption (Langmuir model) 107
Table 6.1. A/V and A/M ratios for some typical solid/solution geometries.
The A/V ratios are used for models where the solution concentrations are
expressed in volume units. The A/M ratios are used in models where the
solution concentrations are expressed in mass units. The units for density, ρ,
are g/cm3.
Cylindrical pipe A = 2R and A = 2
radius (R, m) V M 1000 ρR
Rectangular pipe A = 2 ( H + W ) and A = 2 ( H + W )
height (H, m), width (W, m) V HW M 1000 ρHW
Fracture with parallel walls A = 2 and A = 2
width (W, m) V W M 1000 ρW
Packed bed of spheres
A = 1 − φ 6 and A = 1 − φ 6
diameter (D, m), porosity (φ) V φ D M φ 1000 ρD
(a) 1500
1000
A/M, m2/kg
500
= 0.25
= 0. 5
0
0 0.0005 0.001
grain diameter, m
(b) 3000
1000
= 0.25
A/M, m2/kg
100
= 0.5
Figure 6.1. (a) A/M ratio
for packed beds filled with
a solution with a density
of 1.0 g/cm3 as a function 10
of grain diameter and bed
porosity. (b) Same data
presented on a graph with 3
logarithmic scales. 0.00001 0.0001 0.001
grain diameter, m
dθ
= k+ mΓ (1 − θ ) − k− Γθ (6.13)
dt
k+ mΓ (1 − θ ) = k− Γθ (6.14)
(6.15)
k+ Γ θ
m =
k− Γ 1 − θ
θ
Kbm = (6.16)
1 − θ
Kbm
θ= (6.17)
1 + Kbm
Kbm
Γa = Γ (6.18)
1 + Kbm
1 1 1
= + (6.19)
Γ a Γ ΓK b m
A graph of 1/Γa versus 1/m has a slope of 1/ΓKb and an intercept of 1/Γ.
This fit gives excessive weight to data in the low concentration range and is
sensitive to data error. The Scatchard equation gives more weight to data in
the high concentration range.
Γa
= ΓK b + K b Γ a (6.20)
m
m m 1
= + (6.21)
Γ a Γ ΓK b
1× 10–5
30
8× 10–6
10
6× 10–6
, mol/m2 3
4× 10–6
a
Distribution coefficient
The distribution coefficient is a simplification of the Langmuir equation,
Eq. (6.16). For small amounts of surface coverage, the 1 − θ term in that
equation is approximately equal to 1, so it simplifies to a linear relationship
between the fractional surface coverage and the solution concentration.
Kbm = θ (6.22)
Γ a = Γθ = ΓK b m (6.23)
nsp
KD = (6.25)
m
Surface catalysis: unimolecular decomposition 111
K b mΓAsp
KD = = K b ΓAsp (6.26)
m
In this derivation, Γa and Γt are multiplied by the specific surface area of the
solid in order to convert the concentration units from mol/m2 to mol/g-solid.
Defining distribution coefficients this way makes them functions of the spe-
cific surface area so KD values vary with grain size.
BET surface area
Reaction rates between fluid species and solid surfaces are often normalized
using a surface area value determined by the BET (Brunauer–Emmett–Teller)
method. The theory behind this method was developed by Brunauer et al.
(1938) and is thoroughly described and evaluated in Lowell and Shields
(1991). This model is based on the idea that as the pressure and temperature
of a gas approaches the pressure and temperature where liquid and vapor are
in equilibrium, the molecules will deposit on the surface as a thin, multilayer
film. In most cases N2 is used as the gas and a powdered sample is placed into
an evacuated chamber that is immersed into boiling liquid nitrogen to fix the
temperature. Increasing amounts of N2 gas are introduced into the chamber
and the amount adsorbed at pressures ranging from 5 to 35 kPa is measured.
Usually three to five data points are determined. These data are fit to the BET
isotherm and the amount of gas that comprises a monolayer on the solid sur-
face is computed from the slope and intercept of the isotherm. The surface
area is then calculated assuming a sorption cross-section for an adsorbed N2
molecule. Prior to the BET surface area determination, the surface of the sam-
ple is usually “cleaned” by heating in a vacuum, often to temperatures near
300°C. This treatment has little effect on most minerals but can drastically
alter hydrous phases and sulfide minerals that decompose releasing water or
sulfur vapor into the vacuum system. BET-determined surface areas are usu-
ally five to ten times larger than surface areas calculated using Eq. (6.3). For
example, White and Peterson (1990) report that the average BET surface area
for a relatively large number of samples is about seven times the surface area
estimated from the geometric surface area estimated using Eq. (6.4).
A + S←
→ AS k
k+1
+2
→ products (6.27)
k−1
After a short time, the rate that A adsorbs to the surface equals the rate
that AS breaks down to products. The rate of disappearance of A can be
expressed in terms of the rate of conversion of AS to products.
dmA
= k+2 Γ AS = k+2θΓ (6.28)
dt
Once the reaction reaches steady state, the rate of change of the AS concen-
tration is zero.
d Γ AS
= 0 = k+1AΓ (1 − θ ) − k−1Γθ − k+2 Γ (6.29)
dt
This equation can be rearranged to find the fraction of surface sites covered
by adsorbed A.
k+1mA
θ= (6.30)
k+1mA + k−1 + k+2
Equation (6.30) can be combined with Eq. (6.28) to find the rate of
consumption of A.
If k+2 << k+1mA and k+2 << k−1, the fraction of surface coverage is described
by the Langmuir isotherm.
k+1mA
θ= (6.33)
k+1mA + k−1
This means that the rate depends upon the concentration of AS.
dmA
= K b k+2 ΓmA (6.35)
dt
When mA is large, so the fraction of surface coverage is large, the rate is zer-
oth order in A.
dmA
= k+2 Γ (6.36)
dt
In a detailed study of the dissolution rate of synthetic UO2, de Pablo et al. (1999)
found that instead of being linear, as is typical of many mineral dissolution reac-
tions, the log J versus log MHCO3 graph shows a distinct curvature. They inter-
preted this to mean that bicarbonate ions sorb to a surface site as a dissolution
reaction step and when bicarbonate concentrations are high these sites become
nearly saturated with adsorbed bicarbonate. This means that the reaction order
transitions from first order at low bicarbonate concentration to zeroth order at
high bicarbonate concentration. To model this behavior, they fit the rate data to
–6
–6
–8
–8
–10
log J
–10
–12
Figure 6.3. Graph of the
0
–12 uranium dissolution flux
–1 from oxidizing UO2 as a
lo –2 function of bicarbonate
gm
H –3 0 and dissolved oxygen
CO –1
3 – –4 –2 concentrations at 25°C,
–5 –3 based on data from de
–4 log m O 2
–5 Pablo et al. (1999).
114 Surface kinetics
an equation that assumes that both bicarbonate and oxygen interact with the
dissolving surface via Langmuir isotherms.
A graph of log J as a function of log M HCO3 and log M O2 (Figure 6.3) shows that
for high values of log M O2 there is a linear relationship between log J and log
M HCO3 , but at low values of log M O2 and high values log M HCO3 , the log J values
become independent of log M HCO3 .
Surface retreat rate
For a dissolving flat surface, the linear rate of surface retreat (dz/dt, m/sec)
is the dissolution flux (J, mol/m2sec) multiplied by the molar volume of the
solid (Vm, m3/mol).
dz
= − JVm (6.38)
dt
The rate of forsterite dissolution in Example 6.2 is 1.99 × 10−8 mol/m2sec and the
molar volume of forsterite is 4.365 × 10−5 m3/mol. At pH 3 and 25°C the surface
retreat rate is calculated using Eq. (6.38).
dz mol m3 m
= − 1.99 × 10 −8 2 4.365 × 10 −5 = −8.69 × 10
−13
(6.39)
dt m sec mol sec
Shrinking particle model
Many hydrometallurgy and chemical engineering processes involve the dis-
solution of solid particles. These processes are typically modeled using a
shrinking particle model, which accounts for the change in dissolution rate
due to decreasing grain surface area as they dissolve away (Burkin, 2001;
Levenspiel, 1972a). The classical derivation of the shrinking particle rate
equation given here shows how the dissolution rate constant, k+(mol/m2sec),
can be calculated from the particle rate constant (kp) that is derived from a fit
of the extent of reaction versus time. The shrinking particle model assumes
Shrinking particle model 115
that the dissolution flux (J = k+, mol/m2sec) is constant over the full extent
of the reaction, either because the solution composition is constant, making
the reaction pseudo-zeroth order, or because the reaction is actually zeroth
order. The first step is to express the rate of reaction of a spherical particle
(r, mol/sec) in terms of the particle’s radius (R, m).
The rate of volume loss for the particle (dV/dt, m3/sec) is the rate of reaction
multiplied by the molar volume (Vm, m3/mol).
dV
= rVm = −4π R 2 k+Vm (6.41)
dt
The rate of volume loss of the particle is also given by the time derivative of
the volume formula for a sphere.
dV dR
= 4π R 2 (6.42)
dt dt
Setting Eq. (6.41) equal to Eq. (6.42) gives the rate of change of the radius.
dR
= − k+Vm (6.43)
dt
The fraction of material that is reacted away (α) at any time can be expressed
in terms of the radius at that time (R, m) and the initial radius (Ro, m).
R3
α = 1− 3 (6.44)
Ro
R = Ro (1 − α ) 3
1
(6.45)
The time derivative of Eq. (6.45) gives the rate of change of the fraction
reacted in terms of the rate of change in the particle radius.
dα R 2 dR
= −3 3 (6.46)
dt Ro dt
d α 3k+Vm
(1 − α ) 3
2
= (6.47)
dt Ro
116 Surface kinetics
k+Vm
kp = (6.48)
Ro
dα
= 3k p (1 − α ) 3
2
(6.49)
dt
1 − (1 − α )
1
3 = k pt (6.51)
A graph of 1 − (1 − α ) 3 versus t has the slope of kp. The rate constant (k+) for
1
Ro k p
k+ = (6.52)
Vm
Van Herk et al. (1989) proposed using forsterite (Mg2SiO4) to neutralize waste
industrial acids. In order to provide rate data to support their idea, they per-
formed a series of laboratory tests that determined the fraction of forsterite
reacted away (α) by contact with acidic solutions for various times. Their data
can be used, along with Eqs (6.51) and (6.52), to find the rate constant for for-
sterite dissolution.
The experiment dissolved 105–125 µm forsterite grains in pH 1 HCl solution
at 40°C. The results are shown in Figure 6.4. When 1 − (1 − α)⅓ is graphed versus
t, the slope of the resulting straight line is 1.23 × 10−6 sec−1. The effective diam
eter of the particles can be estimated using Eq. (6.2).
(a) 1
0.8
0.6
0.4
0.2
0
0 2000 4000 6000
time, min
(b) 0.5
0.4
0.3
1 – (1 – )1/3
0.2
0.1
0
0 1× 105 2× 105 3× 105 4 ×105
time, sec
Figure 6.4. (a) Fraction of forsterite reacted away as a function of time. (b) The
reacted fraction data fit to the shrinking particle model. The slope of the line is kp.
Particle lifetime model
Geochemists are more familiar with the closely related particle lifetime
model of Lasaga (1998b), which is based on the same assumptions as the
shrinking particle model. This section shows how these models are related.
When the temperature, solution composition, and other rate-controlling
variables remain nearly constant, the dissolution flux (J = k+, mol/m2sec) for
a dissolving particle is constant and the release rate of the dissolving sub-
stance (r, mol/sec) is proportional to the particle’s surface area (A, m2).
r = Ak+ (6.55)
If the particle is a sphere with a diameter D (m) its area A = πD2, m2. The
rate of volume loss of the particle is the release rate multiplied by the molar
volume of the dissolving substance (Vm, m3/mol).
dV
= − rVm = −π D 2 k+Vm (6.56)
dt
The rate of volume loss of the particle is also given by the time derivative of
the volume change for a sphere (V = πD3/6).
dV π D 2 dD
= (6.57)
dt 2 dt
Setting these two equations equal to each other gives the rate of reduction
of the particle’s diameter.
dD
= −2 k+Vm (6.58)
dt
This equation is integrated and evaluated between the initial diameter (Do,
m) and the diameter (Dt, m) when time = t.
Dt = Do − 2 k+Vmt (6.59)
If the particle dissolves until its diameter reaches zero, it has dissolved away
completely and the particle lifetime (tl, sec) is found by setting Dt = 0 in Eq.
(6.59).
Do
t = (6.60)
2 k+Vm
Because Do/2 = Ro, the lifetime of the grain can also be expressed in terms
of kp as defined by Eq. (6.48). This relationship links the shrinking particle
and particle lifetime models.
Particle lifetime model 119
Ro 1
t = = (6.61)
k+Vm k p
Equation (6.60) can be rearranged to find the ratio of the original diameter
(Do) to the lifetime (tl).
Do
2k+Vm = (6.62)
t
Substituting Eq. (6.62) into Eq. (6.59) gives a relationship between dimen-
sionless diameter and dimensionless time.
Dt t
= 1− (6.63)
Do t
This equation shows that the particle diameter shrinks as a linear function
of elapsed time as shown in Figure 6.5.
The surface area of the particle is πD2, so the diameter can be
expressed in terms of the surface area and that relationship can be sub-
stituted into Eq. (6.63) to find the reduced surface area as a function of
reduced time.
1
A 2
D= (6.64)
π
2
At t
= 1 − (6.65)
Ao t
Equation (6.55) states that the rate of transfer of material (r, mol/sec) from
the particle to the solution is directly proportional to the surface area, so it
can be combined with Eq. (6.65) to give the dimensionless release rate as a
function of dimensionless time.
2
r t
= 1− (6.66)
ro t
1
6V 3
D= (6.67)
π
120 Surface kinetics
0.8
0.6
reduced dimension D/Do
0.4 A/Ao
3
Vt t
= 1 − (6.68)
Vo t
Multiplying the volume of the particle by the molar volume of the sub-
stance gives the number of moles of substance in the particle and the ratio,
nt/no, is the fraction of material remaining (p) at any time and p equals one
minus the fraction of material that has dissolved away (α).
3
nt t
= p = (1 − α ) = 1 − (6.69)
no t
1 − α = (1 − k pt )
3
(6.70)
Taking the cube root of both sides of this equation followed by rearrange-
ment produces the shrinking particle model derived earlier as Eq. (6.51).
Equation (6.63) can also be rearranged and expressed in terms of kp.
D t
1 − t = = k pt (6.71)
Do t
This relationship can be used for dissolving grains that are sampled and their
diameters measured over the course of an experiment. The slope of a graph
Particle lifetime model 121
1
r 2
t
1− t = = k pt (6.72)
ro t
(Robie and Hemmingway, 1995). The lifetime of a 1 mm diameter grain is cal-
culated using Eq. (6.60).
Do 1 × 10 −3 m
t = = = 8.33 × 109 sec (6.73)
2 k+Vm mol m3
2 2.75 × 10 −9 2 4.365 × 10 −5
m sec mol
Do 1 × 10 −3 m
t = = = 8.01 × 1013 sec (6.74)
2 k+Vm mol m3
2 2.75 × 10 −13
2.269 × 10
−5
m 2sec mol
The lifetime of a 1 mm diameter quartz grain is 2.54 million years, which means
that the quartz grain will persist 19,000 times longer than the forsterite grain.
122 Surface kinetics
Equilibrium condition
At equilibrium the chemical potential of each component in the solid equals
its chemical potential in the solution so ∆µr = 0 and the concentration of
each component in solution remains unchanged over time. However, dis-
solution and precipitation reactions continue to occur and their rates are
exactly equal because of the principle of detailed balance. This means that
although there is no net mass transfer between the solid and solution, these
reactions constantly reconstruct the surface. For most solids at room tem-
perature this effect is small but, at the high temperatures or for the long
times characteristic of many geological settings, surface reconstruction can
transfer significant quantities of trace element or isotopic species between
the surface and the solution. At low temperatures the rate of exchange of
components between the surface and the interior of the crystal is limited by
slow solid-state diffusion, so the composition of crystal interiors typically
remains unchanged over geologically significant times. The constant inter-
change of components between the solution and the surface allows the sur-
face defects and heterogeneities to anneal away over time so that the solid
surface evolves toward its equilibrium structure and composition.
The Wulff theorem states that the surface free energy of a crystal face is
directly proportional to the distance from the Wulff point, which lies in the
interior of the crystal (Mutaftschiev, 2001; Wulff, 1977). For amorphous
solids, the equilibrium shape is a sphere. For crystals, the equilibrium shape
can be predicted using the Wulff construction, which relates the surface free
energy of a crystal face to the length of a vector that is perpendicular to
that crystal face and passes through the Wulff point in the interior of the
crystal. Wolff and Gualtieri (1962) describe the methods used to predict
the equilibrium shape of a crystal. They also point out that crystals sel-
dom achieve a perfect equilibrium shape because they grow from a non-
equilibrium, supersaturated solution. Dissolution rate experiments typically
use broken grains produced by crushing. The surface free energy of some
sites on these grains can be quite high and these sites dissolve quickly when
Slightly undersaturated/slightly supersaturated conditions 123
J = J + (1 − e − ∆µr RT
) = J + (1 − e A RT ) (6.75)
124 Surface kinetics
ledge
F-face terrace
kink
S-face
K-face
The affinity model is based on the principle of detailed balance (see the
derivation in Chapter 3). According to this model, when A = 0, J = 0 and
when A > 5RT, J ≈ J+. This model seems to work for simple compounds, such
as quartz, that dissolve congruently and reversibly. It is easier to visualize
the goodness of fit by converting Eq. (6.75) to a linear form (Figure 6.7).
(J J + ) = 1 − e A RT (6.76)
(a) A/RT
0 2 3 4 5 6
4 10–7
3 10–7
J, mol/m2sec
2 10–7
1 10–7
0
0 5 10 15 20 30
A, kJ/mol
(b) A/R
0 1 2 3
1. 2
(b )
1
0. 8
The Kossel crystal model shown in Figure 6.6 provides the conceptual
basis for many models of crystal growth and dissolution. Most of those
models also incorporate the concepts of terrace, ledge, and kink (TLK),
which describe the free energy of attachment of a growth unit as being pro-
portional to the number of bonds (1 through 6) between it and the crystal.
According to the TLK model, the energy needed to form a 2D nucleus,
which initiates the growth of a new layer, is high, making the predicted rates
of crystal growth much lower than is observed. This problem was resolved
by the Burton–Cabrera–Frank (BCF) theory (Burton and Cabrera, 1949;
Burton et al., 1949, 1951; Cabrera and Burton, 1949) that postulates that
new ledges are the result of line dislocations that grow into spiral hil
locks (Figure 6.6). Zhang and Nancollas (1990) and Teng et al. (2000) give
examples of the application of these concepts to the growth of minerals.
Dissolution processes can be visualized as the inverse of the Kossel crystal
growth processes. For example, line defects at the spiral hillocks provide
steps that enhance dissolution rates leading to etch pit formation (Brantley
et al., 1986a; Brantley et al., 1986b). Etch pits and vacancy islands (Dove
et al., 2005) exert a profound control on the surface structure of dissolving
solids (Lasaga and Blum, 1986; MacInnis and Brantley, 1993). Etch pits
and the opening of vacancy islands cause weathered mineral surfaces to be
extremely rough (Berner et al., 1980; Velbel, 1989). Monte Carlo models of
dissolving Kossel crystal surfaces show that surface roughness increases with
the extent of the reaction and the departure from equilibrium (Bandstra and
Brantley, 2008; Wehrli, 1989). Because ledge retreat can initiate at crystal
edges, the edge regions dissolve more quickly than the centers of the faces,
so the dissolving crystal becomes rounded. Eventually the dissolving crystal
becomes bounded by rough and curved surfaces from which growth units
can be removed at any location (Cabrera and Vermiyea, 1958).
Incongruent dissolution changes the surface chemistry of dissolving min-
erals. Incongruent dissolution occurs when there are significant differences
among the bond strengths of the mineral’s components. This means that in
the initial stages of the dissolution process, weakly bonded components are
released quickly and strongly bonded components are retained in the sur-
face. The release rate of the weakly bonded components diminishes as their
concentration declines until their release rate matches the release rate of the
strongly bonded components. Eventually the dissolution becomes congruent
as the mineral surface becomes paved by a “leached layer” that is enriched
in the slowest dissolving components (Hellmann et al., 1990; Schott et al.,
2012). The leached layer can undergo reconstruction reactions that cross-
link the components, which further slows its dissolution rate (Casey et al.,
1993).
Very undersaturated/very supersaturated conditions 127
Diffusion and fluid advection are the most important mass transfer pro-
cesses in low temperature geochemical systems. Whenever mass transfer
processes are slow relative to chemical reaction rates, diffusion and advec-
tion processes will control the overall rate of a geochemical process. This
chapter introduces some simple models that link transport rates to chemical
reaction rates.
Advection
Advection rates are expressed in terms of volume discharge, specific dis-
charge, or surficial velocity depending upon the requirements of the model.
Discharge (also called volumetric flow rate) (Q = dV/dt, m3/sec) is the vol-
ume of fluid passing through a cross-sectional area per unit time. Specific
discharge (also called Darcy velocity) (q = Q/A, m/sec) is the volumetric flow
rate divided by the surface area through which the fluid is flowing. Surficial
velocity (v, m/sec) is the fluid velocity relative to an adjacent solid surface.
For a packed bed, the surficial velocity is the specific discharge divided by
the porosity (v = q/ϕ).
Viscosity is a fluid’s resistance to flow, so advection models must account
for viscosity either as an explicit term or in a scaling constant. The dynamic
viscosity (η, kg/m sec) equals the kinematic viscosity (ν, m2/sec) divided by
the density (ρ, kg/m3). The dynamic viscosity of liquid water in equilibrium
with the vapor phase for temperatures between the freezing and critical
points can be computed from an equation given by Watson et al. (1980).
The density of liquid water in equilibrium with the vapor phase for tem-
peratures between the freezing and critical points can be computed from an
equation given by Wagner and Pruß (2002).
128
Dimensionless numbers 129
k
v= ( A P )2 3 s1 2 (7.3)
n
In this model, k is a flow constant (= 1 m1/3/sec for SI units) that makes the
units work out and n (no units) is an empirically determined channel rough-
ness factor that ranges from ~0.03 for clean and straight natural streams
to ~0.05 for natural streams lined with cobbles and large boulders (Barnes,
1967; Limerinos, 1970). More sophisticated models have been developed
and calibrated (Bjerklie and Dingman, 2005). Darcy’s law predicts the spe-
cific discharge (q, m/sec) in porous media (packed beds) as a function of the
pressure drop across the bed (∆P, Pa) and the length of the bed (L, m).
k ∆P
q= (7.4)
η L
In this model, k (m2) is the permeability of the packed bed and η (kg/m sec)
is the dynamic viscosity of the fluid.
Dimensionless numbers
Dimensionless numbers are commonly used in diffusion and advection
models because they simplify the scaling of the models from laboratory
experiments to practical dimensions. This approach also takes advantage of
the Buckingham π theorem (Barenblatt, 2003), which states that n indepen-
dent variables with k independent dimensions can be expressed in terms of
p independent dimensionless numbers.
130 Diffusion and advection
p=n−k (7.5)
ρvL
Re = (7.6)
η
For this geometry, when Re < ~2300 the flow is laminar and when Re >
~4000 the flow is turbulent.
For flow through a packed bed with a porosity (φ, no units) and consisting
of grains with a diameter (D, m), the Reynolds number is
ρvD
Re = (7.7)
η (1 − φ )
When Re < ~10 the flow is laminar and when Re > ~2000 the flow is fully
turbulent.
For flow in a smooth-walled fracture (Qian et al., 2005) with an aperture
width of L, the Reynolds number is
ρvL
Re = (7.8)
2η
ρND 2
Re = (7.9)
η
Dimensionless numbers 131
η ν
Sc = = (7.10)
ρDi Di
The Schmidt number and Reynolds number can be combined to find the
Péclet number.
Pe = Re Sc (7.11)
The Péclet number is defined as the ratio of mass transfer by advection (v,
m/sec) to mass transfer by dispersion (DL, m2/sec) and diffusion (Di, m2/
sec).
vL
Pe = (7.12)
Di + DL
2 L2v 2
DL = (7.13)
105 Di
DL + Di
αL = (7.14)
v
This makes the Péclet number inversely related to the longitudinal dispersiv-
ity (Knapp, 1989).
1
Pe = (7.15)
αL
The 1 in the numerator of this fraction has units of meters. Empirical stud-
ies of flow in geological media have demonstrated that αL is related to the
scale of the study, L (m) (Table 7.1).
132 Diffusion and advection
Medium c m
α L = cLm (7.16)
( A / M )kL
Da I = (7.17)
meqv
If DaI >1, the reaction produces or consumes aqueous species faster than
they are delivered or removed by advection, so the overall rate is limited by
fluid flow. DaII compares the reaction rate to the rate of diffusive mass trans-
fer across the solid/solution interface.
Da II = Da I Pe =
( A M ) kL2 (7.18)
Di meq
If DaII >1, the reaction produces or consumes aqueous species faster than
they are delivered or removed by diffusion so the reaction rate is limited by
diffusion. A/M (m2/kg) is the ratio of the interfacial area between the solu-
tion and solid to the mass of solution; k (mol/m2sec) is the dissolution flux
rate constant; meq (mol/kg) is the equilibrium concentration; v (m/sec) is the
linear flow velocity; and Di (m2/sec) is the effective diffusion coefficient.
Fick’s laws
Diffusion is caused by the displacement of molecules due to random molecu-
lar motions. If the molecules are not homogeneously distributed in a region,
Fick’s laws 133
x+∆x
x ∆x
Jout
A
more of them will be displaced away from areas of high concentration into
areas of low concentration than vice versa. This causes net transport down a
concentration gradient. Erdey-Grúz (1974) and Robinson and Stokes (1959)
provide theoretical details about how and why diffusion occurs in aqueous
solutions. Cussler (2009), Denny (1993), Levich (1962), Probstein (1989),
and Vogel (1994) are useful resources for models that link diffusion and
advection. Berner (1980) and Lerman (1979) present models that are appli-
cable to sediments.
Fick’s laws are the basis for diffusive transport models. Fick’s first law
states that diffusive flux (J, mol/m2sec), is the product of a concentration
gradient (dc/dx, mol/m4) and a diffusion coefficient (D, m2/sec). Note that
concentration (c) is expressed with units of mol/m3.
dc
J = −D (7.19)
dx
Fick’s second law is related to the first by a simple derivation, which consid-
ers the volume element shown in Figure 7.1. As molecules diffuse into and
out of the volume element, the rate of change in the number of molecules
in the volume is the difference between the flux into and out of the volume
across the A areas at each end of the volume.
∆n
= − A( J x + ∆x − J x ) (7.20)
∆t
The change in the number of molecules divided by the volume of the ele-
ment is the concentration change in the volume.
∆n
∆c = so ∆n = ∆c ( A∆x ) (7.21)
A∆x
Substituting Eq. (7.21) into Eq. (7.20) gives the rate of change of concen-
tration in the volume.
134 Diffusion and advection
∆c ( A∆x )
= − A( J x + ∆x − J x ) (7.22)
∆t
∆c A ( J − Jx ) (J − Jx )
= − x + ∆x = − x + ∆x (7.23)
∆t A ∆x ∆x
The fluxes are defined in terms of the concentration gradients using Fick’s
first law.
∂c ∂c ∂c ∂c
D − D D −
∆c ∂x x + ∆x ∂x x ∂x x + ∆x ∂x x
=− =− (7.24)
∆t ∆x ∆x
dc d 2c
= −D 2 (7.25)
dt dt
Diffusion coefficients
Several experimental designs are used to measure diffusion coefficients
(Robinson and Stokes, 1959). Tracer diffusion coefficients for ions are
typically calculated from the limiting equivalent conductivity of the ions
(Lerman, 1979; Oelkers and Helgeson, 1988; Robinson and Stokes, 1959).
At 25°C, the diffusion coefficients for most ions range between 0.5 × 10−9 and
10 × 10−9 m2/sec (Li and Gregory, 1974; Miller, 1982; Oelkers and Helgeson,
1988). The diffusion coefficients for uncharged organic species are similar to
those for most ions (Oelkers, 1991).
The diffusion coefficients of H+ and OH− ions are much larger than
other ions because their movement in solution involves a kind of hop-
ping mechanism (Grotthuss mechanism) related to the very fast dissocia-
tion and re-association of water molecules (Cukierman, 2006; Tuckerman
et al., 2002).
Temperature, solution composition, and tortuosity 135
D = Ae − Ea RT
(7.26)
For aqueous species, graphs of log D versus 1/T are curved because Ea
decreases with increasing temperature (Figure 7.2). Values of Ea range from
around 20 kJ/mol at 25°C to near 10 kJ/mol at 300°C (Table 7.2). Ea for H+
and OH− are lower and their diffusion coefficients change somewhat less
with increasing temperature, but they are always larger than for other spe-
cies. Diffusion coefficients for aqueous species along the water liquid/vapor
curve are usually fit by a polynomial in 1/T.
b c
log D = a + + 2 (7.27)
T T
d log D
E a = 2.303R (
2c
= 2.303R b + T
d (1 T )
) (7.28)
136 Diffusion and advection
–7.5
H+
–8.0
Cl – OH –
logD
–8.5 Na +
Ca
–9.0 2+ SO
2
4 –
9.5
0.0017 0.0022 0.0027 0.0032 0.0037
1/T, K–1
–8.0 oc
tan oc
e tan
ol
–8.5 methane
logD
methanol
–9.0
acetic acid
be
nz
K BT (7.29)
D=
6πηR
Temperature, solution composition, and tortuosity 137
This equation can be recast to create a function that can be used to extrapo-
late diffusion coefficients to other temperatures.
ηT T
DT2 = 1 2 DT1 (7.30)
T1 ηT2
2.0 ×10–8
D, m2/sec
1.5 ×10–8
1.0 ×10–8
5.0 ×10–9
0
0 50 100 150 200 250 300
T, °C
–8.0
log D
–8.5
–9.0
–9.5
0.0017 0.0022 0.0027 0.0032 0.0037
1/T
The diffusion coefficient of H4SiO4 in pure water at 25°C is 1.17 × 10−9 m2/sec
(Rebreanu et al., 2008). This value can be extrapolated to hydrothermal
temperatures using Eq. (7.30). The easiest way to perform this extrapolation is
to log transform Eq. (7.30).
ηT T
log DT2 = log 1 + log 2 + log DT1 (7.31)
T1 ηT2
Values of η are calculated from Eq. (7.1) and used to find the first and second
terms of Eq. (7.31). The values of log (T2/ηT2) can be fitted to a temperature
function.
2
T 1 1
log 2 = −100039 − 289.64 + 7.6142 (7.32)
ηT2 T2 T2
At T1 = 298°C, log (η/T1) = −5.524 and log DT1 = −8.93. These values along
with Eq. (7.32) are inserted into Eq. (7.31) to produce a temperature function
for log DT2.
2
1 1
log DT2 = −100039 − 289.64 − 6.840 (7.33)
T2 T2
The values of log D and D predicted by Eq. (7.33) are shown in a graph in
Figure 7.3.
2DcationDanion
Dmutual = (7.34)
Dcation + Danion
In addition, ions can interact with each other to form ion pairs and
complexes. The diffusion coefficients of ion pairs can be estimated using a
method described by Applin and Lasaga (1984).
Diffusion rates in porous media are slower than in bulk solution because
the species must diffuse around the solid particles, which increases their dif-
fusion path length. This increase in diffusion path length is expressed as
Mineral dissolution without flow 139
tortuosity (τ), which is defined as the ratio of the actual distance traveled
by the species to the linear distance between the initial and final location of
the species. Tortuosity increases with decreasing porosity (φ) and there are
many theoretical and empirical models for this relationship (Lerman, 1979;
Shen and Chen, 2007). Three examples of empirical fits taken from Shen
and Chen (2007) illustrate the diversity of tortuosity models.
τ 2 = φ −1.14 (7.35)
τ 2 = φ + 3.79 (1 − φ ) (7.36)
τ 2 = 1 − 2.02 ln φ (7.37)
D
Dτ = (7.38)
τ2
x
c( x, t ) = cs erfc (7.39)
2 Dt
15 15 ×10 –3
100 sec
10 ×10 –3
concentration, molal
10 10 sec
concentration, mol/m3
1.0 sec
5 5 ×10 –3
0.1 sec
0 0
0 20 40 60 80 100
distance, µm
DCaSO4 is 0.95 × 10−9 m2/sec and gypsum solubility is 15 × 10−3 mol/kg (cs = 15
mol/m3) (Dutrizac, 2002). Using these values, Eq. (7.39) predicts the concentra-
tion versus distance profiles shown in Figure 7.4.
Transport-limited dissolution
Reactions are limited by the slowest step in the overall process. If the reac-
tion flux (JR, mol/m2sec) for a dissolving solid is fast compared to the dif-
fusion flux (JD, mol/m2sec), the dissolution process is limited by the rate of
transfer of reactants or products between the solution and the solid’s sur-
face. When JR ≈ JD, the overall rate of dissolution is controlled by mixed
kinetics, which is modeled by convolving the models for both the dissolution
and transport process.
The second Damköhler number, Eq. (7.18), is often used to test whether
a process is limited by diffusive transport. DaII is the ratio of the rate of
reaction at the surface (numerator) to the rate of transfer between the sur-
face and the solution (denominator). If DaII >> 1, the surface reaction is
fast relative to the diffusion rate. If DaII << 1, the surface reaction is slow
relative to the transport rate, so the overall rate is reaction limited. If DaII
Transport-limited dissolution 141
≈ 1, the surface reaction rate is nearly the same as the transport rate, so the
overall rate is controlled by mixed kinetics. The DaII test applies only for no
flow situations.
Da II = = 0.99
Di ceq (9.45 × 10 −9 )(15)
This means that for a 1 mm (0.001 m) wide fracture, the rate of gypsum dis-
solution nearly equals the rate of diffusive transport away from the dissolving
surface. For narrower fractures, DaII < 1, which means that the rate of diffusive
transfer is fast compared to the dissolution rate. DaII > 1 for wider fractures,
which means that the reaction rate is faster than the diffusive transfer rate.
Mass transfer coefficients are the basis for models where the dissolved
species are transported by a combination of diffusive and advective pro-
cesses. The diffusive mass transfer coefficient (kD, m/sec) is based on bound-
ary layer theory. The basic premise of boundary layer theory is that, for
laminar flow, the fluid velocity adjacent to a solid surface is zero (the “no
slip condition”) and the velocity increases as a parabolic function of dis-
tance away from the surface until it matches the velocity of the bulk fluid
(Figure 7.5). This means that there is a thin layer of fluid with a thickness
of δD (m) adjacent to the surface that is effectively static. The rate of mass
transport through this layer is limited by the diffusion rate of the dissolved
species. The diffusional boundary layer is much thinner than the velocity
boundary layer. For laminar flow past a flat surface, the thickness of the dif-
fusional boundary layer is related to the thickness of the velocity boundary
layer (δV) by the Schmidt number, which compares the fluid viscosity to the
diffusivity (Probstein, 1989).
142 Diffusion and advection
vf
flow velocity
0.6
δD = δV (7.41)
Sc1 3
Di
JD = (cs − c ) (7.42)
δD
D Lq 2 ν 3
1 1
Forced convection kD = 2.0 + 0.6 L = sphere diameter
around a sphere L ν D
J D = kD (cs − c ) (7.43)
The flux of material into or out of the boundary layer due to reactions at the
surface (JR, mol/m2sec) is driven by the difference between the surface con-
centration (cs, mol/m3) and the equilibrium concentration (ceq, mol/m3) and
is proportional to the dissolution rate constant (k+, mol/m2sec).
ceq − cs
J R = k+ (7.44)
ceq
144 Diffusion and advection
k+
JR =
ceq
(ceq − cs ) (7.45)
J R = kR (ceq − cs ) (7.46)
kRceq + kDc
cs = (7.48)
kR + kD
kR kD
J=
kR + kD
(ceq − c ) (7.49)
–3
JD
–4
JR
log J
–5
–6
10–8 10–6 10–4 10–2
q, m/sec
Figure 7.6. Dissolution flux of gypsum into pure water flowing past a flat plate.
The line labeled JR is the dissolution flux due to surface reaction only and
the line labeled JD is for the diffusion flux only.
1 1
D Lq 2 ν 3
kD = 0.646 i D
L ν i
1 1
9.45 × 10 −9 1 × 10 −3 q 2
8.94 × 10 −7 3
kD = 0.646 (7.50)
1 × 10 −3 8.94 × 10 −7
9.45 × 10 −9
1
kD = 9.52 × 10 −4 q 2
Inserting these values into Eq. (7.49) produces the results shown in Figure 7.6.
It is relatively easy to appreciate the mixed kinetics of this example. When
kR >>kD, kR + kD ≃ kD and J = kDceq and when kD >> kR, ks + kD ≃ ks and
J = kRceq.
146 Diffusion and advection
temperature, °C
100 50 25 0
–2
–3
log J, mol/m2 sec
JR
–4
JD
–5
The 25°C dissolution rate constant, 7 × 10−5 mol/m2 sec (Colombani, 2008), can
be extrapolated to other temperatures using the Arrhenius equation and an acti-
vation energy of 41.8 kJ/mol (Liu and Nancollas, 1971).
Ea 1 1 2183 2183
log kT2 = log kT1 + − = −4.15 + 7.32 − = 3.17 − (7.52)
2.303R T1 T2 T2 T2
Examples 7.4 and 7.5 show that diffusive transport tends to become a
limiting process for geochemical reactions as flow rates decrease and tem-
peratures increase. This means that although chemical reaction rates often
limit geochemical processes at or near the Earth’s surface, diffusional trans-
port becomes more rate limiting at depth because of higher temperatures
and slower flow rates due to lower rock permeability.
Coating growth model
The rate of interaction between minerals and solution species is sometimes
limited by the diffusion of a reactant from the solution through a coating
of insoluble reaction products (Figure 7.8). The diffusion rate through this
coating can be many orders of magnitude lower than the solution diffusion
rate, so that the rate of delivery of the reactant to the mineral surface is
much slower than its rate of reaction with the mineral. The thickness of the
coating (x, m) is a function of the number of moles of coating precipitated
per square meter (np, mol/m2), the molar volume of the coating (Vm, m3/mol),
the surface area of the grains (A, m2), and the porosity of the coating (φ).
The coating is likely to be porous if the molar volume of the coating is less
than the molar volume of the mineral or if some of the dissolving constitu-
ents are lost to the solution. This model is only applicable when x is much
smaller than the grain diameter.
148 Diffusion and advection
cr
concentration
coating
solution
solid
Figure 7.8. Conceptual
model showing how the
rate of reaction between
a mineral with a coating
of reaction products is
controlled by the rate of
diffusion of a reactant from
the surrounding solution 0
through the coating to the 0 x
surface. distance
n pVm
x= (7.53)
( − φ)
1
D
Jr = cr (7.54)
x
ν Dc (1 − φ ) 1
Jp = r r n (7.55)
Vm p
The terms in the parentheses can be combined into a rate constant (kp,
mol2/m4sec).
ν r Dcr (1 − φ )
kp = (7.56)
Vm
Coating growth model 149
dn p k p
= (7.57)
dt np
n 2p
= k pt (7.59)
2
np = ( )
2 k p t1 2 (7.60)
The thickness of the coating increases with the square root of time.
x = Vm n p = Vm ( )
2 k p t1 2 (7.61)
kp
dn p 1
dt
=
2
( )
2 k p t −1 2 =
2
−1 2
t (7.62)
This means that the thickness of the coating increases as a linear function
of t−½.
dx kp
= Vm −1 2
t (7.63)
dt 2
dnr 1 k p −1 2
= t (7.64)
dt ν r 2
This means that a graph of the rate of reactant consumption versus t−½ is
a straight line with a slope of (1/νr)(kp/2)½. The diffusion coefficient for the
reactant in the coating can be calculated from this slope by rearranging Eq.
150 Diffusion and advection
(7.56). This equation is only valid after enough time has passed to allow a
sufficiently thick coating to develop so that the flux of the reactant to the
mineral surface is smaller than its rate of consumption.
Huminicki and Rimstidt (2009) suggested that treating pyritic mine wastes with
bicarbonate solutions would produce an iron oxyhydroxide coating on the pyrite
surfaces, which would inhibit the development of acid mine drainage. The over-
all reaction converts iron from the pyrite into an iron oxyhydroxide coating.
FeS2(pyrite) + 3.75 O2 + 3.5 H2O + 4 HCO3−
= Fe(OH)3 (iron oxyhydroxide) + 4 H2CO3 + 2 SO42−
x = Vm ( )
2 k p t1 2 = (2.72 × 10 −5 ) ( )
2 ( 4.77 × 10 −14 ) t1 2 = 8.40 × 10 −12 t1 2
(a)
3.0×10–6
2.0×10–6
thickness, m
1.0 ×10–6
0
0 10 20 30 40 50
time, weeks
(b)
1.5×10–11
(b)
o×ygen consumption rate, mol/m2 sec
1.0×10–11
0.5 ×10–11
0
0 10 20 30 40 50
time, weeks
dn 4π R 2 dc
=− D (7.65)
dt νr dR
If the diffusion rate of reactant through the coating is fast compared to the
motion of the core/coating interface, a steady-state concentration gradient
will develop so that Eq. (7.64) can be integrated across R to find the rate of
reaction in terms of R and Ro.
dn 4π Dc RRo
= (7.66)
dt ν r Ro − R
coating
unreacted
core
Ro
Figure 7.10. Schematic
diagram showing the
conceptual basis for the
shrinking core model.
Shrinking core models 153
R3
α = 1− (7.67)
Ro3
d α 3R 2 dR
= 3 (7.68)
dt Ro dt
Substituting Eqs (7.67) and (7.68) into Eq. (7.66) gives the rate of trans-
formation in terms of the fraction of the particle reacted.
d α 3VmDc (1 − α 1 3 )
= (7.69)
dt ν r Ro2 1 − (1 − α 1 3 )
2 2V Dc 2VmDc
1 − α − (1 − α ) = m 2 t = k pt
23
3 νr Ro k p = ν R 2 (7.70)
r o
0.8
0.6
0.4
0.2
0
0 20 40 60 80 100
time, yr
Vm ( goe ) 2.08 × 10 −5
φ = 1− = 1− = 0.13
Vm ( py ) 2.39 × 10 −5
DO2 2.3 × 10 −9
Dτ O2 = = = 2.49 × 10 −10 m 2 /sec
τ2 9.24
cr = 0.27 mol/m 3
Ro = 0.005 m = 0.5 cm
Shrinking core models 155
2 2
1 − α − (1 − α ) 1 − α − (1 − α )
23 23
t= 3 = 3
kp 1.05 × 10 −10
This model shows that the rate of conversion of pyrite to goethite is quite fast
during the initial stages of the process so that 50% of the pyrite is converted dur-
ing the first 11 years of exposure to oxidizing conditions. As the layer of goethite
grows thicker, the diffusion distance becomes long and the rate of conversion
slows so that it takes an additional 89 years to convert the remaining pyrite to
goethite.
Chapter 8
Quasi-kinetics
Quasi-kinetic models deal with processes that are controlled by mass trans-
fer rates rather than by chemical reaction rates. These models assume nearly
instantaneous attainment of equilibrium within the region of interest, so
changes in the species distribution are controlled by the rate of transfer of
substances into or out of that region. These models are constrained by con-
tinuity equations making them similar to the chemical reactors models in
Chapter 4.
Da I =
RL
=
{(A M ) J } L (8.1)
meqv meqv
R (molal/sec) is the reaction rate and meq (molal) is the equilibrium con-
centration of the reactant. v (m/sec) is the surficial velocity of the fluid and L
is the characteristic length of the domain. A rule of thumb is that when DaI
< 0.1, less than 10% of the aqueous species will react as the solution transits
L; and when DaI > 10, more than 90% of the aqueous species will react over
156
Local equilibrium assumption 157
this distance. Even if the reaction rate is fast, the rate of reaction between
solution and the mineral might still be limited by the rate of diffusive and
dispersive transport to or from the minerals’ surfaces. The Péclet number,
Pe, compares the rate of species transport to or from the minerals’ surfaces
by diffusion or dispersion to the rate that they are added to or removed from
the domain by advection. The Péclet number is the ratio of mass transfer
by advection to the rate of mass transfer to or from the minerals’ surfaces
expressed by the sum of the diffusion coefficient (Di, m2/sec) and the longi-
tudinal dispersion coefficient (DL, m2/sec).
vL
Pe = (8.2)
Di + DL
When Pe < 50, diffusion and dispersion dominate and aqueous species are
delivered to or removed from the surface faster than they are added to or
removed from the domain by advection. When Pe > 100, the rate of advec-
tion is so fast that the reacting species are transported into or out of the
domain before they can migrate to or from minerals’ surfaces. Figure 8.1
shows a map of Da and Pe. The local equilibrium assumption is valid in the
upper right corner of the map where Da > 10 and Pe > 100 but reactions
will be unlikely to reach equilibrium for the other conditions shown because
they are either transport or reaction limited, or both.
The quartz geothermometer (Fournier, 1977) requires that the dissolved silica
concentration in an ascending hydrothermal solution becomes “quenched in” as
the fluid approaches the surface. On the other hand, the phase transfer model for
the rate of quartz deposition, described in this chapter, requires that the LEA be
valid. The conditions for which each of these models is appropriate can be found
by calculating values of DaI and Pe. Figure 8.1 shows the temperature range for
the transition from LEA to quenched conditions.
The rate of quartz precipitation, Rp (molal/sec), (Rimstidt and Barnes, 1980)
is needed to calculate DaI.
Rp =
A
M (
k+ 1 − Q K ) (8.3)
For a fracture, the surface area to mass of solution ratio (A/M = 2/1000ρW)
depends on the fracture width (W, m) (see Table 6.1). This model considers a
500 µm wide fracture so (A/M) = 4. To simplify the calculations, the solution
density (ρ) is set equal to one over the entire temperature range and the degree
of supersaturation (Q/K) is set equal to two, which makes Rp = −(A/M)k+. The
158 Quasi-kinetics
8
10
6 LEA Valid
10
Fast reaction Fast reaction
4
Slow transport Fast transport
10
2
10 300°C
200°
1
Da
100°
–2
10 25°
–4 Slow reaction
10 Slow transport
Slow reaction
–6 Fast transport
10
–8
10
5
0.01 0.1 11 0 100 1000 10 4 10
Pe
Figure 8.1. Map of DaI versus Pe that shows that the Local Equilibrium
Assumption (LEA) is valid for DaI >~10 and Pe > ~100. The dots are values of the
DaI and Pe values for quartz precipitation in a fracture at different temperatures
calculated in Example 8.1. They show that the LEA is valid for T > ~250°C.
rate constant, k+, for each temperature is computed from Rimstidt and Barnes
(1980) and Rimstidt (1997b).
Da I =
RL
=
(2.77 × 10 −13 )(10) = 1.51 × 10 −3 (8.7)
meqv (1.83 × 10 −4 ) (1 × 10 −5 )
Note that if L = 100 m, DaI = 1.5 × 10−2, then DaI scales directly with L. This
value of DaI falls in the “slow reaction” part of Figure 8.1, indicating that the
solution does not reach equilibrium within a 10 or even 100 m fracture length
because of the very slow precipitation rate of quartz.
The diffusion coefficient for dissolved silica, Di (m2/sec), and the longitudinal
dispersion coefficient, DL (m2/sec), are needed to compute the Péclet number.
Di is calculated from equation (5) in Rebreanu et al. (2008) and DL is calculated
from the fracture width, using the Horne and Rodriguez (1983) model.
DL =
2 2 ( )
W 2 v2
(8.8)
105 Di
vL
Pe = (8.9)
DL + Di
DL =
2 ( )
W 2 v2
2 =
2 ( )
0.0005 2 1 × 10 −5 2
2 ( )
= 4.3 × 10 −9 (8.10)
105 Di 105 1.1 × 10 −9
vL
Pe = = 2.3 × 10 4 (8.11)
DL + Di
This result means that the transport of silica to the quartz surface by dispersion
and diffusion is very fast and is not rate limiting.
Values of DaI and Pe calculated for other temperatures are plotted on
Figure 8.1. The figure shows that the LEA holds for high temperatures (>~250°C)
but fails for T <~200°C. This is consistent with using the quartz geothermom-
eter to find the temperature of geothermal reservoirs. The quartz geothermom-
eter is based on the assumption that the dissolved silica concentration remains
unchanged as the solution flows to the surface (Fournier, 1977).
160 Quasi-kinetics
Box models
Some models assume that a system reaches a steady state rather than equi-
librium. Equilibrium is defined by the principle of detailed balance, which
requires that the forward and reverse rates are equal and that each step along
the reaction path is reversible. The forward and reverse rates of steady-
state processes are equal but the process steps that produce the forward
rate are different from those that produce the reverse rate. At steady state,
the state variables of an open system remain constant even though there is
mass and/or energy flow through the system. The steady-state assumption
is especially useful for processes that occur in a series, because the concen-
trations of intermediates that are formed and subsequently destroyed are
constant. Perturbation of a steady-state system produces a transient state
where the state variables evolve over time and approach a new steady state
asymptotically.
Box models are steady-state models that are used to follow the flow of
a conserved substance between reservoirs in the geological environment
(Lasaga, 1980b; Lasaga and Berner, 1998; Lerman and Wu, 2008). Box
models divide the environment into various reservoirs, represented by boxes
(Figure 8.2), each of which represents a major part of the overall system. If
the model considers cycling of a substance throughout the entire Earth, the
reservoirs might be the oceans, the atmosphere, the crust, etc. For example,
Figure 8.3 shows the main reservoirs and fluxes for the global water cycle
(Berner and Berner, 1987). It is presumed that each of these reservoirs is
well mixed, so the concentration of the substance is approximately the same
throughout. Material is transferred from one reservoir to another via fluxes
(F) of the substance, which are represented by arrows. Simply constructing a
diagram shows how a substance flows through the system (Figures 8.2a and
8.3). In addition, box models are used to quantify how the overall system
might respond to a perturbation.
The simplest box model consists of one reservoir that is fed by a con-
stant flux (F0) (Figure 8.2). The flux out of a reservoir is the product of the
mass of the substance (M) in the reservoir and a mass transfer constant
(k). Box model reservoirs are simply gigantic ideal mixed flow reactors (see
Chapter 4) where there is no net generation or consumption of the sub-
stance. Over geologic time spans these reactors tend to attain a steady state
so that the flux into a reservoir is matched by the flux out, which means that
the rate of accumulation in the reservoir is zero and M is constant (Mss). The
flux out of the reservoir equals a mass transfer constant (k, yr−1) times the
mass of the substance in the reservoir.
F0 = F = kM ss (8.12)
Box models 161
(a)
reservoir
flux in flux out
F0 F1 = k 1 M 1
M1
(b) 12
Figure 8.2. (a) Schematic
illustration of a box model
Mss with a constant flow (Fo) of
10
material into the reservoir,
which contains M amount
8 of material. The flux out
of the reservoir (F1) is
a function of the mass
6 transfer constant (k1) and
M
F0 F
k= = (8.13)
M ss M ss
Once the mass transfer constant is known, it is possible to model the effects
of a perturbation on the reservoir contents and output. The rate of change
of mass of material in the reservoir is the difference between the flux into
and the flux out of the reservoir.
dM
= F0 − kM (8.14)
dt
dM
dt = (8.15)
F0 − kM
162 Quasi-kinetics
atmosphere
13 10 9 kg
vapor transport
110 10 9 kg/yr
386 10 9 kg/yr
423 10 9 kg/yr
precipitation
73 10 9 kg/yr
precipitation
evaporation
37 10 9 kg/yr
evaporation
ice
29000 10 9 kg
t M
dM
∫ dt = ∫ F0 − kM
(8.16)
0 M0
1 F − kM
t = − ln 0 (8.17)
k F0 − kM 0
F0 − kM
e − kt = (8.18)
F0 − kM 0
F0 F0
M= − − M 0 e − kt (8.19)
k k
Figure 8.2b, which is based on Eq. (8.19), shows how an empty reservoir
returns to the steady state. During the initial stage of the process the res-
ervoir contains very little of the substance, so the rate of removal by kM
is low and M grows quickly. However, increasing M over time causes kM
to increase and that causes the rate of growth of M to slow. Eventually, M
becomes large enough to cause kM to equal F0 and M has returned to Mss.
Because M approaches Mss as an asymptote, it is not possible to mathem-
atically define a time for the return to steady-state conditions. However, it
is possible to find a practical time by defining a time constant (tc) for the
system.
Box models 163
1
tc = (8.20)
k
If 1tc is substituted into Eq. (8.19), the calculated value of M is 63% of Mss.
At 3tc, M is 95% of Mss and at 5tc, M is 99% of Mss. Somewhere between 3tc
and 5tc the difference between M and Mss becomes smaller than the uncer-
tainty in determining M so that there is no longer a way to distinguish M
from Mss and the system has effectively reached steady state.
Although the box models treat substances as a continuum, they actually
consist of discrete particles (atoms, molecules, clusters, etc.) and it is useful
to consider the average time that a single particle might remain in a reser-
voir. Some particles will persist in the reservoir for longer times and some
for shorter but the average time spent in the reservoir is defined by the resi-
dence time (tres).
M ss
tres = (8.21)
F0
For a one-reservoir model, tres = tc; see Eq. (8.13). Residence-time analysis is
a useful way to understand the behavior of tracers, such as tritiated water, in
geochemical cycles. Residence-time analysis methods for more complicated
situations are described in many chemical engineering textbooks (Hill Jr.,
1977).
Although one-box models are informative, most geochemical cycles con-
sist of multiple reservoirs coupled by many fluxes; see Figure 8.3 for exam-
ple. Models for multiple reservoir–multiple flux situations and for coupled
geochemical cycles are somewhat more complicated but are well under-
stood (Lasaga, 1980a, 1980b; Lasaga and Berner, 1998; Lerman and Wu,
2008). One important insight that comes from these more complex models
is that increasing complexity seems to increase the stability of the system.
This means that environmental perturbations caused by natural processes or
human activities tend to be damped by complex natural systems.
Figure 8.3 shows that for the global water cycle, the sum of the water evaporated
from the oceans and the continents is 496 × 109 kg/yr. At steady state the pre-
cipitation flux must have the same value. The water content of the atmosphere
can be modeled as a single reservoir with a water content of 13 × 109 kg, so k =
(496 × 109 kg/yr/13 × 109 kg) = 38.2 yr−1. For this simple one-box model, tc = tres
= 0.026 yr = 9.6 days. This means that the average residence time of a water mol-
ecule in the atmosphere is slightly less than 10 days, so we expect that the rainfall
flux should respond quickly to a change in the evaporation flux. For example, if
164 Quasi-kinetics
time, days
0 20 40 60 80
9
13.0 10
9
12.0 10
9
11.5 10
0 0.05 0.1 0.15 0.2 0.25
time, yr
the temperature and pressure can be found by taking the derivative of the
mineral’s temperature and pressure solubility function.
∂m ∂m
dm = dT + dP (8.22)
∂T P ∂P T
∂m ∂m
αT = and α P = (8.23)
∂T P ∂P T
dm
= αT ZT′ + α P ZP′ (8.24)
dZ
If the solution has a Darcy velocity (q = dZ/dt) along the direction of the
pressure and temperature gradients, the rate of removal of the dissolved
constituent from the solution (molal/sec) is
dm
R= = q (αT ZT′ + α P ZP′ ) (8.25)
dt
Dividing R by (A/M) gives the flux (J, mol/m2sec) of the dissolved compo-
nent onto a depositional surface and multiplying J by the molar volume
(Vm, m3/mol) of the depositing mineral gives the rate of increase in thickness
(dz/dt, m/sec) of a growing layer of that mineral (see Chapter 6).
dz M
= qVm (αT ZT′ + α P ZP′ ) (8.26)
dt A
Most quartz veins are the result of the cooling and decompression of flowing
hydrothermal solutions. Because quartz solubility increases with increasing tem-
perature and pressure, the presence of quartz veins in an outcrop indicates that
silica-bearing solutions were flowing along a fracture system away from a hot,
high-pressure source toward a cooler, lower pressure discharge point.
166 Quasi-kinetics
0.3
0.25
0. 2
dz/dt, cm/yr
0.15 LEA valid
0. 1
0.05
0
0 50 100 150 200 250 300
T, °C
Figure 8.5. Growth rate of a layer of quartz on the wall of a fracture with a 1 cm
aperture predicted by the phase transfer model. The quartz is precipitating from
an aqueous solution flowing with a Darcy velocity of 1 m/day along a geothermal
gradient of 25°/km.
This example estimates the rate of growth of quartz on the wall of a 10 cm
wide (W = 0.1 m) fracture as the result of the upward flow of a hydrothermal
solution at the rate of 4.3 m/day along a geothermal gradient of 25°/km.
A 2
= (Table 6.1) ρ for liquid water is from Wagner and Pruß (2002)
M ρW
For relatively shallow earth conditions, we can assume that pressure has little
or no effect on solubility so that αP = 0.
When these parameters are inserted into Eq. (8.26) they predict that the rate of
growth of quartz on the walls of the fracture (Figure 8.5) ranges from 0.25 cm/yr
Reaction path models 167
ni = nio + ν i ξ (8.27)
Reaction path models are local equilibrium models that assume that all the
reaction rates are so fast that the system is completely equilibrated at each
step. This means that the rates of formation or destruction of all the species
in the model are tied by the extent of reaction to the rate of addition, or
reaction, of the titrant.
dni dξ
= νi (8.28)
dt dt
If the rate of addition of the titrant is specified, then the rates of formation
or destruction of all the other phases can be calculated from Eq. (8.28); for
example, see Zhu and Lu (2009).
Reaction path models are unreliable whenever the local equilibrium
assumption fails. This occurs whenever slow reaction kinetics does not allow
a phase to nucleate and grow or when a phase dissolves away too slowly.
This means that predictions of reaction path models must be evaluated by
a knowledgeable geochemist to identify these refractory phases. The typ-
ical remedy for this problem is to remove the uncooperative phase from the
model and allow a more reactive, but metastable, phase to take its place.
168 Quasi-kinetics
When rainwater interacts with a rock that contains potassium feldspar, the feld-
spar dissolves and new minerals grow until the aqueous solution comes to equi-
librium with the potassium feldspar. The reactions add K+ and H4SiO4 to the
solution and consume H+. This reaction path model tracks the changing solu-
tion composition on an activity diagram (Figure 8.6) to show how the solution
composition traverses the stability fields of gibbsite (gib), kaolinite (kaol), and
muscovite (mu) until it reaches equilibrium with potassium feldspar (Kf). The
phase boundaries are defined by five reactions.
1. kaol + 5 H2O = 2 gib +2 H4SiO4 log K1 = −8.50
2. mu + 9 H2O + H+ = gib + K+ + 3 H4SiO4 log K2 = −9.38
3. 2 mu + 3 H2O + 2H = 3 kaol + 2 K
+ +
log K3 = 6.74
4. 3 Kf + 12 H2O + 2H+ = mu + 2 K+ + 6 H4SiO4 log K4 = −14.41
5. 2 Kf + 8 H2O + 2H+ = kaol + 2 K+ + 4 H4SiO4 log K5 = −7.36
This model tracks the changes in the concentrations of H+, K+, and H4SiO4 as
small amounts of potassium feldspar are added to 1 kg of rainwater solution.
The feldspar first alters to gibbsite, then to kaolinite and then to muscovite. The
reaction for each path segment is:
mu Kf
4
D
3
2
log aK /aH
1 gib kaol
0 B
A
–1
–2
–5.5 –5 –4.5 –4 –3.5 –3
log aH4SiO4
The stoichiometric coefficients, νi, are moles of each species produced or con-
sumed per mole of potassium feldspar reacted.
Slightly acidic rainwater has a H+ concentration of about 3.2 × 10−5 mol/kg
(pH ~4.5) and a K+ concentration of about 3.2 × 10−6 mol/kg (0.12 ppm) and
contains about 5.8 × 10−6 mol/kg (0.1 ppm) dissolved silica.
The reaction path begins at log aK/aH = −1.0 and log aH4SiO4 = −5.2. As the
amount of K and Si in the system is incremented by the addition of potassium
feldspar, reaction A occurs followed by reactions B, C, and D. The stepwise
increase in the extent of reaction, ξ, is reflected by changes in the concentration
(activity) of the aqueous species and by changes in the number of moles of each
mineral.
mH+ = ( mH+ )o + ν H+ ξ
mK + = ( mK + )o + ν K + ξ
The size of the ξ step must be slightly smaller than the change in hydrogen ion
concentration at the end of the reaction path otherwise its concentration will
become negative. This choice can be made by trial and error.
The resulting diagram shows that the first increments of potassium feldspar
are converted to gibbsite as the solution tracks along path A. Then the gibbsite
formed along path A dissolves and its components along with those from the
added potassium feldspar are converted to kaolinite along path B. After the
gibbsite is consumed, potassium feldspar is converted to kaolinite along path
C until the solution composition reaches saturation with respect to muscovite.
Finally, potassium feldspar is converted to kaolinite and muscovite until the
solution composition reaches the muscovite–kaolinite–potassium feldspar triple
point. At that point the solution is in equilibrium with potassium feldspar and
no further reaction takes place.
Garrels and MacKenzie (1967) published the classic example of the application of
a mass balance model to account for chemical weathering in the granitic rocks of
the Sierra Nevada batholith. The model, which uses water analyses from Feth et al.
(1964), takes the chemical analysis for snowmelt as an initial water composition
and the average ephemeral spring water as a final water composition. It is presumed
that the infiltrating snowmelt dissolved CO2 from soil gases to make H2CO3 that
reacted with plagioclase, potassium feldspar, biotite, and quartz to form kaolinite.
The difference between the snowmelt and spring water compositions (∆m) is the
result of the reaction of the water with the minerals in the batholith.
The minerals in the batholith that are available for weathering reactions are pla-
gioclase, which consists of albite (NaAlSi3O8) and anorthite (CaAl2Si2O8); bio-
tite (KMg3AlSi3O10(OH)2); potassium feldspar (KAlSi3O8); and quartz (SiO2).
The reactions produce kaolinite (Al2Si2O5(OH)4).
There are five postulated weathering reactions. They are written as a series of
linear equations. This is the key feature of the mass balance approach.
There are six equations and five unknowns. The extent of each reaction can be
estimated using multiple linear regression. Multiplying the estimated extent of
reaction, ∆ξi (values in parentheses are the standard error of the estimate) by the
reaction coefficient for the mineral gives the number of moles of that mineral
consumed per kilogram of water. If the contact time, ∆t, is known, the rate of
each reaction can be estimated from ∆ξi /∆t.
Garrels and MacKenzie (1967) solved for these values using a sequential sub-
traction process. The more general multiple linear regression method used here
distributes the rounding and analytical errors over the ξi estimates and is amen-
able to expansion to a much larger set of reactions.
Partition coefficients
The trace element concentrations in geologic materials are strongly linked to
major element behaviors. This means that trace element distributions offer
useful clues about processes that control the distribution and redistribution
of major elements. The distribution of a trace element between a solid and
Mixed flow reactors 173
X mTr
λ = Tr m (8.29)
XM M
This definition incorporates the activity coefficients for the trace element
in the solid and solution into λ; see McIntire (1963) for an approach that
explicitly accounts for activity coefficients.
The model assumes that during the reaction process the solid’s surface
is always in equilibrium with the solution. Then it becomes completely iso-
lated when it is overgrown by the next increment of precipitating solid. At
low temperatures, the diffusion rate of trace elements in most solids is very
slow, so the trace element concentration in each layer remains unchanged
once that layer is overgrown.
If λ > 1, the trace element is enriched in the solid relative to the solution
during precipitation, so mTr/mM is reduced as precipitation proceeds. If λ <
1, the trace element is enriched in the solution relative to the solid during
precipitation, so mTr/mM increases as precipitation progresses. In a mixed
flow reactor operating at a steady state, mTr/mM remains constant and the
solid has a constant XTr/XM ratio throughout. If the solid precipitates in a
batch reactor, mTr/mM changes during the precipitation process causing the
XTr/XM ratio to change from the center to the surface of the precipitating
solid.
X Tr RTr
= =
(( m ) − ( mTr )out )Q ( mTr )in − ( mTr )out
Tr in
= (8.30)
X M RM (( m
M )in − ( mM )out )Q ( mM )in − ( mM )out
Combining Eqs (8.29) and (8.30) and realizing that the effluent solution is a
sample of the solution in the reactor leads to
174 Quasi-kinetics
λ=
( mTr )in − ( mTr )out ( mTr )out
=
(( m ) ( mTr )out ) − 1
Tr in
(8.31)
( mM )in − ( mM )out ( mM )out (( m
M )in ( mM )out ) − 1
This means that λ can be computed from the difference between the con-
centration of Tr and M in the feed and effluent solutions regardless of the
flow or precipitation rate. These same data along with the flow rate can be
used to find the precipitation rate of the major element, so that a mixed flow
reactor experiment is well suited to quantify the effect of the precipitation
rate on λ.
λ=
(1 pTr ) − 1 (8.32)
(1 pM ) − 1
Equation (8.32) is then written in terms of αTr and αM.
λ=
(1 (1 − α )) − 1
Tr
(8.33)
(1 (1 − α )) − 1
M
1 1
= λ − 1 + 1 (8.34)
1 − αTr 1 − αM
Taking the reciprocal of both sides of Eq. (8.34) allows it to be solved for αTr.
Batch reactors 175
0.8
0.6 Co ( =3.7 )
Tr
=1
0.4
0.2
Sr ( =0.04)
0
0 0.2 0.4 0.6 0.8 1
M
1
αTr = 1 − 1 λ − 1 + 1 (8.35)
1 − αM
This equation and λ values from Curti (1999) can be used to model the incorpor-
ation of Co and Sr into calcite precipitating in a mixed flow reactor.
λCo = 3.7
λSr = 0.04
Graphs of Eq. (8.35) shown in Figure 8.7 demonstrate that this scheme could
work for Co because precipitating 20% of the Ca would remove nearly 50% of
the Co. Using a series of mixed flow reactors would further improve the process
because removal of 50% of the remaining Co in a second reactor would leave
25% and a third reactor would further reduce the Co to 12.5%. On the other
hand, this process would be inefficient for Sr because precipitating 80% of the
Ca would only remove about 10% of the Sr. Hence, co-precipitation in a mixed
flow reactor is only effective for trapping trace elements if λ > 1.
Batch reactors
Co-precipitation in batch reactors is modeled using the Doerner–Hoskins
approach (Doerner and Hoskins, 1925). Because batch reactors are closed
176 Quasi-kinetics
systems, differences in the removal rate of the trace and major elements from
the solution cause the solution to become relatively enriched or depleted
in the trace element as the precipitation proceeds. This simple model from
Curti (1997) considers how removing a small amount of Tr and M by pre-
cipitating a small mass of solid affects the composition of the solution.
X Tr − ∆mTr m + ∆mTr
X = = λ Tr (8.36)
M surface − ∆mM mM + ∆mM
As the amount of solid precipitated tends to zero, ∆mTr → dmTr and ∆mM
→ dmM so
dmTr m
= λ Tr (8.37)
dmM mM
m m
ln Tr
o
= λ ln Mo
(8.38)
mTr mM
This expression allows the calculation of λ from the initial and final con-
centrations of the trace and major element.
m m
λ = ln Tr
o
ln Mo
(8.39)
mTr mM
Applying the antilog transform to Eq. (8.38) gives a relationship that pre-
dicts the concentration of the trace element in the solution as a function of
the amount of major element precipitated.
λ
mTr mM
= o (8.40)
mTr
o
mM
The concentrations of Fe, Mn, and Mg in calcite cements are often used to inter-
pret their depositional environment. These interpretations use distribution coef-
ficients to estimate the composition of the solution based on an analysis of the
calcite cement. Laboratory determinations of these distribution coefficients are
typically carried out in batch reactors. Creating a model of the precipitation pro-
cess in these experiments is helpful for understanding the uncertainties associated
with the measurements. If the calcite is precipitated in a batch reactor, Eq. (8.40)
Batch reactors 177
Mn ( =14.8)
0.8
Fe ( =2.7)
0.6
Tr
=1.0
0.4
0.2
Mg ( =0.097)
0
0 0.2 0.4 0.6 0.8 1
M
pTr = ( pM )
λ
(8.41)
This equation can be further recast in terms of the fraction of the element
precipitated.
1 − αTr = (1 − α M )
λ
(8.42)
Equation (8.42) can be rearranged to find the fraction of Tr that has precipitated
as a function of the fraction of M that has precipitated using the λ values listed
below.
αTr = 1 − (1 − α M )
λ
(8.43)
Figure 8.8 shows that within analytical resolution all the Mn is removed from
solution by the time 30% of the Ca has precipitated. This means that to deter-
mine λMn, only small amounts of Ca can be allowed to precipitate and the value
of λMn will be strongly affected by the Ca analysis. On the other hand, deter-
mining λMg requires the precipitation of most of the Ca and the accuracy of
λMg depends upon accurate Mg determinations. Values of λ near 1 are easiest to
determine and most accurate. Note that these same constraints apply to deter-
mining distribution coefficients using mixed flow reactors (see Figure 8.7).
For near surface conditions, mineral reaction rates slow and the flow rates
are relatively faster so the local equilibrium assumption often fails and chro-
matography models are less useful for modeling mineral transformations.
However, many adsorption–desorption processes are fast enough to main-
tain local equilibrium so that chromatography models are useful for mod-
eling the retardation of species that adsorb to aquifer minerals. The central
feature of these models is the retardation factor Rf (no units), which is the
ratio of the velocity of transport of the adsorbing species to the velocity of
a non-adsorbing, unreactive tracer (Higgins, 1959).
ρb
Rf = 1 + Kd (8.44)
φ
ρb (g/cm3) is the bulk density of the packed bed, φ (no units) is the por-
osity, and Kd (cm3/g) is the distribution coefficient. See Zheng and Bennett
(1995) for a derivation of this equation. It is related to the KD (mol/g-solid/
mol/kg-water) derived in Chapter 6 by a units transformation involving the
molecular weight of the adsorbing species and the density of water. Because
Kd is related to KD, its value is also a function of the specific surface area
of the solid.
mol i
g solid 1 = 1000 K = K cm soln
3
KD × (8.45)
d
mol i g solid
D
kg soln
kg soln 1000 cm 3soln
Rf − 1 26 − 1
Kd = = = 17.9 (8.46)
ρb φ 0.92 0.66
(a)
c/co
0.5
0
0 1 2 3 4 5
pore volumes
(b)
1
c/co
0.5
0
0 20 40 60 80 100
pore volumes
Figure 8.9. (a) Breakthrough curve for three Br− tracer experiments showing that
the inert tracer breakthrough occurs at one pore volume. The dashed line shows
the pattern expected for ideal plug flow and the deviation of concentration from
this line is the result of dispersion and diffusion. (b) Breakthrough curve for the
Zn2+ experiment showing that c/co occurs after 26 pore volumes have eluted.
Chromatography models and retardation factors 181
Polymerization rates
Polymers such as DNA, proteins, polysaccharides, and polyphenols are
essential components of organisms; and synthetic polymers such as Nylon,
polyethylene, styrene, and Teflon are basic raw materials for modern tech-
nology, so it is not surprising that there is a rich literature about the for-
mation of polymers. Most of that literature focuses on the formation of
biopolymers and plastic but the polymerization models (Dotson et al., 1996)
for those cases are potentially useful to geochemists.
Many biopolymers and synthetic polymers are linear, which means
that their constituent monomers are linked up in a head to tail fashion.
For example, cellulose, (C6H10O5)n, is a linear array of glucose molecules,
and polyethylene, (C2H2)nH2, is a linear grouping of ethylene molecules.
However, polymers linked in three dimensions play a more important role
in geochemical situations. For example, humic acids are randomly linked
three-dimensional polymers consisting of various degradation products of
biochemical substances.
Models for the formation of linear polymers are relatively simple (Dotson
et al., 1996). The chemical reaction for the growth of an (AB)n linear poly-
mer involves hooking up additional (AB)m molecules head to tail.
182
Polymerization rates 183
The overall rate of this reaction depends upon the rate of interaction between
the A and B functional groups with concentrations A and B.
dA dB
= = − kAB (9.2)
dt dt
dP
= − kP 2 (9.3)
dt
Po
P= (9.4)
1 + ktPo
dP1
= −2 kPP
1 (9.5)
dt
The factor of two in this equation arises because there are two distinguishable
reactions, one for each end of the polymer. Additional equations can be writ-
ten to find the rate of change in the concentration of polymers of each chain
length (Dotson et al., 1996). This model has didactic value, but most geo-
chemical polymers are three-dimensional so they have more than two attach-
ment points. Models of three-dimensional polymerization networks, such as
occur with silica polymerization (Jin et al., 2011), can become quite complex.
K
H 4SiO 4 + H 4SiO 4
1
H6Si2 O7 + H2 O (9.6)
K
H6Si2 O7 + H 4SiO 4
2
H8Si3O10 + H2 O (9.7)
The rate of monomer loss from solution (R, molal/sec) tends to be proportional
to the fourth power of monomer concentration (Icopini et al., 2005; Rothbaum
and Rohde, 1979). Icopini et al. (2005) explain this by assuming that the rate-
determining step is a reaction converting the trimer, formed by reaction (9.7),
into a cyclical tetramer.
If reactions (9.6) and (9.7) reach equilibrium quickly, the overall rate is propor-
tional to the concentration of H8Si3O10, which is related to the monomer con-
centration by the K1 and K2 equilibrium constants. The f(γ) term is the product
of activity coefficients needed to convert the activity of the silica species into
concentrations.
mH4 4 SiO4
R = f ( g )k3 mH8Si3 O10 mH4 SiO4 = f ( g )k3 K1K 2
aH2 2O (9.9)
dm
− = k ′m 4 (9.10)
dt
Equation (9.10) can be rearranged and integrated for the boundary condition
that when t = 0, m = mo.
m t
dm
∫m 4
= − k ∫ dt (9.11)
mo 0
1 −3
− k ′t = ( m − mo−3 ) (9.12)
3
25
monomer concentration, millimolal
20
pH 11
pH 3
15
10
5
pH 7
0
0 1000 2000 3000
time, hr
3
1
m= −1 / 3
(9.13)
3k ′t + mo
This equation, along with values of mo and k´ (Table 2 in Icopini et al., 2005),
can be used to forward model the monomer concentration (m, millimolal) as a
function of time (Figure 9.1) for an ionic strength of 0.01.
mo = 20.8 mm (millimolal) = 2.08 × 10−2 m (same as for experiments in Icopini
et al., 2005)
Nucleation rates
When a homogeneous aqueous phase becomes supersaturated with respect
to a solid phase, the rate of appearance of the new solid phase is con-
trolled by nucleation kinetics. Nucleation occurs whenever enough of the
solute species come together to create a particle that can grow spontan-
eously. Classical models of nucleation kinetics are based on a combination
186 Accretion and transformation kinetics
(2 3) bσVm 2
log K c = log K b + (9.15)
2.303RT D
The surface free energy of a cluster plays a critical role in determining its
stability. Figure 9.2a shows that surface free energy of ionic solids correlates
with the substance’s bulk solubility (Söhnel, 1982) and Figure 9.2b shows
how the solubility product of a cluster decreases with increasing cluster
diameter until it approaches the solubility product of the bulk phase. A
larger value of σ results in a higher solubility product for a particular cluster
diameter. For clusters that are larger than about 1 µm the contribution of
∆Gs to the overall free energy is negligible, so their solubility is effectively the
same as the bulk phase.
Nucleation 187
(a) 160
120
, mJ/m 2
80
40
0
–6 –4 –2 0 2
log s
(b) 25
20
Figure 9.2. (a) Correlation
barite between the interfacial free
energy (σ, mJ/m2) of ionic
15 solids and the logarithm of
their molar solubility (Ms,
Kc/Kb
Nucleation
Nucleation models postulate that clusters must attain a critical size before
they can irreversibly grow to a macroscopic size. Although details differ,
all the classical nucleation models develop the same general relationship
between the nucleation barrier, surface free energy, and degree of saturation.
The models mostly differ in how the rate of addition of a monomer to the
188 Accretion and transformation kinetics
A + A = A2 K 2 = aA2 aA2
A 2 + A = A3 K 3 = aA3 aA3
(9.16)
A n −1+ A = A n K n = aAn aAn
∆G °f ( An ) = − nRT ln aA − RT ln K n = − RT ln M An (9.18)
d ∆Gc
= − kBT ln S + (2 3)(36 π ) nA−1 3V12 3σ
13
(9.22)
dn
When d∆Gc/dn = 0, n = n* so Eq. (9.22) can be set equal to zero and solved
for n*.
kBT ln S = (2 3)(36π )
13
( n* )−1 3 V12 3σ (9.23)
n* =
(32π 3)V12σ 3 (9.24)
( kBT ln S )3
The free energy of formation of the critical nucleus (∆Gc*, J/cluster) is found
by substituting Eq. (9.24) into Eq. (9.21).
23
(32π 3)V12σ3 k T ln S + 36π 1 3 (32π 3)V12σ3
∆Gc* = − ( B ) ( ) 3
V12 3σ (9.25)
( kBT ln S )3 ( kBT ln S )
23 23
48π
(36π )1 3 (32π 3)2 3 = (6)2 3
32 192
π = (64 ) π = 16π =
23
π = (9.26)
3 3 3
∆Gc* = −
(32π 3)V12σ3 + ( 48π 3)V12σ3 = (16π 3)V12σ3 (9.27)
( kBT ln S )2 ( kBT ln S )2 ( kBT ln S )2
The nucleation rate (N, nuclei/L sec) is the product of the flux of monomers
(JA, molecules/m2 sec) to each critical nucleus and the surface area of each
critical nucleus (Ac*, m2/cluster) multiplied by the number of critical nuclei
per unit volume of solution (c*, nuclei/L).
N = J AAc*c* (9.28)
Determining the flux of monomers to the critical nuclei (JA, molecules/L sec)
is a challenge. This flux has been modeled using various statistical mechan-
ics approaches. These models estimate the flux of molecules from the solu-
tion to the surface of the critical nuclei but they include various adjustments
such as the Zeldovich factor (Markov, 2003), which accounts for the possi-
bility that some critical-sized clusters will dissolve to a smaller size instead
of growing into crystallites. Presently there is no widely accepted model for
JA, so Eq. (9.28) is usually recast into a semi-empirical form, which com-
bines the monomer flux and cluster surface area into a pre-exponential term
(Ω, crystallites/L sec).
N = Ωc* (9.29)
c* = e − ∆Gc (9.30)
* kB T
Combining Eqs (9.27), (9.29), and (9.30) produces an equation that expresses
the nucleation rate in terms of the degree of saturation.
Barite scales can clog water treatment, petroleum production, and geothermal
facilities (Boerlage et al., 2002; He et al., 1994). An additional problem is that
the barite that precipitates in these and other systems often incorporates signifi-
cant amounts of radioactive elements, especially radium. Removal and disposal
of these radioactive barite precipitates is complicated and costly (Ceccarello
et al., 2004). Schemes to mitigate barite scale formation typically use additives
to inhibit barite nucleation. The first step toward evaluating the efficacy of these
additives is to develop a model for barite nucleation rates for untreated condi-
tions. Most of these additives are chosen because they increase the surface free
energy of the barite clusters thus slowing the nucleation rate.
The following data are needed to construct a model of barite nucleation.
The first step toward developing a conceptual model of the nucleation process
uses Eq. (9.21) to map out ∆Gc as a function of log S. It is convenient to recast
this equation by inserting the appropriate numerical constants.
∆Gc = − n (1.38 × 10 −23 ) (298) 2.303 log S + 4.84 n 2 3 (8.65 × 10 −29 ) (136 × 10 −3 )
23
The maximum on each graph of ∆Gc versus n (Figure 9.3a) occurs where n = n*.
These graphs show that the ∆Gc* values associated with each maximum decrease
192 Accretion and transformation kinetics
(a) (b)
2×10 –18 1000
1.5
800
1×10–18
Gc, J/cluster
2.0 600
0
n*
2.5 400
–1×10–18 3.0
200
–2×10 –18 0
0 100 200 300 400 0 1 2 3 4 5
n, molecules log S
(c)
10
4
log N
–2
0 1 2 3 4 5 6 7
log S
Figure 9.3. (a) Free energy of formation of BaSO4 clusters as a function of the number of BaSO4 molecules in
the cluster contoured for values of log S ranging from 1.5 to 3.0. (b) Number of BaSO4 molecules in the critical
cluster as a function of log S. (c) Nucleation rate (nuclei/L sec) of BaSO4 as a function of log S.
as log S increases, which means that the barrier to nucleation decreases as the
solution becomes more supersaturated.
The size of the critical nucleus as a function of the degree of saturation is
found using Eq. (9.34).
Aggregation rates 193
n* =
(33.5)(8.65 × 10 −29 )2 (0.136)3 =
742.7
(9.34)
(1.38 × 10 ) (298) (2.303) ( log S ) ( S )3
−23 3 3 3 3
log
This equation shows that the size of the critical nucleus grows very large when
the degree of saturation drops below 10 (log S < 1) and becomes very small
when the degree of saturation exceeds 10,000 (log S > 4) (Figure 9.3b).
The barite nucleation rate as a function of the degree of saturation is calcu-
lated using Eq. (9.32).
Aggregation rates
If the concentration of polymers or crystallites is high enough they can
aggregate into larger clusters. For the simplest case, each cluster (A) moves
around in the solution as the result of Brownian motion (diffusive behav-
ior) and whenever two clusters come into contact they stick together.
This is called diffusion-limited aggregation (DLA). For clusters of equal
diameter, the concentration of A decreases following a second-order rate
equation.
dcA k
= − S cA2 (9.36)
dt 2
The rate constant for this equation can be calculated from the Smoluchowski
model.
8kBT
kS = (9.37)
3η
The Boltzmann constant (kB) equals 1.38 × 10−23 J/molecule K and at 25°C
the dynamic viscosity of water (η) equals 8.91 × 104 Pa sec (1 Pa = 1 J/m3).
194 Accretion and transformation kinetics
This constant can be inserted into Eq. (9.36) to predict the rate of disap-
pearance of A due to diffusion-limited aggregation at 25°C.
dcA clusters 1 1
= − kS cA = − 1.23 × 10 cA
−17 2
3
2
(9.39)
dt m sec 2 2
Equation (9.36) can be integrated from the initial cluster concentration (cA°)
to the time when cA = ½cA°.
2 2 2 2
− + o =− + o = − kS t1 2 (9.40)
cA cA 1 2 cA cA
Equation (9.40) can be rearranged to find the time needed to reduce the
cluster concentration to half its original value.
2
t1/ 2 = (9.41)
kS cAo
Equation (9.41) shows that cluster aggregation is very sensitive to the cluster
concentration, so that for low concentrations very long times are required
for significant aggregation to occur (Figure 9.4a) but aggregation times are
very short for suspensions with high cluster concentrations. Because dif-
fusion rates increase with increasing temperature, the aggregation times
decrease significantly with increasing temperature (Figure 9.4b). When log
kS is graphed versus 1/T (Figure 9.4c), the activation energy for aggregation
is found to be ~15.5 kJ/mol, which is consistent with a diffusion-limited
process.
Diffusion-limited aggregation (DLA) is often referred to as a fast process
because the above model predicts the maximum aggregation rate. Reaction-
limited aggregation (RLA) rates are slower than predicted by the DLA
model. The effect of the slow attachment reaction rate on the aggregation
rate is expressed as a stability ratio (W) and is defined as the ratio of the
Smoluchowski rate constant (kS) to the observed rate constant (kO).
kS
W= (9.42)
kO
(a) (b)
10
8 7000
10
7 6000
1 yr
10 6 5000
1 wk
10 5 1 day 4000
t1/2, sec
t1/2, sec
10
4 3000
1 hr
10 3 2000
2
10 1000
10 1 0
10 –8 10 –6 10 –4 10 –2 0 50 100 150 200 250 300
c°, mol/L temperature, °C
(c)
–15.5
–16.0
log kS
–16.5
–17.0
–17.5
0.002 0.003 0.004
1/T, K–1
Figure 9.4. (a) DLA model prediction of the time needed for half the clusters to aggregate as a function of clus-
ime needed for half the clusters to aggregate as a function of temperature for c° = 1 × 10−4
ter concentration. (b) T
mol/L. (c) Arrhenius plot showing log kS as a function of 1/T for a DLA process.
Axford (1997) showed that the W value for silica cluster aggregation is affected
by the ionic strength and cluster size. This information can be used to model the
effect of ionic strength and cluster size on the aggregation time.
The first step in building this model uses multiple linear regression to fit the W
data in Figure 9.5 as a function of ionic strength (I) and cluster size (N).
196 Accretion and transformation kinetics
7
N
6 1
5 10
Figure 9.5. Logarithm
of experimentally
determined stability 4 100
log W
The next step substitutes log kS into the logarithmic transform of Eq. (9.42) to
find an empirical rate constant (ko).
Combining Eqs (9.43), (9.44), and (9.45) gives a function for t½ as a function of
cluster size, ionic strength, and initial cluster concentration.
Figure 9.6 shows that t½ is very small for large clusters at high ionic strength and
very large for small clusters at low ionic strength. The silica cluster aggregation
rate is also strongly affected by pH, which controls the surface charge of the
clusters.
Solid–solid transformation rates 197
–0.1
log t = 1
2
–0.2
3
log I
4
–0.3
5
–0.4
0 1 2 3 4
log N
Figure 9.6. Logarithm of the time needed for one half of the silica clusters to
aggregate (log t½, time in seconds) contoured in terms of the logarithm of
the ionic strength (I) and logarithm of the cluster size (N = number of 12 nm
diameter units in the cluster).
Solid–solid transformation rates
The rate of transformation of a metastable solid (parent) phase (A) to form
a more stable solid (product) phase (B) is usually modeled using the Avrami
equation (Avrami, 1939, 1940), which is also known as the Johnson–Mehl–
Avrami–Kolmogorov (JMAK) equation. This equation is based on a model
that assumes that the transformation involves the nucleation of the prod-
uct phase followed by its growth until the parent phase is replaced by the
198 Accretion and transformation kinetics
d α = d α e (1 − α ) (9.47)
For the case where abundant nuclei are present in the initial stages of the
transformation, which is termed site saturation, the extended volume of one
sphere of product phase is modeled by changing the sign for the rate con-
stant in the shrinking particle model (Chapter 6) to make it a growing par-
ticle model. The linear growth rate of each sphere of the product phase is G
(m/sec) and its volume is V (m3).
4π
Vi = (Gt )3 (9.48)
3
αe =
∑V i
=N
4π
(Gt )3 (9.49)
Vt 3
d α e = N 4π (Gt ) dt
2
(9.50)
d α = (1 − α ) N 4π (Gt ) dt
2
(9.51)
α t
dα
∫ (1 − α ) = N 4πG ∫ t dt
2 2
(9.52)
0 0
N 4πG 2 3
− ln (1 − α ) = t = kt3 (9.53)
3
α = 1 − e − kt3 (9.54)
α = 1 − e − ktn (9.55)
ln ( − ln (1 − α )) = ln k + n ln t (9.56)
α = 1 − e − (kt)
n
(9.57)
ln ( − ln (1 − α )) = n ln k + n ln t (9.58)
for interpreting reaction processes from values of n found using the modi-
fied Avrami equation.
It is often useful to express the conversion in terms of a characteristic
time. The classical Avrami equation (9.55) can be solved to express t as a
function of α.
ln (1 − α )
1n
t = − (9.59)
k
Substituting 0.5 for α in Eq. (9.59) gives the time to 50% conversion.
ln (1 − 0.5)
1n 1n
0.693
t1 2 = − = (9.60)
k k
1
(− ln (1 − α ))
1n
t= (9.61)
k
Substituting 0.5 for α in Eq. (9.61) gives the time to 50% conversion.
1 1
( − ln (1 − α )) = (0.693)
1n 1n
t= (9.62)
k k
TTT diagrams
Time–temperature–transformation (TTT) diagrams are extensively used in
materials sciences to map out the characteristic time for solid–solid trans-
formations. They are sometimes called kinetic phase diagrams. These dia-
grams show the extent of transformation, α, contours on graphs with log
t on the x-axis and either T or 1/T on the y-axis. They are often prepared
by direct plotting of experimental data but they can be constructed for any
type of reaction using the integrated rate equation along with a temperature
function (Arrhenius equation) for the rate constant.
The simplest case is a TTT diagram based on a zeroth-order rate
equation.
dn
= −k (9.63)
dt
The integrated form of this equation gives the number of moles remaining
at a given time.
TTT diagrams 201
(a) 1.00
k = 10–4
n=3
0.80
k = 10–5
n=3
0.60
k = 10–4
n=2
0.40
k = 10–5
0.20 n=2
n − no = − kt (9.64)
n k
1− = 1− p = α = t (9.65)
no no
202 Accretion and transformation kinetics
ln α = ln k − ln no + ln t (9.66)
E 1
ln α = ln A − a − ln no + ln t
RT
A E 1
= ln − a + ln t (9.67)
no R T
A Ea 1
logα = log − + log t
no 2.303R T
A Ea 1
log t = − log + logα + (9.68)
no 2.303R T
E 1
ln ( − ln (1 − α )) = ln A − a + n ln t (9.69)
RT
1 E 1
ln t = − ln A + ln ( − ln (1 − α )) + a
n R T
(9.70)
1 Ea 1
log t = − log A + log ( −2.303 log (1 − α )) +
n 2.303R T
The TTT equation for the modified Avrami model (Eq. (9.58)) is slightly
different.
(1 n) ln ( − ln (1 − α )) = ln A −
Ea 1
+ ln t (9.71)
RT
E 1
ln t = − ln A + (1 n ) ln ( − ln (1 − α )) + a
RT
(9.72)
Ea 1
log t = − log A + (1 n ) log ( −2.303 log (1 − α )) +
2.303R T
(a) 200
150
goethite
= 0. 1
0.5
temperature, °C
0.9
100
50
schwertmannite
0
10 100 1000 10000 100000
time, sec
(b) 0.0022
0.0024
0.0026 goethite
= 0. 1
100° C
0. 5
0.0028 0. 9
1/T, K–1
75°
0.0030
50°
0.0032 schwertmannite
25°
0.0034
0.0036
10 100 1000 10000 100000
time, sec
Davidson et al. (2008) measured the rates of this transformation over the tem-
perature range of 60 to 190°C and fit the results to the modified Avrami equa-
tion. The average value of n from their fits was 2.1(0.25). Fitting their reported
rate constants to the Arrhenius equation gives an estimated value of Ea of 32.4
kJ/mol with log A = 1.31. These values are substituted into Eq. (9.72) to create
the TTT diagram shown in Figure 9.8.
Chapter 10
Pattern formation
So far this book has offered appetizers. This chapter deals with the prepa-
ration of the main course by introducing concepts for linking the simple
models discussed in this book to understand the complex processes that
lead to pattern formation in geological settings. Geoscientists spend much
of their time and effort identifying and explaining naturally occurring spa-
tial and temporal patterns with the goal of interpreting those patterns to
understand the processes and conditions that formed them. They are chal-
lenged by the need to identify meaningful patterns in situations that often
appear to be chaotic (Crutchfield, 2012). Meaningful patterns are frequently
subtle and recognizing them often requires clues provided by process mod-
els. Most pattern-forming systems are too complex to interpret in a holis-
tic way so our strategy is to first parse them into simple processes, which
can be accurately modeled using methods like those explained in this book.
The resulting models of discrete processes are then linked to simulate the
overall pattern-forming scenario. This strategy is widely used in science and
technology. For example, engineers design processing plants by dividing the
overall process into unit operations, each of which is responsible for a single
chemical or physical transformation of a feedstock (Gupta and Yan, 2006;
Hendricks, 2006; McCabe et al., 1993). These unit operations are modeled
separately and the models are combined to simulate the entire processing
plant. This strategy is especially effective when the unit operations occur in
a linear array of steps so that the product of one step is the feed for the next.
Natural processes are often more complicated because they can switch from
one path to another in a stochastic manner or because there is one or more
feedback loops in the overall process. Regardless of the complexity of the
situation, the divide and analyze strategy is the most effective way to under-
stand how observed patterns are related to unit processes. The challenge of
interpreting pattern-forming processes is a very exciting scientific frontier
(Ball, 1999; Nicolis and Prigogine, 1989).
Geochemists observe patterns of element and mineral distribution
in nature and then use forward and inverse models of unit processes to
understand how these patterns developed. The patterns result from the
205
206 Pattern formation
Non-equilibrium thermodynamics
Non-equilibrium thermodynamics (de Groot and Mazur, 1984; Prigogine,
1977) affords an abstract, and therefore very general, foundation for under-
standing pattern formation. The abstract nature of this approach makes
direct application to geochemical problems difficult, but non-equilibrium
thermodynamics does provide a powerful conceptual basis for thinking
about pattern formation.
Non-equilibrium thermodynamics uses the notion of rate of entropy pro-
duction (Kondepudi and Prigogine, 1998) to show that physical and chemi-
cal transformations can all be described in terms of forces and fluxes. For
example, the force driving a chemical reaction is the difference between the
chemical potential of the reactants and products. The greater this difference,
the faster the flux (rate) of conversion of reactants into products. The force
driving diffusion is the difference between a species’ chemical potentials
in two spatial regions. The diffusion flux from the region of high chemical
potential to the region of low chemical potential is directly proportional to
this force. Similar relationships exist between the rate of flow of electrical
current and voltage differences and the rate of heat flow and temperature
differences. In all these cases, a resistance limits the rate of transformation
or mass transfer. This resistance shows up in rate equations as the rate con-
stant and in diffusion equations as the diffusion coefficient divided by dis-
tance. This means that all these entropy-producing processes are analogous
at a very fundamental level.
An even more potent concept comes from the Onsager reciprocal relations,
which states that there is a coupling between conjugate force–flux pairs. For
example, mass transfer of species by diffusion in an aqueous solution causes
a change in concentration, which is accompanied by heat consumption or
release due to the heat of dilution. This sets up a thermal gradient, which
causes heat flow. The resulting link between heat flux and isothermal diffu-
sion is the Dufour effect. Its conjugate is the Soret effect, which is the diffu-
sional mass flux linked to heat flow. The Soret effect has been coupled with
Turing model 207
Turing model
The coupling of individual processes, either due to Onsager reciprocal type
relationships or due to continuity requirements, produces feedbacks between
unit processes. Feedback effects control the location and structure of many
patterns. Feedback occurs when the output of a succeeding unit process
somehow affects the behavior of an antecedent unit process. Turing (1953)
quantified the relationship between feedback and pattern formation in bio-
logical systems and the Turing model has proved useful for explaining many
kinds of biological patterns (Kondo and Miura, 2010). It has great potential
for explaining many other kinds of pattern formation such as oscillating
chemical reactions and Liesegang bands. The Briggs–Rauscher reaction,
the Bray–Liebhafsky reaction, and the Belousov–Zhabotinsky reaction are
well-known examples of a growing number of recognized oscillating chemi-
cal reactions (Epstein et al., 1983; Field and Schneider, 1989) that exhibit
cyclical temporal patterns. These same reactions can also produce combined
temporal and spatial patterns, which appear as waves that move through the
reaction medium (Ross et al., 1988; Winfree, 1972). Liesegang (1913) rec-
ognized the occurrence of spatial patterns caused by precipitation linked to
diffusion and used this idea to explain band formation in agates (Liesegang,
1915). Although agate bands are probably not a good example of this phe-
nomenon (Heaney, 1993), Liesegang bands do occur in many geological set-
tings (Augustithis et al., 1980; Fu et al., 1994) and they are relatively easy
to produce in laboratory settings (Henisch, 1988; Krug et al., 1996; Müller
208 Pattern formation
et al., 1982). Feedback processes that link mass and energy fluxes are a key
feature of pattern formation.
Geological examples
Porphyry copper ore deposits (John, 2010) are a good example of how
energy and mass flux linked by feedback produce recognizable geologic pat-
terns. Heat from the mantle produces copper-enriched magma by partial
melting in the upper mantle or at the base of the crust. Because the magma
has a lower density than the surrounding rocks, it rises carrying along this
thermal energy until it is emplaced near the Earth’s surface where cooling
and crystallization release heat to the surroundings. The crystallizing magma
also releases copper-rich hydrothermal fluids with sufficient pressure to
hydrofracture a well-defined region near the crystallizing magma. The shape
of this region is controlled by a dynamic interaction between the magmatic
fluids and the surrounding meteoric water (Weiss et al., 2012). The cool-
ing hydrothermal solutions permeate the fractures where changes in pres-
sure and temperature cause the precipitation of copper-bearing minerals
in enough abundance to form an economic deposit. The porphyry copper
process is driven by energy transport from the upper mantle to the Earth’s
surface and this energy flow does work, some of which concentrates copper.
The location of the copper-enriched zone is fixed by feedback between the
fluid flow, permeability generation, and chemical precipitation. These linked
processes produce a distinct pattern of igneous features, rock alteration, iso-
topic distribution, etc., which form the basis for the exploration models used
to discover new deposits (John, 2010). Geothermal energy drives this aspect
of the porphyry copper pattern.
A subsequent phase of porphyry copper deposit pattern formation is
driven by solar energy, via the hydrologic cycle. The hydrologic cycle trans-
ports water from the oceans to the location of a porphyry copper deposit.
The resulting erosion exposes the sulfide-rich porphyry copper ore to the
atmosphere allowing the copper sulfides to oxidize and release copper into
solution. Some of the soluble copper is dispersed by flowing surface and
groundwater to produce copper dispersal patterns that are used to locate
deposits using geochemical exploration techniques (Rose et al., 1979). Some
of the copper-bearing groundwater infiltrates to the water table, where the
redox conditions change from oxidizing to reducing so the dissolved copper
reacts with pyrite and other sulfides preserved by the reducing conditions to
precipitate minerals that form an enrichment blanket. The enrichment blan-
ket often contains very high copper values and is a key exploration target.
Similar energy + mass flow with feedback processes are responsible for
the formation of other types of ore deposits. Ore deposits are a subset of
Geological examples 209
the large number of geological and geochemical processes that are driven
by energy and mass fluxes leading to pattern formation. Examples of a few
that have been modeled include patterns in stylolites (Merino, 1992), agates
(Wang and Merino, 1995), terra rosa (Merino and Banerjee, 2008), and bur-
ial dolomites (Merino and Canals, 2011) with many more examples given in
Jamtveit and Hammer (2012), Jamtveit and Meakin (1999), and Ortoleva
(1994). Creating general models that link pattern formation and unit proc-
esses like those described in this book is an exciting scientific frontier.
References
Adamson, A.W., Gast, A.P. (1997). Physical Chemistry of Surfaces. Wiley, New York.
Alkattan, M., Oelkers, E.H., Dandurand, J.L., Schott, J. (1998). An experimental study
of calcite and limestone dissolution rates as a function of pH from −1 to 3 and
temperature from 25 to 80°C. Chemical Geology, 151, 199–214.
Amis, E.S. (1966). Solvent Effects on Reaction Rates and Mechanisms. Academic Press,
New York.
Appelo, C.A.J. (1996). Multicomponent ion exchange and chromatography in natural
systems. In Reactive Transport in Porous Media, eds. Lichtner, P.C., Steefel, C.I.
Mineralogical Society of America, Washington D.C., pp. 193–227.
Applin, K.R., Lasaga, A.C. (1984). The determination of SO42−, NaSO4−, and MgSO4o
tracer diffusion coefficients and their application to diagenetic flux calculations.
Geochimica et Cosmochimica Acta, 48, 2151–2162.
Aris, R. (1956). On the dispersion of a solute in a fluid flowing through a tube.
Proceedings of the Royal Society A, 235, 67–77.
Aris, R. (1989). Elementary Chemical Reactor Analysis. Dover Publications,
New York.
Aris, R. (1994). Mathematical Modeling Techniques. Dover Publications, New York.
Arking, A. (1996). Absorption of solar energy in the atmosphere: Discrepancy between
model and observations. Science, 273, 779–782.
Arnold, F.H., Schofield, S.A., Blanch, H.W. (1986). Analytical affinity chromatography
I. Local equilibrium and the measurement of association and inhibition constants.
Journal of Chromatography, 355, 1–12.
Arvidson, R.S., Luttgbe, A. (2010). Mineral dissolution kinetics as a function of
distance from equilibrium: New experimental results. Chemical Geology, 269,
79–88.
Asano, T., le Nobel, W.J. (1978). Activation and reaction volumes in solution. Chemical
Reviews, 78, 407–489.
Ašperger, S. (2003). Chemical Kinetics and Inorganic Reaction Mechanisms, 2nd edn.
Kluwer Academic, New York.
Astruc, D. (1995). Electron transfer and radical processes in transition-metal chemistry.
VCH Publishers, Inc., New York.
Augustithis, S.S., Mposkos, E., Vgenopoulos, A. (1980). Diffusion rings (sphaeroids) in
bauxite. Chemical Geology, 30, 351–362.
Avrahami, M., Golding, R.M. (1968). The oxidation of the sulphide ion at very low
concentrations in aqueous solutions. Journal of the Chemical Society A, 647–651.
Avrami, M. (1939). Kinetics of phase change. I. Journal of Chemical Physics, 7,
1103–1112.
210
References 211
Brezonik, P.L. (1994). Chemical Kinetics and Process Dynamics in Aquatic Systems.
Lewis Publishers, Boca Raton F.L.
Bridgman, P.W. (1931). Dimensional Analysis. Yale University Press, New Haven C.T.
Bruckner, R. (2002). Advanced Organic Chemistry Reaction Mechanisms. Harcourt
Academic Press, San Diego C.A.
Brunauer, S., Emmett, P.H., Teller, E. (1938). Adsorption of gases in multimolecular
layers. Journal of the American Chemical Society, 60, 309–319.
Bunge, H.J. (1997). Some remarks on modelling and simulation of physical phenomena.
Textures and Microstructures, 28, 151–165.
Burd, A.B., Jackson, G.A. (2009). Particle aggregation. Annual Reviews of Marine
Science, 1, 65–90.
Burkin, A.R. (2001). Chemical Hydrometallurgy. Imperial College Press, London.
Burrows, N.D., Yuwono, V.M., Penn, R.L. (2010). Quantifying the kinetics of crystal
growth by oriented aggregation. MRS Bulletin, 35, 133–137.
Burton, W.K., Cabrera, N. (1949). Crystal growth and surface structure: Part 1.
Discussions of the Faraday Society, 5, 33–39.
Burton, W.K., Cabrera, N., Frank, F.C. (1949). Role of dislocations in crystal growth.
Nature, 163, 398–399.
Burton, W.K., Cabrera, N., Frank, F.C. (1951). The growth of crystals and the
equilibrium structure of their surfaces. Philosophical Transactions Royal Society of
London, 243, 299–358.
Buxton, G.V., Greenstock, C.L., Helman, W.P., Ross, A.B. (1988). Critical review of rate
constants for reactions of hydrated electrons, hydrogen atoms and hydroxyl radicals
(•OH/•O−) in aqueous solution. Journal of Physical and Chemical Reference Data, 17,
513–886.
Buxton, G.V., Mulazzani, Q.G., Ross, A.B. (1995). Critical review of rate constants for
reactions of transients from metal ions and metal complexes in aqueous solution.
Journal of Chemical and Physical Reference Data, 24, 1055–1349.
Cabrera, N., Burton, W.K. (1949). Crystal growth and surface structure: Part 2.
Discussions of the Faraday Society, 5, 40–48.
Cabrera, N., Vermiyea, D.A. (1958). The growth of crystals from solution. In Growth
and Perfection of Crystals, eds. Doremus, R.H., Roberts, B.W., Turnbull, D. John
Wiley & Sons, New York, pp. 393–410.
Campanario, J.M. (1995). Automatic ‘balancing’ of chemical equations. Computers in
Chemistry, 19, 85–90.
Capellos, C., Beielski, B.H.J. (1972). Kinetic Systems. Wiley-Interscience, New York.
Carslaw, H.S., Jaeger, J.C. (1959). Conduction of Heat in Solids. Oxford University Press,
New York.
Casado, J., López-Quintela, M.A., Lorenzo-Barral, F.M. (1986). The initial rate method
in chemical kinetics. Journal of Chemical Education, 63, 450–452.
Casey, W.H. (2001). A view of reactions at mineral surfaces from the aqueous phase.
Mineralogical Magazine, 65, 323–337.
Casey, W.H., Swaddle, T.W. (2003). Why small? The use of small inorganic clusters to
understand mineral surface and dissolution reactions in geochemistry. Reviews in
Geophysics, 41, 1008.
References 213
Casey, W.H., Westrich, H.R., Banfield, J.F., Ferruzzi, G., Arnold, G.W. (1993). Leaching
and reconstruction at the surface of dissolving chain-silicate minerals. Nature, 366,
253–256.
Cazes, J., Scott, R.P.W. (2002). Chromatography Theory. Marcel Dekker, Inc.,
New York.
Ceccarello, S., Black, S., Read, D., Hodson, M.E. (2004). Industrial radioactive barite
scale: Suppression of radium uptake by introduction of competing ions. Minerals
Engineering, 17, 323–330.
Chaïrat, C., Schott, J., Oelkers, E.H., Lartigue, J.-E., Harouiya, N. (2007). Kinetics
and mechanism of natural fluroapatite dissolution at 25°C and pH from 3 to 12.
Geochimica et Cosmochimica Acta, 71, 5901–5912.
Chaudhry, M.H. (2008). Open-Channel Flow. Springer, New York.
Chermak, J.A., Rimstidt, J.D. (1990). Hydrothermal transformation rate of kaolinite to
muscovite/illite. Geochimica et Cosmochimica Acta, 54, 2979–2990.
Chotantarat, S., Ong, S.K., Sutthirat, C., Osathaphan, K. (2011). Effect of pH on
transport of Pb2+, Mn2+, Zn2+, and Ni2+ through lateritic soil: Column experiments
and transport modeling. Journal of Environmental Sciences, 23, 640–648.
Cölfen, H., Antonietti, M. (2008). Mesocrystals and Nonclassical Crystallization. John
Wiley & Sons, Chichester, U.K.
Colombani, J. (2008). Measurement of the pure dissolution rate constant of a mineral in
water. Geochimica et Cosmochimica Acta, 72, 5634–5640.
Crank, J. (1975). The Mathematics of Diffusion. Oxford University Press, New York.
Crutchfield, J.P. (2012). Between order and chaos. Nature Physics, 8, 17–24.
Cukierman, S. (2006). Et tu, Grotthuss! and other unfinished stories. Biochimica et
Biophysica Acta, 1725, 876–885.
Curti, E. (1997). Coprecipitation of Radionucleides: basic concepts, literature review, and
first applications. Wettingen, Switzerland, p. 107.
Curti, E. (1999). Coprecipitation of radionuclides with calcite: Estimation of partition
coefficients based on a review of laboratory investigations and geochemical data.
Applied Geochemistry, 14, 433–445.
Cussler, E.L. (2009). Diffusion, 3rd edn. Cambridge University Press, Cambridge.
Cygan, R.T., Carrigan, C.R. (1992). Time-dependent Soret transport: Applications to
brine and magma. Chemical Geology, 95, 201–212.
Damasceno, P.F., Engel, M., Glotzer, S.C. (2012). Predictive self-assembly of polyhedra
into complex structures. Science, 337, 453–457.
Danckwerts, P.V. (1953). Continuous flow systems. Distribution of residence times.
Chemical Engineering Science, 2, 1–13.
Davidson, L.E., Shaw, S., Benning, L.G. (2008). The kinetics and mechanism of
schwertmannite transformation to goethite and hematite under alkaline conditions.
American Mineralogist, 93, 1326–1337.
de Groot, S.R., Mazur, P. (1984). Non-equilibrium Thermodynamics. Dover Publications,
New York.
de Pablo, J., Casas, I., Giménez, J., Molera, M., Rovira, M., Duro, L., Bruno, J. (1999).
The oxidative dissolution mechanism of uranium dioxide. I. The effect of temperature
in hydrogen carbonate medium. Geochimica et Cosmochimica Acta, 63, 3079–3103.
214 References
Debessy, J., Pagel, M., Beny, J.-M., Christensen, H., Hickel, B., Kosztolanyi, C., Poty,
B. (1988). Radiolysis evidenced by H2–O2 and H2-bearing fluid inclusions in three
uranium deposits. Geochimica et Cosmochimica Acta, 52, 1155–1167.
Deming, W.E. (1943). Statistical Adjustment of Data. Dover Publications, New York.
Denbigh, K.G., Turner, J.C.R. (1984). Residence-time distributions, mixing, and
dispersion. In Chemical Reactor Theory: An Introduction, 3rd edn. Cambridge
University Press, New York, pp. 81–110.
Denny, M.W. (1993). Air and Water. Princeton University Press, Princeton.
Derome, D., Cathelineau, M., Lhomme, T., Cuney, M. (2003). Fluid inclusion evidence
of the differential migration of H2 and O2 in the McArthur River unconformity-type
uranium deposit (Saskatchewan, Canada). Possible role on post-ore modifications of
the host rocks. Journal of Geochemical Exploration, 78–79, 525–530.
Dietzel, M. (2000). Dissolution of silicates and the stability of polysilicic acid.
Geochimica et Cosmochimica Acta, 64, 3275–3281.
Doerner, H.A., Hoskins, W.M. (1925). Co-precipitation of radium and barium sulfates.
Journal of the American Chemical Society, 47, 662–675.
Dominé, F., Dounaceur, R., Scacchi, G., Marquaire, P.-M., Dessort, D., Pradier, B.,
Brevart, O. (2002). Up to what temperature is petroleum stable? New insights from a
5200 free radical reactions model. Organic Geochemistry, 33, 1487–1499.
Dotson, N.A., Galván, R., Laurence, R.L., Tirrel, M. (1996). Polymerization Process
Modeling. VCH Publishers, Inc., New York.
Douglas, J.F. (1969). Dimensional Analysis for Engineers. Sir Isaac Pitman & Sons Ltd.,
London.
Dove, P.M. (1994). The dissolution kinetics of quartz in sodium chloride solutions at
25°C to 300°C. American Journal of Science, 294, 665–712.
Dove, P.M., Han, N., De Yoreo, J.J. (2005). Mechanisms of classical crystal growth
theory explain quartz and silicate dissolution behavior. Proceedings of the National
Academy of Sciences, 102, 15357–15362.
Draganic, I.G. (2005). Radiolysis of water: A look at its origin and occurrence in the
nature. Radiation Physics and Chemistry, 72, 181–186.
Dria, M.A., Bryant, S.L., Schechter, R.S., Lake, L.W. (1987). Interacting precipitation/
dissolution waves: The movement of inorganic contaminants in groundwater. Water
Resources Research, 23, 2076–2090.
Drljaca, A., Hubbard, C.D., van Dldik, R., Asano, T., Basilevsky, M.V., le Nobel,
W.J. (1988). Activation and reaction volumes in solution. 3. Chemical Reviews, 98,
2167–2298.
Druschel, G.K., Hamers, R.J., Luthe IIIr, G.W., Banfield, J.F. (2003). Kinetics and
mechanism of trithionate and tetrathionate oxidation at low pH by hydroxyl
radicals. Aquatic Geochemistry, 9, 145–164.
Dunning, W.J. (1955). Theory of crystal nucleation from vapor, liquid, and solid
systems. In Chemistry of the Solid State, ed. Garner, W.E. Academic Press, Inc., New
York, pp. 159–183.
Dunning, W.J. (1969). General and theoretical introduction. In Nucleation, ed.
Zettlemoyer, A.C. Marcel Dekker, Inc., New York, pp. 1–67.
Dutrizac, J.E. (2002). Calcium sulphate solubilities in simulated zinc processing
solutions. Hydrometallurgy, 2002, 109–135.
References 215
Eckert, C.A. (1972). High pressure kinetics in solution. Annual Review of Physical
Chemistry, 239–264.
Edwards, J.O., Green, E.F., Ross, J. (1968). From stoichiometry and rate law to
mechanism. Journal of Chemical Education, 45, 381–385.
Edzwald, J.K., Upchurch, J.B., O’Melia, C.R. (1974). Coagulation in estuaries.
Environmental Science and Technology, 8, 58–63.
El-Kadi, A., Plummer, L.N., Aggarwal, P. (2011). NETPATH-WIN: An interactive user
version of the mass-balance model, NETPATH. Groundwater, 49, 593–599.
Epstein, I.R., Kustin, K., DeKepper, P., Orbán, M. (1983). Oscillating chemical
reactions. Scientific American, 248, 112–123.
Erdemir, D., Lee, A.Y., Myerson, A.S. (2009). Nucleation of crystals from solution:
Classical and two-step models. Accounts of Chemical Research, 42, 621–629.
Erdey-Grúz, T. (1974). Transport Phenomena in Aqueous Solutions. John Wiley & Sons,
New York.
Ershov, B.G., Gordeev, A.V. (2008). A model for radiolysis of water and aqueous
solutions. Radiation Physics and Chemistry, 77, 928–935.
Evans, M.G., Polanyi, M. (1935). Some applications of the transition state method
to the calculation of reaction velocities, especially in solution. Transactions of the
Faraday Society, 31 875–894.
Eyring, H. (1935). The activated complex in chemical reactions. Journal of Chemical
Physics, 3, 107–115.
Feth, J.H., Robertson, C.E., Polzer, W.L. (1964). Sources of mineral constituents in water
from granitic rocks, Sierra Nevada, California and Nevada. U.S. Geological Survey.
Field, R.J., Schneider, F.W. (1989). Oscillating chemical reactions and nonlinear
dynamics. Journal of Chemical Education, 66, 195–204.
Fokin, V.M., Zanotto, E.D., Schmelzer, J.W.P. (2010). On the thermodynamic driving
force for interpretation of nucleation experiments. Journal of Non-Crystalline Solids,
356, 2185–2191.
Ford, I.J. (2004). Statistical mechanics of nucleation: A review. Journal of Mechanical
Engineering Science, 218, 883–899.
Fournier, R.O. (1977). Chemical geothermometers and mixing models for geothermal
systems. Geothermics, 5, 41–50.
Fu, L., Milliken, K.L., Sharp Jr., J.M. (1994). Porosity and permeability variations
in fractured and Liesegang-banded Breathitt sandstones (Middle Pennsylvanian),
eastern Kentucky: Diagenetic controls and implications for modeling dual porosity
systems. Journal of Hydrology, 154, 351–381.
Fubini, B. (1998). Surface chemistry and quartz hazard. Annals of Occupational
Hygiene, 42, 521–530.
Fueno, T. (1999). The Transition State. Gordon & Breach, Amsterdam, p. 329.
Furukawa, Y., Shimada, W. (1993). Three-dimensional pattern formation during growth
of ice dendrites – its relation to universal law of dendritic growth. Journal of Crystal
Growth, 128, 234–239.
Garg, L.C., Maren, T.H. (1972). The rates of hydration of carbon dioxide and
dehydration of carbonic acid at 37°C. Biochimica et Biophysical Acta, 261, 70–76.
Garn, P.D. (1975). An examination of the kinetic compensation effect. Journal of
Thermal Analysis, 7, 475–478.
216 References
Garrels, R.M., MacKenzie, F.T. (1967). Origin of the chemical composition of some
springs and lakes. In Equilibrium Concept in Natural Water Systems, ed. Gould, R.F.
American Chemical Society, Washington D.C., pp. 222–242.
Garten, V.A., Head, R.B. (1970). Homogeneous nucleation in aqueous solution. Journal
of Crystal Growth, 6, 349–351.
Gauch Jr, H.G. (1993). Prediction, parsimony and noise. American Scientist, 81,
468–478.
Gebauer, D., Völkel, A., Cölfen, H. (2008). Stable prenucleation calcium carbonate
clusters. Science, 322, 1819–1822.
Gelhar, L.W., Welty, C., Rehfeldt, K.R. (1992). A critical review of data on field-scale
dispersion in aquifers. Water Resources Research, 28, 1955–1974.
Gibbs, G.V., Cox, D.F., Ross, N.L., Crawford, T.D., Burt, J.B., Rosso, K.M. (2005).
A mapping of the electron localization function for earth materials. Physics and
Chemistry of Minerals, 32, 208–221.
Gillespie, R.J., Robinson, E.A. (2006). Gilbert N. Lewis and the chemical bond:
The electron pair and the octet rule from 1916 to the present day. Journal of
Computational Chemistry, 28, 87–97.
Gimarc, B.M. (1974). Applications of qualitative molecular orbital theory. Accounts of
Chemical Research, 7, 384–392.
Goldich, S.S. (1938). A study in rock weathering. Geology, 46, 17–58.
Goldschmidt, V.M. (1958). Geochemistry. Clarendon Press, Oxford.
Gordon, L., Rowley, K. (1957). Coprecipitation of radium with barium sulfate.
Analytical Chemistry, 29, 34–37.
Greenberg, A.E., Clesceri, L.S., Eaton, A.D. (1992). Standard Methods for the
Examination of Water and Wastewater, 18th edn. American Public Health
Association, Washington D.C.
Grossman, R.B. (1999). The Art of Writing Reasonable Reaction Mechanisms. Springer,
New York.
Gunnarson, I., Arnórsson, S. (2003). Silica scaling: The main obstacle in efficient use of
high-temperature geothermal fluids. International Geothermal Conference, Reykjavík,
Iceland, pp. 30–36.
Gupta, A., Yan, D.S. (2006). Mineral Processing Design and Operation. Elsevier,
Amsterdam.
Gutmann, V. (1978). The Donor–Acceptor Approach to Molecular Interactions. Plenum
Press, New York.
Guy, B. (1993). Mathematical revision of Korzhinskii’s theory of infiltration
metasomatic zoning. European Journal of Mineralogy, 5, 317–339.
Hall, C., Day, J. (1977). Systems and Models: terms and basic principles. In Ecosystem
Modeling in Theory and Practice, eds Hall, C., Day, J. Wiley-Interscience, New York,
pp. 6–36.
Hall, C., Day, J., Odum, H. (1977). A circuit language for energy and matter. In
Ecosystem Modeling in Theory and Practice, eds Hall, C., Day, J. Wiley-Interscience,
New York, pp. 37–49.
Harris, R.L. (1999). Information Graphics. Oxford University Press, Oxford.
Harte, J. (1988). Consider a Spherical Cow. University Science Books, Mill
Valley, C.A.
References 217
Hartman, P., Perdok, W.G. (1955a). On the relations between structure and morphology
of crystals. I. Acta Crystallographica, 8, 49–52.
Hartman, P., Perdok, W.G. (1955b). On the relations between structure and morphology
of crystals. II. Acta Crystallographica, 8, 521–524.
Hartman, P., Perdok, W.G. (1955c). On the relations between structure and morphology
of crystals. III. Acta Crystallographica, 8, 525–529.
Harvey, M.C., Schreiber, M.E., Rimstidt, J.D., Griffith, M.A. (2006). Scorodite
dissolution kinetics: Implications for arsenic release. Environmental Science and
Technology, 40, 6709–6714.
He, S., Oddo, J.E., Tomson, M.B. (1994). The inhibition of gypsum and barite nucleation
in NaCl brines at temperatures from 25 to 90°C. Applied Geochemistry, 9, 561–567.
Heaney, P.J. (1993). A proposed mechanism for the growth of chalcedony. Contributions
to Mineralogy and Petrology, 115, 66–74.
Heizmann, J.J., Bessieres, J., Bessieres, A. (1986). Advances in kinetic models. Journal de
chimie physique, 83, 725–732.
Helfferich, F.G. (1989). The theory of precipitation/dissolution waves. AIChE Journal,
35, 75–87.
Helgeson, H.C. (1968). Evaluation of irreversible reactions in geochemical processes
involving minerals and aqueous solutions – I. Thermodynamic relations. Geochimica
et Cosmochimica Acta, 32, 853–877.
Helgeson, H.C. (1979). Mass transfer among minerals and hydrothermal solution. In
Geochemistry of Hydrothermal Ore Deposits, second edn., ed. Barnes, H.L. John
Wiley & Sons, New York, pp. 568–610.
Helgeson, H.C., Garrels, R.M., MacKenzie, F.T. (1969). Evaluation of irreversible
reactions in geochemical processes involving minerals and aqueous solutions – II.
Applications. Geochimica et Cosmochimica Acta, 33, 455–481.
Hellmann, R., Eggleston, C.M., Hochella Jr., M.F., Crerar, D.A. (1990). The formation
of leached layers on albite surfaces during dissolution under hydrothermal
conditions. Geochimica et Cosmochimica Acta, 54, 1267–1281.
Hem, J.D. (1985). Study and Interpretation of the Chemical Characteristics of Natural
Water. U.S. Geological Survey, Washington D.C.
Hendricks, D. (2006). Water Treatment Unit Processes. Taylor & Francis/CRC Press,
Boca Raton, F.L.
Henisch, H.K. (1988). Crystals in Gels and Liesegang Rings. Cambridge University
Press, Cambridge.
Higgins, G.H. (1959). Evaluation of the ground-water contamination hazard from
underground nuclear explosions. Journal of Geophysical Research, 64,
1509–1519.
Hill Jr., C.G. (1977). An Introduction to Chemical Engineering Kinetics and Reactor
Design. John Wiley & Sons, New York.
Hofmann, A. (1972). Chromatographic theory of infiltration metasomatism and its
application to feldspars. American Journal of Science, 272, 69–90.
Horne, R.L., Rodriguez, F. (1983). Dispersion in tracer flow in fractured geothermal
systems. Geophysical Research Letters, 10, 289–292.
Houston, P.L. (2001). Chemical Kinetics and Reaction Dynamics. Dover Publications,
New York.
218 References
Huminicki, D.M.C., Rimstidt, J.D. (2009). Iron oxyhydroxide coating of pyrite for acid
mine drainage control. Applied Geochemistry, 24, 1626–1634.
Huntley, H.E. (1967). Dimensional Analysis. Dover Publications, New York.
Huss Jr., A., Lim, P.K., Eckert, C.A. (1982). Oxidation of aqueous sulfur dioxide.
2. High-pressure studies and proposed reaction mechanisms. Journal of Physical
Chemistry, 86, 4229–4233.
Icopini, G.A., Brantley, S.L., Heaney, P.J. (2005). Kinetics of silica oligomerization and
nanocolloid formation as a function of pH and ionic strength at 25°C. Geochimica et
Cosmochimica Acta, 69, 293–303.
Ingold, C.K. (1969). Structure and Mechanism in Organic Chemistry. Cornell University
Press, Ithaca.
Jamtveit, B., Hammer, Ø. (2012). Sculpting of rocks by reactive fluids. Geochemical
Perspectives, 1, 341–480.
Jamtveit, B., Meakin, P. (1999). Growth, Dissolution and Pattern Formation in
Geosystems. Kluwer Academic Publishers, Dordrecht.
Jensen, W.B. (1980). The Lewis Acid-Base Concepts. John Wiley & Sons, New York.
Jin, L., Auerbach, S.M., Monson, P.A. (2011). Modeling three-dimensional network
formation with an atomic lattice model: Application to silicic acid polymerization.
Journal of Chemical Physics, 134, 134703.
John, D.A. (2010). Porphyry Copper Deposit Model. U.S. Geological Survey,
Washington D.C., p. 169.
Kashchiev, D., van Rosmalen, G.M. (2003). Review: Nucleation in solution revisited.
Crystal Research and Technology, 38, 555–574.
Kirby, C.S., Cravotta III, C.A. (2005a). Net alkalinity and net acidity 1: Theoretical
considerations. Applied Geochemistry, 20, 1920–1940.
Kirby, C.S., Cravotta III, C.A. (2005b). Net alkalinity and net acidity 2: Practical
considerations. Applied Geochemistry, 20, 1941–1964.
Kline, S.J. (1965). Similitude and Approximation Theory. McGraw-Hill Book Company,
New York.
Knapp, R.B. (1989). Spatial and temporal scales of local equilibrium in dynamic fluid-
rock systems. Geochimica et Cosmochimica Acta, 53, 1955–1964.
Kondepudi, D., Prigogine, I. (1998). Modern Thermodynamics. John Wiley & Sons,
Chichester, U.K.
Kondo, S., Miura, T. (2010). Reaction–diffusion model as a framework for
understanding biological pattern formation. Science, 329, 1616–1620.
Korzhinskii, D.S. (1970). Theory of Metasomatic Zoning. Clarendon Press, Oxford, London.
Krishnamurthy, E.V. (1978). Generalized matrix inverse approach for automatic
balancing of chemical equations. International Journal of Mathematical Education in
Science and Technology, 9, 323–328.
Krug, H.-J., Brandstädter, H., Jacob, K.H. (1996). Morphological instabilities in pattern
formation by precipitation and crystallization processes. Geologische rundschau, 85, 19–28.
Laidler, K.J. (1987a). Theories of Reaction Rates. Harper & Row, New York, pp. 80–99.
Laidler, K.J. (1987b). Chemical Kinetics, 3rd edn. Harper & Row, New York.
Lamberto, D.J., Alverez, M.M., Muzzio, F.J. (1999). Experimental and computational
investigation of the laminar flow structure in a stirred tank. Chemical Engineering
Science, 54, 919–942.
References 219
Larson, M.A., Garside, J. (1986). Solute clustering and interfacial tension. Journal of
Crystal Growth, 76, 88–92.
Lasaga, A.C. (1980a). Dynamic treatment of geochemical cycles: global kinetics. In:
Kinetics of Geochemical Processes, eds Lasaga, A.C., Kirkpatrick, R.J. Mineralogical
Society of America, Washington D.C., pp. 69–105.
Lasaga, A.C. (1980b). The kinetic treatment of geochemical cycles. Geochimica et
Cosmochimica Acta, 44, 815–828.
Lasaga, A.C. (1998a). Geochemical Kinetics. Princeton University Press,
Princeton N.J.
Lasaga, A.C. (1998b). Kinetic Theory in Earth Sciences. Princeton, University Press,
Princeton, N.J.
Lasaga, A.C., Berner, R.A. (1998). Fundamental aspects of quantitative models for
geochemical cycles. Chemical Geology, 145, 161–175.
Lasaga, A.C., Blum, A.E. (1986). Surface chemistry, etch pits and mineral-water
reactions. Geochimica et Cosmochimica Acta, 50, 2363–2379.
Lasaga, A.C., Gibbs, G. (1990). Ab-initio quantum mechanical calculations of water-
rock interactions: Adsorbtion and hydrolysis reactions. American Journal of Science,
290, 263–295.
Laub, R.J. (1985). Theory of chromatography. In Inorganic Chromatographic Analysis,
ed. Macdonald, J.C. John Wiley & Sons, New York, pp. 13–186.
Le Caër, S. (2011). Water radiolysis: Influence of oxide surfaces on H2 production under
ionizing radiation. Water, 3, 235–253.
Lefticariu, L., Pratt, L.A., LaVern, J.A., Schimmelmann, A. (2010). Anoxic pyrite
oxidation by water radiolysis products – A potential source of biosustaining energy.
Earth and Planetary Science Letters, 292, 57–67.
Leifer, A. (1988). The Kinetics of Envrionmental Aquatic Photochemistry. American
Chemical Society, Washington D.C.
Lerman, A. (1979). Geochemical Processes: Water and Sediment Environments. Wiley-
Interscience, New York.
Lerman, A., Wu, L. (2008). Kinetics of global geochemical cycles. In Kinetics of Water-
Rock Interaction, eds Brantley, S.L., Kubicki, J.D., White, A.F. Springer, New York,
pp. 655–736.
Leubner, I.H. (2000). Particle nucleation and growth models. Current Opinion in Colloid
and Interface Science, 5, 151–159.
Levenspiel, O. (1972a). Chemical Reaction Engineering. John Wiley & Sons,
New York.
Levenspiel, O. (1972b). Single Ideal Reactors, Chemical Reaction Engineering, 2nd edn.
John Wiley & Sons, New York, pp. 91–117.
Levich, V.G. (1962). Physicochemical Hydrodynamics. Prentice-Hall, Inc., Englewood
Cliffs, N.J.
Lewis, G.N. (1916). The atom and the molecule. Journal of the American Chemical
Society, 38, 762–785.
Li, Y.-H., Gregory, S. (1974). Diffusion of ions in sea water and in deep-sea sediments.
Geochimica et Cosmochimica Acta, 38, 703–714.
Liesegang, R.E. (1913). Geologische Diffusionen. Verlag von Theodor Steinkopff,
Dresden.
220 References
Liesegang, R.E. (1915). Die Achate. Verlag von Theodor Steinkopff, Dresden.
Limerinos, J.T. (1970). Manning Coefficient from Measured Bed Roughness in Natural
Channels, U.S. Geological Survey Water Supply Paper 1898-B.
Lin, L.-H., Slater, G.F., Lollar, B.S., Lacrampe-Couloume, G., Onstott, T.C. (2005). The
yield and isotopic composition of radiolytic H2, a potential energy source for the
deep subsurface biosphere. Geochimica et Cosmochimica Acta, 69, 893–903.
Liu, J., Aruguete, D.M., Jinschek, J.R., Rimstidt, J.D., Hochella Jr., M.F. (2008). The
non-oxidative dissolution of galena nanocrystals: Insights into mineral dissolution
rates as a function of grain size, shape, and aggregation state. Geochimica et
Cosmochimica Acta, 72, 5984–5996.
Liu, L., Guo, Q.-X. (2001). Isokinetic relationship, isoequilibrium relationship, and
enthalpy–entropy compensation. Chemical Reviews, 101, 673–695.
Liu, S.-T., Nancollas, G.H. (1971). The kinetics of dissolution of calcium sulfate
dihydrate. Journal of Inorganic and Nuclear Chemistry, 33, 2311–2316.
Lowell, S., Shields, J.E. (1991). Powder Surface Area and Porosity, 3rd edn. Chapman &
Hall, London.
MacInnis, I.N., Brantley, S.L. (1993). Development of etch pit size distributions on
dissolving minerals. Chemical Geology, 105, 31–49.
Maffezzoli, A., Kenny, J.M., Torre, L. (1995). On the physical dimensions of the Avrami
constant. Thermochimica Acta, 269/270, 185–190.
Mandel, J. (1964). The Statistical Analysis of Experimental Data. Dover Publications,
New York.
Marangoni, A.G. (1998). On the use and misuse of the Avrami equation in
characterization of the kinetics of fat crystallization. Journal of the American Oil
Chemists’ Society, 75, 1465–1467.
Marcus, R.A. (1964). Chemical and electrochemical electron-transfer theory. Annual
Reviews of Physical Chemistry, 155–196.
Marcus, R.A. (1968). Theoretical relations among rate constants, barriers, and Brønsted
slopes of chemical reactions. Journal of Physical Chemistry, 72, 891–898.
Marcus, R.A. (1985). Electron transfers in chemistry and biology. Biochimica et
Biophysica Acta, 811, 265–322.
Marcus, R.A. (2000). Tutorial on rate constants and reorganization energies. Journal of
Electroanalytical Chemistry, 483, 2–6.
Marin, G.B., Yablonsky, G.S. (2011). Kinetics of Chemical Reactions. Wiley-VCH,
Weinheim, Germany.
Markov, I.V. (2003). Crystal Growth for Beginners. World Scientific, London.
McCabe, W.L., Smith, J.C., Harriott, P. (1993). Unit Operations of Chemical
Engineering. McGraw-Hill, Inc., New York.
McClamroch, N.H. (1980). State Models of Dynamic Systems. Springer-Verlag,
New York.
McIntire, W.L. (1963). Trace element partition coefficients: A review of theory and
applications to geology. Geochimica et Cosmochimica Acta, 27, 1209–1264.
Meldrum, F.C., Sear, R.P. (2008). Now you see them. Science, 322, 1802–1103.
Merino, E. (1992). Self-organization in stylolites. American Scientist, 80, 466–473.
Merino, E., Banerjee, A. (2008). Terra rosa genesis, implications for karst, and eolian
dust: A geodynamic thread. Journal of Geology, 116, 62–75.
References 221
Ohlin, C.A., Villa, E.M., Rustad, J.R., Casey, W.H. (2010). Dissolution of insulating
oxide materials at the molecular scale. Nature Materials, 9, 11–19.
Oreskes, N., Shrader-Frechette, K., Belitz, K. (1994). Verification, validation, and
confirmation of numerical models in the Earth Sciences. Science, 263, 641–646.
Ortoleva, P.J. (1994). Geochemical Self-Organization. Oxford University Press, New York.
Overton, W. (1977). A strategy of model construction. In Ecosystem Modeling in
Theory and Practice: An introduction with case histories, eds, Hall, C., Day, J. Wiley-
Interscience, New York, pp. 50–73.
Palandri, J.L., Kharaka, Y.K. (2004). A compilation of rate parameters of water-mineral
interaction kinetics for application to geochemical modeling. U.S. Geological Survey,
Menlo Park, C.A., p. 64.
Pauling, L. (1960). The Nature of the Chemical Bond and the Structure of Molecules and
Crystals: An Introduction to Modern Structural Chemistry. Cornell University Press,
Ithaca, N.Y. (3rd edition).
Penn, R.L. (2004). Kinetics of oriented aggregation. Journal of Physical Chemistry B,
108, 12707–12712.
Penn, R.L., Tanaka, K., Erbs, J. (2007). Size dependent kinetics of oriented aggregation.
Journal of Crystal Growth, 309, 97–102.
Petrovich, R. (1981a). Kinetics of dissolution of mechanically comminuted rock-
forming oxides and silicates – I. Deformation and dissolution of quartz under
laboratory conditions. Geochimica et Cosmochimica Acta, 45, 1665–1674.
Petrovich, R. (1981b). Kinetics of dissolution of mechanically comminuted rock-
forming oxides and silicates – II. Deformation and dissolution of oxides and silicates
in the laboratory and at the Earth’s surface. Geochimica et Cosmochimica Acta, 45,
1675–1686.
Pina, C.M., Putnis, A. (2002). The kinetics of nucleation of solid solutions from
aqueous solutions: A new model for calculating non-equilibrium distribution
coefficients. Geochimica et Cosmochimica Acta, 66, 185–192.
Piscitelle, L.J. (1990). Determination of initial rates for a general class of chemical
reactions: A methodology. International Journal of Chemical Kinetics, 22, 683–688.
Platten, J.K., Bou-Ali, M.M., Dutrieux, J.F. (2003). Enhanced molecular separation
in inclined thermogravitational columns. Journal of Physical Chemistry B, 107,
11763–11767.
Plummer, L.N., Parkhurst, D.L., Thorstenson, D.C. (1983). Development of reaction
models for ground-water systems. Geochimica et Cosmochimica Acta, 47, 665–686.
Plummer, L.N., Prestemon, E.C., Parkhurst, D.L. (1991). An interactive code
(NETPATH) for modeling net geochemical reactions along a flow path. U.S.
Geological Survey, Reston, V.A., p. 227.
Plummer, L.N., Prestemon, E.C., Parkhurst, D.L. (1992). NETPATH: An interactive
code for interpreting NET geochemical reactions from chemical and isotopic data
along a flow PATH. In 7th International Symposium on Water-Rock Interaction–
WRI-7, eds Kharaka, Y.F., Maest, A.S. A.A. Balkema, Park City, U.T., pp.
239–242.
Pollard, J.H. (1977). A Handbook of Numerical and Statistical Techniques. Cambridge
University Press, Cambridge.
Powers, J.E., Wilke, C.R. (1957). Separation of liquids by thermal diffusion. A.I.Ch.E.
Journal, 3, 231–222.
References 223
Wagner, W., Pruß, A. (2002). The IAPWS Formulation 1995 for the thermodynamic
properties of ordinary water substance for general and scientific use. Journal of
Physical and Chemical Reference Data, 31, 387–535.
Wainer, H. (2005). Graphic Discovery. Princeton University Press, Princeton, NJ.
Waley, S.G. (1981). An easy method for the determination of initial rates. Biochemistry
Journal, 193, 1009–1012.
Walton, A.G. (1969). Nucleation in liquids and solutions. In Nucleation, ed.
Zettlemoyer, A.C. Marcel Dekker, Inc., New York, pp. 225–307.
Wang, Y., Merino, E. (1995). Origin of fibrosity and banding in agates from flood
basalts. American Journal of Science, 295, 49–77.
Watson, J.T.R., Basu, R.S., Sengers, J.V. (1980). An improved representation equation
for the dynamic viscosity of water substance. Journal of Physical and Chemical
Reference Data, 9, 1255–1290.
Weber Jr., W.J., DiGiano, F.A. (1996). Process Dynamics in Environmental Systems. John
Wiley & Sons, Inc., New York.
Wechsler, J. (1988). On aesthetics in science. In Design Science Collection, ed. Loeb, A.L.
Birkhäuser, Boston, p. 180.
Wehrli, B. (1989). Monte Carlo simulations of surface morphologies during mineral
dissolution. Journal of Colloid and Interface Science, 132, 230–242.
Weinstein, L., Adam, J.A. (2008). Guesstimation. Princeton University Press, Princeton, N.J.
Weiss, P., Driesner, T., Heinrich, C.A. (2012). Porphyry-copper ore shells form at stable
pressure-temperature fronts within dynamic fluid plumes. Science, 338, 1613–1616.
Weissbart, E.J., Rimstidt, J.D. (2000). Wollastonite: Incongruent dissolution and leached
layer formation. Geochimica et Cosmochimica Acta, 64, 4007–4016.
Wen, C.Y. (1968). Noncatalytic heterogeneous solid fluid reaction models. Industrial and
Engineering Chemistry, 60, 34–54.
Weng, P.F. (1995). Silica scale inhibition and colloidal silica dispersion for reverse
osmosis systems. Desalination, 103, 59–67.
White, A.F., Peterson, M.L. (1990). Role of reactive-surface area characterization in
geochemical kinetic models. In ACS Symposium Series 416, Chemical Modeling of
Aqueous Systems II, eds Melchior, D.C., Bassett, R.L. American Chemical Society,
Los Angeles, CA, pp. 461–475.
Williamson, M.A., Rimstidt, J.D. (1993). The rate of decomposition of the ferric-
thiosulfate complex in acidic aqueous solutions. Geochimica et Cosmochimica Acta,
57, 3555–3561.
Williamson, M.A., Rimstidt, J.D. (1994). The kinetics and electrochemical rate-
determining step of aqueous pyrite oxidation. Geochimica et Cosmochimica Acta, 58,
5443–5454.
Winfree, A.T. (1972). Spiral waves of chemical activity. Science, 175, 634–636.
Wolff, G.A., Gualtieri, J.G. (1962). PBC vector, critical bond energy ratio and crystal
equilibrium form. American Mineralogist, 47, 562–584.
Wulff, G. (1977). On the question of the rate of growth and dissolution of crystal faces.
In Crystal Form and Structure, ed. Schneer, C.J. Dowden, Hutchinson & Ross, Inc,
Stroudsburg, P.A, pp. 43–52.
Yoneawa, C., Tanaka, Y., Kamioka, H. (1996). Water-rock reactions during gamma-ray
irradiation. Applied Geochemistry, 11, 461–469.
References 227
Zhang, J.-W., Nancollas, G.H. (1990). Mechanisms of growth and dissolution of sparing
soluble salts. In Mineral-Water Interface Geochemistry, eds Hochella, M.F., White,
A.F. Mineralogical Society of America, Washington D.C., pp. 365–396.
Zhang, J.-Z., Millero, F.J. (1993). The products from the oxidation of H2S in seawater.
Geochimica et Cosmochimica Acta, 57, 1705–1718.
Zheng, C., Bennett, G.D. (1995). Applied Contaminant Transport Modeling. Van
Nostrand Reinhold, New York.
Zhu, C., Anderson, G. (2002). Environmental Applications of Geochemical Modeling.
Cambridge University Press, Cambridge.
Zhu, C., Lu, P. (2009). Alkali feldspar dissolution and secondary mineral precipitation
in batch systems: 3. Saturation states of product minerals and reaction paths.
Geochimica et Cosmochimica Acta, 73, 3171–3200.
Zlokarnik, M. (1991). Dimensional Analysis and Scale-up in Chemical Engineering.
Springer-Verlag, New York.
Index
228
Index 229