A comprehensive dissolution rate theory that integrates individual surface reactions into an overall rate is developed. The dissolution theory is based on the movement of dissolution stepwaves stemming from surface defects. The net bulk rate associated with dissolution stepwaves arises quite naturally from the equations describing the spreading of the train of steps from surface defects. The overall rate can be shown to approach a simple linear rate or transition-state theory-like equation far from equilibrium. However, one of the most important results is the strong nonlinear decrease in the rate as equilibrium conditions are approached, as is the case in most natural processes. The model is validated by extensive Monte Carlo simulations of crystal dissolution, which include a detailed treatment of surface defect energetics, adsorption, surface diffusion, transport of elements from solution, and the bonding dependence of detachment processes from the surface. Monte Carlo results show the generation of dissolution stepwaves and the nonlinear dependence of the overall rate on the saturation state. The final rate equations are consistent with both the far-from-equilibrium experimental work and several recent studies that approached equilibrium. The decrease in the rate as equilibrium is approached has far-reaching implications for both man-made problems (e.g., radioactive waste disposal, pollution, etc.) and natural processes from ground-water to metamorphic systems.

crystal dissolution, model, Monte Carlo simulation, vertical scanning interferometry, transition-state theory.

Recent concern with global change, pollution, and ground-water management has led to an intense interest in the science of crystal dissolution. Licenses for radioactive waste repository sites require extrapolation of kinetics to time scales of up to 10,000 years. In general, any quantitative treatment involving fluid flow and reactions with mineral surfaces needs to extract a rate law from experimental data, which should include the possibility of approaching mineral-fluid equilibrium. Although many computer “codes” routinely input such “rate laws” (Lichtner et al., 1996), the experimental studies of mineral-fluid reactions have not provided definitive results in most cases except for conditions that are quite far from equilibrium. The existing data on bulk phases usually have been obtained by using macroscopic techniques such as batch, stirred fluid-flow, or column reactors (e.g., Brantley & White, 1995).

In recent years, bulk kinetic data have fueled a formulation of the dissolution rate based on the formation of surface molecular complexes by surface reaction with solution species (Brady & Walther, 1989; Stumm & Wieland, 1990; Stumm & Wollast, 1990; Blum & Lasaga, 1991; Gautier et al., 1994; Oelkers et al., 1994; Lasaga, 1995; Casey & Ludwig, 1995; Dove, 1995; Blum & Stillings, 1995; Xiao & Lasaga, 1994). This shift to the molecular scale has been very fruitful and important.

The molecular approach highlights the chemical bonding properties of mineral surfaces as well as the role of adsorption on the overall mineral dissolution rate. The pivotal role of surface groups such as Si-OH, Si-OH2+, Al-OH, Al-OH2+, Si-O− in alumino-silicate minerals, for example, has been recognized and used to predict the pH dependence of the rate (Carroll-Webb & Walther, 1988; Brady & Walther, 1989; Lasaga & Gibbs, 1990; Oelkers et al., 1994; Dove, 1995; Ganor et al., 1995). However, the molecular approach needs to be linked with the experimental measurement of bulk kinetic data by many-body statistical mechanical formulation of the dissolution (or growth) process. Such a link has been largely lacking in the mineral-fluid reaction field.

Many of the reaction rates for common minerals in contact with solutions are quite slow. The mean lifetimes of 1-mm mineral grains at 25°C and pH 5 can range from 103 to 107 years (Lasaga et al., 1994). The dissolution or growth of these minerals is usually determined by the dynamics of fairly euhedral mineral surfaces (Eggleston & Hochella, 1992; Hochella & Banfield, 1995; Hochella, 1995). SEM studies of mineral surfaces from the highly reactive environment of modern soils show crystallographically controlled etch pits as well as relatively flat crystal faces (Berner & Schott, 1982; White, 1995; Brantley & Chen, 1995). In general, one does not see amorphous and structure-less mineral surfaces. The existence of well-defined mineral faces, as well as steps and spirals on the surface, has been verified also using atomic force microscopy (Gratz & Bird, 1993; Dove & Hochella, 1993; Bosbach, & Rammensee, 1994). If crystallographic features have not been demolished in the dissolution process, then a molecular approach is too limited to explain the overall breakdown of a crystal structure or a bulk dissolution rate constant.

