## Abstract

A formulation based on defect-generated dissolution stepwaves of the variation of dissolution rate with the degree of undersaturation is validated by near-atomic-scale observations of surfaces, Monte Carlo simulations, and experimental bulk dissolution rates. The dissolution stepwaves emanating from etch pits provide a train of steps similar to those of a spiral but with different behavior. Their role in accounting for the bulk dissolution rate of crystals provides a conceptual framework for mineral dissolution far from equilibrium. Furthermore, the law extends research to conditions closer to equilibrium and predicts a nonlinear decrease in the rate of dissolution as equilibrium is approached, which has implications for understanding artificial and natural processes involving solid-fluid reactions.

The dissolution of minerals is central to a wide range of engineering, pollution, corrosion, and Earth science phenomena. Whether it be the weathering of rocks to form soils and control CO_{2} in the atmosphere, the possible reaction of radioactive containment materials, the breakdown of cement in buildings and dams, the leaching of aluminum or heavy metals in drinking water, or the flow of fluids deep in Earth's crust, no quantitative treatment is possible without an understanding of the kinetics of fluids reacting with solid surfaces. The lack of a general dissolution model means that the overall crystal dissolution rate and its variation with environmental conditions cannot be fully determined.

Historically, the science of crystal growth has advanced much more than that of crystal dissolution. The theory of surface-controlled crystal growth is dominated by an emphasis on the movement of steps and by two conceptual mechanisms: the nucleation of stable clusters and the spiral-generating Burton-Cabrera-Frank (BCF) theory (1, 2). Both mechanisms are related to the generation of growth steps that can then advance, resulting in overall crystal growth. With the use of atomic force microscopy and related technologies, many studies have described the presence of steps in a wide range of materials and monitored the motion of these steps during dissolution and growth (3–7). However, the scale limitations of these techniques have not yet allowed a quantitative integration of the movement of steps to produce a model for the overall dissolution rate.

For most practical crystal growth studies, the system is close to equilibrium. Therefore, the growth rate law plays a pivotal role in predicting and controlling crystal size, crystal shape, impurity content, or defect concentration in the product. In crystal dissolution studies, the system is often far from equilibrium conditions. Therefore, the treatments of crystal dissolution have taken three different approaches: (i) treat the crystal dissolution rate as independent of undersaturation (the distilled water case), (ii) treat systems as reaching near-equilibrium conditions quickly (assumed in many natural studies), or (iii) treat the kinetic rate law as a linear relation between rate and deviation from equilibrium. In the first case, the research has then limited the scope of investigation to the variation of the dissolution rate with variables such as ionic strength, Al^{3+} content, pH, and temperature (8). Consequently, these studies have ignored the ultimate approach of artificial or natural systems to equilibrium. In the third case, the most often invoked relation has been based on the principle of detailed balancing or transition state theory and leads to the rate law (2, 9)(1)where *A* and α are general constants, which could vary with pH, temperature *T*, or inhibitor molecules. Δ*G* (<0) is the Gibbs free energy of the dissolution reaction, and *R* is the gas constant. This equation most closely applies to the dissolution of rough surfaces and conditions that are far from equilibrium.

The formation of etch pits during crystal dissolution has been noted for decades (10, 11), and recent papers have described their formation at the nanometer scale (3, 12–19). Our model has arisen, in part, from the availability of surface topography data using vertical scanning interferometry [MicroXAM (ADE Phase Shift, Tucson, Arizona)] (20, 21). Although these studies have quantified the ubiquitous etch pits, the introduction of an absolute height reference for the surfaces has shown a far more important phenomenon (Fig. 1). During reaction, etch pits do form; however, the reference height indicates that the entire crystal surface is retreating in a global dissolution. This retreat has agreed closely with measured bulk dissolution rates for minerals such as anorthite (20). Therefore, the localized “growth” of etch pits is masking the widespread dissolution of the entire surface. This dissolution cannot be attributed to any spiral dissolution or to a coalescence of etch pits. Furthermore, as shown with dolomite (21), the global dissolution is fairly uniform and systematic over the entire crystal surface. Therefore, a surface-wide phenomenon must govern the bulk dissolution rate.

Burton, Cabrera, and Frank worked out the details of dislocation-driven growth in crystals. Cabrera and Levine (22) applied the spiral-growth theory to dissolution. However, if the undersaturation exceeded some critical value (Δ*G*
_{crit}), they showed that the strain field of a dislocation would open up the hollow core to form an etch pit (23, 24). This concept was also developed by Lasaga (25) and Lasaga and Blum (26) in their treatment of mineral dissolution. Experimental work on etch pits in quartz (27) agreed with theory. Therefore, beyond the critical undersaturation, spiral dissolution was not applicable.

