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.
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.
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).
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.
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 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.
Full rate law
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 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.
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).