Because crystal dissolution data are obtained far from equilibrium, the treatments of crystal dissolution have taken three different approaches: (1) treat the dissolution rate as independent of the saturation state, ΔG (the distilled water case), (2) treat systems as reaching near-equilibrium conditions quickly (assumed in many field studies), or (3) treat the rate law as a simple linear relationship between rate and deviation from equilibrium (e.g., ΔG), at least close to equilibrium. In case (1), the research has then limited the scope of investigation to the variation of the dissolution rate with variables such as ionic strength, individual components such as H4SiO4 or Al3+,pH and temperature. This research has emphasized the formation of molecular “surface complexes”, which are part of the elementary reactions that would influence the overall statistical dissolution process (Furrer & Stumm, 1986; Schindler & Stumm, 1987; Wieland et al., 1988; Blum & Lasaga, 1988; Xiao & Lasaga, 1994; Oelkers, 2001). Nonetheless, these studies have ignored both the crystal (surface topography, defects, steps, crystallographic variables, etc.) and the ultimate approach (or lack thereof) of man-made or natural systems to equilibrium.

In case (3), the most often invoked relationship has been based on the principle of detailed balancing or a transition-state theory (TST) approach and leads to the rate law (Aagaard & Helgeson, 1982; Casey, 1995; Lasaga, 1998)
where A is a general constant, which could vary with pH, T, inhibitor molecules, etc., and σ should be 1 if ΔG is based on 1 mol of the rate-limiting component. Equation (1) is derived for elementary reactions, where the same activated complex controls the kinetics of the forward and reverse reactions (e.g., microscopic reversibility). Almost no justification for equation (1) has been attempted in the mineral dissolution literature on the basis of crystal properties and a surface reaction mechanism.

To apply this equation to mineral dissolution is tantamount to treating the entire crystal structure as one molecule, ignoring all other crystal surface features. At best, this formula makes the stringent assumption that the dissolution involves very rough surfaces and far-from-equilibrium conditions (Lasaga, 1998). Yet, the formulation is almost always assumed especially as near-equilibrium conditions are approached. However, mineral phases in laboratory experiments as well as in soils and weathering rock rarely exhibit this featureless and amorphous characteristic, such as demanded by a truly TST model.

The dissolution stepwave model (DSM) has arisen from new surface topography data for minerals like anorthite, dolomite, and calcite measured with vertical scanning interferometry (VSI) (Lüttge et al., 1999, 2002; Lasaga & Lüttge, 2001; Arvidson et al., 2003). Current VSI techniques provide surface topography data with an angstrom-to-nanometer vertical resolution and 0.5-microns lateral resolution. At the same time, the noninvasive optical techniques provide a large field of view of almost 1 mm2 per scan. For any set of experimental conditions, intrinsic rate data of the dissolution or growth process can be obtained by measuring the surface normal retreat as a function of time.

Major progress was made when it became possible to introduce reference surfaces within the optical field of view that allow us to conduct absolute height measurements (Lüttge et al., 1999). The application of this technique produced a surprising result (see Fig. 1): The VSI studies showed not only the formation and development of ubiquitous etch pits but also a far more important phenomenon. After some reaction time, the surface spots that had been protected against the dissolution attack formed plateaus, indicating that the entire mineral surface had retreated significantly, leading to what we may term “global dissolution” (Lüttge et al., 1999, 2003; Lasaga & Lüttge, 2001; Arvidson et al., 2003). This global dissolution can be quantified by measuring the height of the plateaus as a function of time. The VSI studies cited above have indicated that the global dissolution is fairly uniform and systematic over the entire mineral surface. For example, measurements of the surface drop for 60 different regions of three different surfaces in dolomite produced the same drop within 25% (Lüttge et al., 2003). The development of etch pits locally, or even the coalescence of etch pits, cannot account for the dissolution of the entire surface. Therefore, a universal process must govern the bulk crystal dissolution.