If the undersaturation is low (i.e., |Δ*G*| < |Δ*G*
_{crit}|), etch pits will not open up and defect sites will not lead to stepwaves forming far from the localized etch pit. The result is a slow dissolution rate (e.g., spiral dissolution). For larger undersaturations, etch pits not only open up, but they are the source of a train of steps. These steps moving away from the defect into the surface lead to dissolution stepwaves, which can travel throughout the mineral surface and eventually control the bulk dissolution rate. This mechanism creates a train of steps in a manner reminiscent of BCF theory. However, the origin, spacing, and movement of the steps are different from the situation of a developing growth or dissolution spiral.

We analyzed the details of dissolution stepwaves in the presence of etch pits using a Monte Carlo (28–30) simulation (Fig. 2A). Nearly perfect crystals can dissolve by nucleating holes; however, this mechanism is slow at conditions close to equilibrium because the critical hole size becomes too large. An etch pit–forming dislocation defect is an efficient nucleating agent for holes and becomes a constant source of stepwaves, even for conditions not far from equilibrium. In a sense, the overall dissolution stems from waves much like those from a rock thrown into a still pond. Each etch pit would, in fact, represent a location where a rock has hit the water surface.

Another method of analyzing dissolution waves is to follow the expansion of the waves themselves by using contouring programs. Figure 2B shows an example of the close spacing of the contours close to the etch pit and the expansion of the waves as they move away from the center. This process continues inexorably and produces dissolution stepwaves that can travel across the crystal surface.

Natural crystals contain numerous dislocations, and the train of dissolution stepwaves from each of the etch pits will combine. However, unlike the waves in a pond, the dissolution stepwaves will annihilate each other. The overall conceptual model then is as shown in Fig. 2C. If one were to analyze the crystal surface at a later time, the etch pits would be visible surrounded by fairly flat regions of the surface (e.g., Fig. 1). However, throughout the dissolution, stepwaves have been emanating from each of the etch pits and combining to lower the entire mineral surface. Therefore, the actual topography would look like that in Fig. 2C.

A consequence of the dissolution stepwave model is that the bulk dissolution rate will depend only slightly on the etch pit density (Fig. 2C). The overall rate for various etch pit densities converges to the single etch pit rate, except for an initial short transient time. This prediction is in accord with experimental investigations (31–34) and resolves a long-standing problem in water-rock kinetics.

The dissolution waves emanating from an etch pit will be closely approximated by circular stepwaves (Fig. 2B). In 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, the Gibbs-Thomson equation needs to be altered to include the strain field of dislocation defects, *u*(*r*)(2)where μ is the bulk shear modulus, *b*is the size of the Burgers vector, *K* (≈1) depends on the dislocation, and the radius *r*
_{h} limits the magnitude of the strain energy close to the dislocation center (35, 36). For a circular dissolution step with radius *r*, centered on a dislocation defect, the equilibrium solubility, modified by both surface area considerations as well as the strain energy, will be given by(3)where σ is the surface free energy,*v* is the molecular volume, and *k* is Boltzmann's constant (2). This modification affects the equilibrium concentration of dissolved molecules around the dissolution stepwave and is a function of *r*. On the basis of Eq. 3, the step velocity becomes(4)where *v*
_{step}is the velocity of a series of straight steps.*v*
_{step} can be obtained from treatments similar to the BCF theory [e.g., (2)]. Equation 4 gives the velocity of dissolution stepwaves at a distance *r* from the dislocation defect.

As a check on Eq. 4, the strain energy term can be set equal to 0. In this case, if we also assume that the σ and Δ*G* exponent terms are small, then(5)which is the usual result from standard crystal growth physics [*r*
_{crit} is a critical radius for the formation of a two-dimensional “hole,” e.g., (2)].

The evolution of dissolution stepwaves can be quantified with Eq. 4 for*v*(*r*). Typical velocity profiles show that the faster velocities are close to the defect center, and a minimum velocity occurs at a distance labeled as *r*
_{pit}(Fig. 3). The fast velocities near *r* = 0 are a dynamic assertion that there is a dislocation core. The position of the minimum is determined by(6)The dissolution wave will expand quickly to the velocity minimum at *r*
_{pit}, circumventing any nucleation problems. These stepwaves then emanate from the core at position *r*
_{pit} and continue to expand into the rest of the crystal surface with a velocity *v*(*r*) that steadily increases until it reaches the limiting velocity*v*
_{step} far from the dislocation center.