As presented first in Lasaga & Lüttge (2001) and developed quantitatively below, the VSI data can be understood by looking for etch pit effects that are far removed from the localized formation of dissolution holes. The importance of etch pits in the overall dissolution process stems not from the localized formation of a hole but rather from their ability to generate a continual sequence of steps. One can think of the outskirts of an etch pit as consisting of the locations where the mineral surface first drops significantly from the pre-etch pit level. At the outskirts of the etch pit, therefore, one always finds steps. These steps, however, are not static but can move into the rest of the mineral surface. After one step moves far from the center of the dislocation-generated etch pit, another step can then be generated at the outskirts of the etch pit. The continuous movements of these steps into the crystal surface produces a train of steps, which we have termed dissolution stepwaves (Lasaga & Lüttge, 2001). These dissolution stepwaves can travel throughout the mineral surface and eventually control the overall bulk dissolution rate. The stair-like train of steps in the dissolution stepwave mechanism is reminiscent of the spiral-based BCF theory (Burton et al., 1951). However, the origin, spacing and movement of the steps are completely different from the situation of a developing growth or dissolution spiral.

The overall dissolution stepwave model, then, is as shown in Fig. 2 The surface is initially flat and pitless, as in many experimental studies. If one were to analyze the crystal surface at a later time, etch pits would be visible surrounded by fairly flat regions of the surface. Without an absolute reference frame, this topography would suggest that the dissolution is all coming from the localized etch pits. However, throughout the dissolution process, stepwaves from each of the etch pits have been combining (Fig. 2), lowering the entire mineral surface. As a result, the actual topography would look like that in Fig. 2, exhibiting a uniform surface drop if a reference surface were available (as in the interferometry investigations).

In natural crystals, many dislocations are outcropping on the surface, and they produce many etch pits during dissolution (Berner & Schott, 1982; White, 1995). The general formulation of the rate constant is based on the velocity and separation of the various steps generated by each etch pit. If the mineral surface had only one etch pit site, the steps would move out like waves on the surface. The dissolution at any point on the surface would depend on the passage of each wave. If we think of each step as having a step height of Δh, then after the passage of the first wave, the surface will move down by Δh. After the second wave passes our spot, the surface will move down another Δh. In general, at time t, the surface will be at position:
where Nwaves is the number of stepwaves that have passed through. In cases of multiple etch pits, the convergence of two dissolution stepwaves merely annihilates both waves (i.e., Nwaves still depends on one etch pit). The only role, therefore, that etch pit density would play is a transient one and arises from its influence on the time needed for the first wave to travel to any one site on the surface and establish a steady-state condition there. As a consequence, the dissolution stepwave model predicts that the bulk dissolution rate will depend only very slightly on the etch pit density but quite strongly on the saturation state, as discussed below. This prediction is quite in accord with several experimental investigations and resolves a longstanding problem in water-rock kinetics. Several earlier studies have shown that for a wide variety of crystal structures, the dissolution rate is independent or at most very slightly depends on the etch pit density (Casey et al., 1988; Cygan et al., 1989; Schott et al., 1989; Blum et al., 1990). The etch pit density in these studies was varied by five orders of magnitude, yet the overall rate of dissolution varied by less than a factor of two. In a simple model based directly on the dissolution and growth of local etch pits, such a result would not be predicted until the etch pits had thoroughly coalesced to produce a very rough surface, unlike the surfaces seen after experiments or in soils.

The Monte Carlo (MC) method is an excellent tool to investigate details involving complex statistical processes of many-body phenomena such as dissolution or growth (Gilmer, 1976, 1980; Lasaga & Blum, 1986; Blum & Lasaga, 1987). We have carried out MC calculations of growth and dissolution, including the possible presence of steps and of defect-related stress fields using techniques similar to those in Blum and Lasaga (1987). Fig. 3 shows a MC simulation that outlines both the formation of an etch pit and the movement of associated dissolution stepwaves (see Appendix for discussion of MC). The bond strength, ϕ/kT has a value 6 (typical values of ϕ/kT are 4 to 8; see Blum & Lasaga, 1987, and Appendix). Δμ/kT in the figure corresponds to the free energy of the dissolution reaction (normalized to 1 mol of molecular units), so that Δμ = 0 at equilibrium.

Dislocation defects induce a strain energy in crystals, which is the energetic driving force for etch pit formation. The dislocation strain energy for most crystals can be accurately obtained from equation (3) (Heinisch et al., 1975):
where μ is the isotropic bulk shear modulus (in ergs/cm3), b is the size of the Burgers vector (Å), rh (Å), the size of the hollow core, fixes the upper bound on the strain energy (e.g., u at r = 0), and r is the distance from the dislocation center (Å) (Lasaga & Lüttge, 2001). For typical conditions as given in Fig. 3, the resulting etch pit exhibits both the formation of a deep hole and the spreading of stepwaves from the boundaries of the etch pit. The spacing of the stepwaves increases as they move away from the center of the etch pit. Note that, as in the VSI experiments, these calculations can reference the original surface. As shown in Fig. 3a, the entire surface has moved down due to the stepwave-driven dissolution, just as observed for the natural materials.

Fig. 3b displays the release of molecular units versus time (in 1/v units); see Appendix. The amount released excludes the region around the hollow core; thus, a region 19 × 19 molecular units centered on the dislocation is ignored in computing the amount dissolved. There is an initial period that exhibits a quadratic increase of the amount dissolved with time. This period is actually another manifestation of the phenomenon of reactive surface area. Initially, the whole surface is flat and of low reactivity. However, once the first dissolution front is initialized, it takes some time before the wave reaches the further regions of the crystal (just like a wave in a pond; see also Fig. 4). In our case, the surface region consists of an area 80 × 80 molecular units in size, with periodic boundary conditions. The “reactive surface area” is the area inside the envelope of the outermost dissolution front (which now contains a sequence of steps). This reactive area will increase as the square of time, if the velocity of the front is reasonably constant. Although this is only an approximation, as we will show, the overall effect is still a fairly quadratic behavior. Once the dissolution waves reach those from the periodic etch pits, the system can settle into a steady state. At that point, a fairly straight curve is achieved, as shown in Fig. 3b The slope of the curve gives directly the dissolution rate of 3.741 × 10−7 blocks/site/(1/v) (see Appendix for units).