Once a dissolution stepwave has commenced, the next stepwave has to wait until the first wave has moved a molecular distance before it can begin. Therefore, if the first wave begins at time*t* = 0, the second wave will begin its expansion at time*t*
_{w} =*a*/*v*(*r*
_{pit}), where*a* is a molecular distance. As time proceeds, both stepwaves (and all later ones) continue to move out into the mineral surface with a velocity given by Eq. 4. On the basis of Eq. 4, for long times (large *r*), *v*(*r*) just becomes *v*
_{step}. Therefore, after the dissolution stepwaves have traveled sufficiently far from the etch pit, the separation of the waves Δ*r* is given by Δ*r* = *a*[*v*
_{step}/*v*(*r*
_{pit})]. The spacing is initially small and increases to a steady-state value as seen in the Monte Carlo results. Therefore, the global dissolution rate is given by(7)The velocity of the dissolution waves reaches 0 once the undersaturation drops below Δ*G*
_{crit} (Fig. 3), where(8)At this point, there is a problem with the expansion of the dissolution waves, and the stepwave mechanism shuts down. Once the undersaturation is reduced below the critical undersaturation, the rate does not necessarily go to 0; rather it will be limited by the much slower nucleation or spiral dissolution mechanisms. Therefore, the rate will decrease substantially.

With the theoretical developments of the previous section, we can analyze the expected rate law for solid-fluid interactions during dissolution. On the basis of the usual treatments of*v*
_{step} (1, 2), we can rewrite the overall dissolution rate law as(9)where the modifying function*f*(*r*
_{pit}) follows from Eq. 4
(10)and where *x*
_{s} is a surface diffusion distance expressed in molecular units. Equation 9 is a general rate equation for dissolution and is the counterpart to the BCF equation for step-controlled crystal growth.

One of the most prominent features of the rate law based on the dissolution stepwaves is the nonlinear variation in the rate as equilibrium is approached. This nonlinearity is in contrast to the usual assumptions made by many dissolution studies. BCF theory also predicts nonlinear behavior as near-equilibrium conditions are reached. However, unlike BCF theory, the dissolution stepwave theory predicts a substantial reduction in the rate once Δ*G* ∼ Δ*G*
_{crit}, which for many crystals occurs far from equilibrium.

Equation 9 can be tested with detailed kinetic models such as the Monte Carlo method. In Monte Carlo simulations, σ,*u*(*r*), and *x*_{s} are known exactly. The only adjustable parameter in the comparison, therefore, is the dissolution rate at infinite dilution [the so-called “dissolution plateau” (37)]. The shape of the dissolution rate law is otherwise completely determined by the parameters in Eq. 9. The agreement of Monte Carlo simulations and the general rate law in Fig. 4 validates the dissolution theory.

It is important to recast the rate law in a form that is readily applied to experimental data. The equation can be rewritten as(11)where *A* and *B* are constants and where(12)This equation has three parameters: *A*,*B*, and Δ*G*
_{crit}. *A* is determined from the value of the dissolution rate far from equilibrium (the dissolution plateau). Therefore, only two parameters, *B*and Δ*G*
_{crit}, determine the rate variation in the region up to the shutdown of the dissolution wave mechanism.

Recent experiments (38–42) have shown a nonlinear response of the dissolution rate to undersaturation as predicted by the general dissolution stepwave theory. 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 was not developed.

The quantitative treatment of the various nonlinear experimental results can now be extracted with the dissolution rate law based on the dissolution stepwaves. In applying Eq. 11 to the experimental data, values of Δ*G*
_{crit} can be calculated from known mineral properties (35), leaving just*A* and *B* to be determined from dissolution experiments. The general dissolution rate law explains the dissolution behavior of albite, gibbsite, labradorite, and smectite (Fig. 5).

The general dissolution rate law predicts a different behavior of reacting solids near equilibrium and far from equilibrium. In fact, for conditions far from equilibrium, Eq. 11 approaches Eq. 1, validating its use far from equilibrium. The strong nonlinearities of the rate law, however, can lead to interesting pattern formation, oscillatory behavior, and metastable reactions (43).

This difference in behavior as equilibrium is approached has important implications in many areas of Earth sciences, where equilibrium is often assumed. Values of Δ*G*
_{crit}for many minerals are rather large (i.e., from 0.5 to several kcal/mol) (26). Mineral dissolution rates near equilibrium have been applied to natural groundwaters, petrologic systems, radioactive waste disposal simulations, or geothermal systems, where the data have been linearly connected with the equilibrium point (i.e., rate = 0 at Δ*G* = 0). This treatment will substantially overestimate the size of the dissolution rates as equilibrium conditions are approached (38). The error can be more than an order of magnitude and more important, therefore, than variations in pH or temperature in many cases. The new theory, in fact, helps explain the long-standing discrepancy between measured laboratory rates and slower estimates of rates from field observations (44).

If the Δ*G* dependence of mineral reactions is extended to that of overall whole-rock reactions such as those used in petrologic analyses, the nonlinearity predicts sizable deviations from equilibrium. Even a deviation of 0.5 kcal/mol for an overall reaction with an entropy change of Δ*S* = 20 cal/mol per K would yield a temperature overstep of Δ*T* = 25°C. This result suggests that many metastable reactions will play a role in the changing mineral assemblages in nature. This result is in strong contradistinction with the often-quoted deviations of 1° or 2°C in the literature (45). In addition, any measurements of rates in the field will yield rates that are orders of magnitude slower than laboratory rates measured far from equilibrium (46).

↵* To whom correspondence should be addressed. E-mail: aluttge{at}rice.edu