The MC results can be used to further analyze the dissolution stepwaves by using contouring programs. Fig. 4a shows some contour plots for dissolution under different conditions: ϕ/kT = 4, Δμ/kT = -0.5 and defect conditions as given in the figure legend. Notice the close spacing of the waves (contours) close to the etch pit and the expansion of the waves as they move away from the dislocation center. This process continues inexorably and produces dissolution waves that can travel across the whole crystal surface. This continuous production of stepwaves is readily seen in the detailed dynamics shown in Fig. 4b,Fig. 4c, and d give the position of some stepwaves as a function of time for conditions corresponding to those used in Fig. 3b Note that the velocity increases initially and then approaches a fairly constant value, 3.68 × 10−7 blocks/(1/v). Fig. 4e shows the spacing of the dissolution stepwaves. Initially, the stepwaves are close together, but after some time, the spacing increases to a steady-state level around 5.2 blocks. Therefore, far from the dislocation center, the dissolution waves approach a constant velocity and spacing. Just as with the BCF theory, these results can be converted into an overall rate. The time it takes the dissolution waves to traverse the spacing is 5.2/(3.68 × 10−7) or 1.41 × 107 (1/v). Therefore, because the spacing of the stepwave contours is 5 molecular units, in this case,
\[Rate\ =\ 5/1.41\ {\times}\ 10^{{-}7}\ =\ 3.55\ {\times}\ 10^{{-}7}\ \mathrm{blocks}/\mathrm{sit}e/(1/v),\]
which agrees closely with the previous statistically averaged rate, 3.74 × 10−7 blocks/site/ (1/v) deduced from Fig. 3b 
The dissolution waves originating at crystal defects, such as etch pits, will be closely approximated by circular stepwaves. The essential treatment of these dissolution waves, therefore, must begin with the statistical physics of curved steps. In classical crystal growth theory, a curved step moves slower than a straight step, because the extra energy of the curve gives rise to a higher solubility around the curved step. A similar situation should arise in the dissolution of steps. However, certain important modifications must be included. First the Gibbs-Thomson equation, which classically relates the size of a grain to its solubility (Lasaga, 1998), needs to be altered to include the strain field of dislocation defects, u(r). Consider a circular dissolution step with depth Δh and radius r. Imagine that the radius of the dissolution step is increased by dissolving a molecular unit from the crystal. The free energy change for the dissolution of one molecule at a distance r from a surface defect is given by
where al, as are the activities of the component in solution and in the solid, respectively, σ is the surface free energy, u(r) is the strain energy at radius r, and ΔA, ΔV are the changes in area and volume on dissolution. For our case, ΔV = v̄, the molecular volume, and ΔA = v̄/r, thus, at equilibrium
or
The variation of the concentration of adatoms (adsorbed atoms, i.e., atoms with one bond to the surface) on the surface, Cs, with distance y from the step is obtained by solving the appropiate diffusion equation. Recalling that at any distance y there are attachment and detachment reactions taking place in addition to diffusion, the equation to solve is
where CB is the concentration in solution, and steady state has been assumed in equation (7). KBS and kSB are the attachment and detachment rate constants for adatoms, respectively. In the case of a series of steps with spacing L, the solution of equation (7) is given by
where CS0 and are the concentrations at the step (y = 0) and far from a step (y = L/2), respectively, and x2s = Ds/kSB, is the mean surface diffusion distance. From equation (8), we can calculate the adatom flux at a step using the gradient of CS at the step and thereby the velocity of the step
The molecular volume v̄ is heeded to convert the molar flux (mol/cm/sec) to a volume flux (cm3/cm/sec). For a unit length of step, this volume change represents a volume of 1 Δh Δy, where Δh is the step height and Δy is the horizontal movement of the step. Hence, the movement of the step is given by 1/Δh v̄Flux. Finally, the factor of two arises because the diffusion flux arrives or leaves at each step from both directions.
The term (CS∞ - CS0) is also proportional to the undersaturation. We can rewrite this term as
\[C_{S{\infty}}\ -\ C_{S0}\ =\ \frac{k_{BS}}{k_{SB}}(C_{B}\ -\ C_{B0})\]
where CBO is the concentration of solute in solution adjacent to the steps and kinks. One may expect that CBO will be small and close to the equilibrium value, CBeq. Therefore, we may replace the above by
The main modification in the case of the circular dissolution stepwaves is the change in solubility associated with both the curvature and the strain energy. This modification affects the equilibrium concentration, CeqB, and is a function of the radius, r. Labeling, Ceq∞B the straight edge (usual) solubility, CeqB(r), the new equilibrium concentration at radius r, and vstep the velocity of a straight step, then the velocity of a spherical dissolution step in our dissolution wave is:
Inserting the expression for Keq(r) (equn. (6) and using the thermodynamic relation CB/Ceq,∞B = exp(ΔG/kT), this last expression can be further rewritten as:
Equation (13) is the relevant equation for the velocity of the dissolution waves at a distance r from the dislocation defect.

The evolution of the dissolution waves can be now quantified using equation (12) for v(r). Fig. 5a gives a typical velocity profile. The most salient features of the figure are the very fast velocities very close to the center of a defect and the existence of a minimum right after. The fast velocities near r = 0 are consistent with the postulation of a dislocation core, which had been proposed decades earlier just on the basis of thermodynamic energetic calculations (Cabrera & Levine, 1956). The dissolution wave will expand very quickly to the velocity minimum. The location of the velocity minimum will be labeled rpit. From rpit, the dissolution wave will then continue to expand into the rest of the crystal surface with a velocity that steadily increases until it reaches the limiting velocity vstep far from the dislocation center.

Using the equation for u(r), the equation for the minimum becomes
Of course, the solution to equation (13) also includes the position of the maximum in the velocity inside the dislocation core. Fig. 5b gives values of rpit for various surface free energies and for various values of rh. For large values of rh, there is no minimum, and rpit is not defined. Large values of rh reduce the strain energy [equation (13)] so much that no etch pits form on the crystal surface. On the other hand, for rh = 0, the size of rpit becomes exactly twice the classical hollow core radius (Cabrera & Levine, 1956; Cabrera et al., 1954).
Each dissolution stepwave begins at rpit and then proceeds to expand with velocity given by v(r) in equation (12). Each subsequent dissolution stepwave has to wait until the preceeding stepwave has moved a molecular distance, a, before it can begin to expand. Thus, the stepwaves are initially close together but eventually spread, reaching a limiting spacing given by (Lasaga & Lüttge, 2001)
The global dissolution rate is determined by the velocity and spacing of the dissolution waves. Therefore, the global dissolution rate is given by (Lasaga & Lüttge, 2001)
Equation (15) is the fundamental equation for the dissolution rate. Note, however, that for small ΔG values, v(rpit) approaches 0. The value of ΔG at which v(rpit) = 0 for the first time is a critical free energy, ΔGcrit. At this point, there is a problem with the expansion of the dissolution waves, and the rate goes essentially to zero (Lasaga & Lüttge, 2001) (i.e., it is limited by the usual nucleation or step dissolution mechanisms, which are much slower). This ΔGcrit will be important in quantifying the dissolution rate law of many different materials. The exact solution for ΔGcrit includes an rh dependence and from equation (12) is given by
with rpit determined by equation (13). Fig. 5c gives the exact solution of ΔGcrit as a function of rh for various surface free energies. It can be seen that the dependence of ΔGcrit on σ is significant, but the dependence on rh is not very strong.
With the theoretical developments of the previous section, we are now in a position to analyze the expected rate law for solid-fluid interactions during dissolution (as opposed to growth). First, the equation for the velocity of the step far from the etch pit, vstep, needs to be recast using the new results. We are interested here in the ΔG dependence of the overall rate, thus we need to extract the ΔG dependence of vstep. Coalescing the constants in equation (10), we can rewrite
\[v_{step}\ =\ {-}v_{{\infty}}(\frac{C_{B}}{C_{B}^{eq,{\infty}}}{-}1)Tanh(\frac{L}{2x_{s}})\]
where equation (10) has been divided by CBeq, and the various constants in equation (10) have been combined into the factor, v∞. Note that v∞ is defined as a negative number by convention. Using the thermodynamic relation for CB/Ceq,∞B (= exp(ΔG/kT)) and equation (14) for L (= Δr), we obtain
Combining equation (17) with equations (12) and (15), the overall rate law is given by
where the key function f(rpit) follows from equation (12):
Equation (18) predicts the overall dissolution rate for minerals based on the expansion of etch-pit induced stepwaves. In a sense, this equation is the counterpart to the BCF equation for dislocation-step-controlled crystal growth.

Fig. 6 compares our theory with results from lengthy Monte Carlo simulations with varying surface free energies, σ or bonding parameters, ϕ/kT. In the case of the MC simulations, all parameters needed for calculation of the rate in the dissolution rate theory are known. Hence, no fitting is done in the comparison. The excellent comparison between the dissolution theory and the MC simulations is a strong validation of the new stepwave theory developed here and in Lasaga & Lüttge (2001).

It is important to recast the rate law in a form that can be readily applied to experimental data in the future. Equation (18) can be rewritten as:
where Ro and B are constants, and where
A is determined from the value of the dissolution rate far from equilibrium (the so-called “dissolution plateau”). Therefore, only two parameters can determine the rate variation in the region up to the shutdown of the dissolution wave mechanism.
The final form of equation (20) is intriguing. Of major importance is that if ΔG becomes sufficiently negative (far from equilibrium), f(ΔG) approaches a constant value, and the rate law reduces to:
with constant R’o given by
\[R_{0}^{{^\prime}}\ {\equiv}\ R_{0}Tanh\left(\frac{B}{e^{{\Delta}G_{crit}/kT}}\right)e^{{\Delta}G_{crit}/kT}\]
Therefore, the overall rate equation acts just like the usual TST-like expressions far from equilibrium. Equation (22) provides a justification for the TST-like expressions under far-from-equilibrium conditions. In fact, experimental work carried out far from equilibrium would “extrapolate” nicely to zero rate at ΔG = 0 by using a TST-like linear approach. However, unless the kinetics were studied sufficiently close to equilibrium, the experimental data would miss the very significant nonlinear behavior of the rate. For example, Fig. 7 shows some of the work on silica dissolution (Gratz & Bird, 1993; Rimstidt & Barnes, 1980). Silica dissolution exhibits a well-defined linear dependence on deviation from equilibrium. However, the rate data are also perfectly consistent with the stepwave model except for the region near equilibrium, which is usually avoided. The region where the rate would drop significantly depends, of course, on the value of ΔGcrit. As Fig. 5 has shown, these values can be sizeable, especially for minerals which have a higher surface free energy than silica.

It is interesting that the BCF theory, although quantitatively of a different form from equation (20), also predicts a “linear” dependence of the growth rate on supersaturation far from equilibrium, which changes to a quadratic dependence on supersaturation as equilibrium is approached (Lasaga, 1998). Furthermore, the linear region does not have arbitrary intercept but extrapolates to zero rate at ΔG = 0.

Recent experimental data (Nagy et al., 1991; Nagy & Lasaga, 1992; Mogollon et al., 1996; Burch et al., 1993; Taylor et al., 2000; Cama et al., 2001) have indeed shown exactly the kind of strongly nonlinear response of the dissolution rate to undersaturation as predicted by the general theory (Lasaga & Lüttge, 2001). Note that it is impossible to “fit” equation (1) to these data. In these studies, the energetics of dislocation defects were correctly singled out as the source of the nonlinearities, but the general kinetic theory of the dissolution stepwaves had not been developed.

There are numerous implications of the general rate law to both natural and anthropogenic processes. Lasaga & Blum (1986) calculated the critical free energy. ΔGcrit for a variety of common minerals. These ΔGcrit values are rather large, i.e., from 0.5 kcal/mol to several kcal/mol (1 kcal = 4.18 kJ). However, most of the data on dissolution kinetics have been gathered at conditions very far from equilibrium. The result has been that in applications of these kinetic data to natural groundwaters, petrologic systems, radioactive waste disposal simulations, or geothermal systems, the data have been simply and linearly connected with the equilibrium point (i.e., Rate = 0 at ΔG = 0). As pointed out before (Nagy & Lasaga, 1992) and now generalized by the dissolution rate law, this treatment will significantly overestimate the size of the kinetic parameters as equilibrium conditions are approached.

Experimental data on several alumino-silicate minerals have provided verification of the overall rate law (Lasaga & Lüttge, 2001). It is interesting to consider the generality of the results discussed here. The clay mineral kaolinite, for example, has a rate that depends relatively simply on the saturation state (Fig. 8a). On the other hand, other layer minerals, e.g., gibbsite and smectite, do show a strong nonlinearity, as expected from the dissolution wave model (Fig. 8b and c). Layer silicates are an intriguing extension of the stepwave model, because the dissolution may occur along the edges of the layers rather than on the basal planes (Nagy, 1995). Although it is not clear whether the dominant dissolution is indeed perpendicular to the prism planes (i.e., dissolving along the sides of the sheets), it is an important possibility. One interesting question is whether the dissolution perpendicular to the prism planes can be controlled by the same defect mechanism as discussed here. Certainly, the effect of defects on a prismatic mineral surface in generating stepwaves will not be different than for any other mineral surface.

The implications of the overall stepwave mechanism for the treatment of fluid-rock reactions is also worthy of comment. All predictive models of reactive fluid flow must consider the possibility that the solutions reach near-equilibrium with at least some of the mineral phases. As a result, the models require a full treatment of the variation in rate as the saturation state varies. The only treatment available, and used of necessity by research workers, is the so-called transition state approach, which assumes a simple relation as given in equation (1). Our overall rate law, in fact, justifies at a deeper level equation (1), but only far from equilibrium. As equilibrium is approached and the rate decreases, the variation should be much more nonlinear. Some recent work (e.g., Gautier et al., 1994; Oelkers et al., 1994) has tried to ascribe all nonlinearity to “inhibition” effects from the same elements (e.g., aluminum or silicon) that are involved in defining ΔG. The problem is that, even if some complicated inhibition effect were to take place involving some particular element, the ΔG variation cannot be ignored, and none of these works have carefully measured the necessary ΔG variation as equilibrium is approached (rather the work is usually relatively far from equilibrium). In some cases, a lack of dependence of the rate on ΔG has been invoked. Yet that is precisely the problem, because a ΔG dependence has to be observed eventually, as demanded by thermodynamic requirements. Therefore, the possibility of “inhibition” does not probe behavior with respect to ΔG (especially based on the form of the rate law derived here). Rather than measuring the ΔG dependence, the full many-body (dissolution) process is treated as a simple elementary reaction, and the transition-state formula is once again invoked for the ΔG dependence. The formulation of the dissolution kinetics, which incorporates the crystal structure in contact with a solution into a full reaction mechanism such as done here, can account for both the ΔG dependence as well as any inhibition effects at the same time (Lasaga & Lüttge, 2003). The important non-linear ΔG behavior analyzed in this paper can then be contrasted with and placed in the context of these other factors, which play a role even further from equilibrium. The recent VSI observations and the systematic development of an overall dissolution rate theory show that only by incorporating the crystal into a treatment of water-rock kinetics can a consistent kinetic theory of use to society be available.

The MC method incorporates the elementary steps of bond formation and bond breaking in surface processes using microscopic reversibility. However, the method can follow the complex crystallographic changes in surface topography and bonding involved in any kind of crystal growth or dissolution. In particular, we can analyze the details of stepwave formation in the presence of defects and etch pits using such techniques.

The Monte Carlo method is based on molecular-scale kinetic and structural parameters, in addition to thermodynamic variables. An important parameter is the bonding strength between two neighboring atomic units, which is labeled ϕ.ϕ really represents the relevant energy change up on breaking the bond (e.g., ΔE for the reaction Si-O-Si + H2O = 2SiOH). The activation energy for dissolution is directly related to ϕ and for activation energies in the 10 to 20 Kcal/mol range, ϕ/kT∼4 to 8. The rate of detaching an individual surface unit with n bonds will be given by the rate constant,
Note that the time scale for the atomic detachment process will be determined by the size of the pre-exponential v, i.e., time is measured in I/v units, v also enters in the unit for dissolution rate of the MC results; that is the dissolution rate is given as atomic units/site/(i/v), which corresponds to the usual mol/cm2/sec used by experimental studies. Another atomic-scale parameter used in the MC calculations is the interatomic distance, a, between neighboring atoms. There are several other quantities that are obtained using v, a, and ϕ. If surface diffusion is included, the surface diffusion coefficient, Ds, is given by
where Ead is the activation energy for diffusion (Lasaga & Lüttge, 2003). The mean surface diffusion distance is then given by
where distance now is given in a units. Another important quantity derived from these parameters, is the surface free energy, o.
ϕ/kT can be related to the surface free energy of the crystal in the presence of an aqueous medium:
where a is the size of the molecular “unit” used in the MC simulation. For example, if a = 2 Å, T= 25°C, ϕ/kT= 6 corresponds to a σ of 308 ergs/cm2, a typical order of magnitude for mineral surface free energy (see Lasaga & Lüttge, 2003 for discussion).
Finally, the attachment of units from solution depends on both the bonding to the surface (i.e., ϕ) and the saturation state of the solution. If Δμ, is the difference in chemical potential of a species in solution and in the solid, then the attachment rate of i is given by
Δμ, is such that Δμi, = 0 at equilibrium; Δμ, is positive for supersaturated solutions and negative for undersaturated solutions. nkink is the number of bonds for kinks involving species i in the solid (for more details, see Blum & Lasaga, 1987; Lasaga & Lüttge, 2001, 2003). Note that in our notation Δμ/kT is the same as ΔG/RT, where ΔG is the free energy change of the dissolution reaction.

In summary, the Monte Carlo method incorporates information at the atomic level, e.g., the number of bonds, the bond strengths, the surface mobility of various species, the concentration of species near the surface, and carries out the kinetic evolution of a huge number of these particles, thereby producing the measurable macroscopic properties such as the surface structure, non-stoichiometric variations, and, of course, the overall dissolution rate.

The authors wish to thank R.S. Arvidson, A. Seilacher for thoughtful discussion. The authors also acknowledge gratefully support for this study from the US Department of Energy (grant # DE-FG02–90ER14153 to A.C. Lasaga and grant#’s DE-FG03–02ER63427 and DE-FG07–01R63295 to A. Lüttge), the National Science Foundation (grant # EAR-9628238 and EAR-9526794 to A.C. Lasaga and grant#’s EAR-0125667 to A. Lüttge), and the Alexander von Humboldt Foundation (to A. Lüttge).