Transport in porous media is quite complex, and still yields occasional surprises. In geological porous media, the rate at which chemical reactions (e.g., weathering and dissolution) occur is found to diminish by orders of magnitude with increasing time or distance. The temporal rates of laboratory experiments and field observations differ, and extrapolating from laboratory experiments (in months) to field rates (in millions of years) can lead to order-of-magnitude errors. The reactions are transport-limited, but characterizing them using standard solute transport expressions can yield results in agreement with experiment only if spurious assumptions and parameters are introduced. We previously developed a theory of non-reactive solute transport based on applying critical path analysis to the cluster statistics of percolation. The fractal structure of the clusters can be used to generate solute distributions in both time and space. Solute velocities calculated from the temporal evolution of that distribution have the same time dependence as reaction-rate scaling in a wide range of field studies and laboratory experiments, covering some 10 decades in time. The present theory thus both explains a wide range of experiments, and also predicts changes in the scaling behavior in individual systems with increasing time and/or length scales. No other theory captures these variations in scaling by invoking a single physical mechanism. Because the successfully predicted chemical reactions include known results for silicate weathering rates, our theory provides a framework for understanding changes in the global carbon cycle, including its effects on extinctions, climate change, soil production, and denudation rates. It further provides a basis for understanding the fundamental time scales of hydrology and shallow geochemistry, as well as the basis of industrial agriculture.

Chemical reactions at the surfaces of porous media drive the weathering of silicate minerals in the earth's crust. This weathering affects surface denudation rates and the global carbon cycle, and thereby also climate change and extinctions through geologic time. But reactions within porous media exhibit non-trivial dependences on time, space, and measurement scales that cannot be interpreted as transport-limited in the context of Gaussian transport. Thus, the role of solute transport in controlling such reaction rates has been controversial. Because our transport theory, which is based on concepts from percolation theory, is the first theory that can actually predict non-Gaussian solute transport in porous media, we are able to resolve the controversial aspects of chemical reaction rate scaling. In particular, we show that our calculated functional dependence of the solute velocity acts as a proxy for chemical reaction rates, proving that they are limited by transport. Importantly, the solute velocity is not based on diffusion; rather, it is generated by the flow of fluids in the medium, and is nearly a power law in both space and time. A wide range of experimental results and field observations are united by this theory, including experiments on the dissolution of radioactive elements, chemical weathering rates, weathering rind formation on surface clasts, and seafloor weathering rates. The agreement between theory and experiment extends over 10 orders of magnitude of time scales.

## INTRODUCTION

Predicting the scaling of chemical reaction rates in porous media is important but difficult. In this study, we focus on the temporal and spatial dependences of surface reactions, rather than on their initial values. The relevant physics addresses a complex of interacting processes, scales, and materials, and researchers must ultimately choose between differential equations for transport in homogeneous media (Bear, 1972), and difference equations for transport in heterogeneous media (Berkowitz and Scher, 1995; Berkowitz *et al.*, 2006; and Hunt *et al.*, 2011). The issue has practical significance, as it concerns (for example) weathering of silicate minerals, asserted (Algeo and Scheckler, 1998; Berner, 1992; Dixon *et al.*, 2009; Raymo, 1994; Sheldon, 2006; Vance *et al.*, 2009; and Anderson and Anderson, 2010) to be a driver for variations in the global carbon cycle, with impacts on episodes of extinction, climate change, soil production, and surface denudation rates. Yet reaction rate scaling is even more broadly important, informing the transport and fate of subsurface radioactive contaminants, the time scale of soil formation (Loague and Corwin, 2006), and even the role of fertilization in agricultural time scales.

Many factors influence rates of chemical reactions in natural porous media (Brantley *et al.*, 1986; Dentz *et al.*, 2011; Li *et al.*, 2008; Noniel *et al.*, 2012; Navarre-Sitchler *et al.*, 2007; and Raoof and Hassanizadeh, 2010), including (1) concentration gradients, (2) solute advection or diffusion, and (3) physical and chemical changes in mineral surfaces. If we are considering natural processes relevant to chemical weathering, then we must also include changes in (4) flow rates, 5) surface conditions, (6) climate, and even (7) tectonics. We postulate that, for chemical reactions that are transport limited, the reaction rate is proportional to the solute velocity. In that case, reaction rates can be predicted from the solute transport velocity and its scaling with time or distance. As we show here, the velocity scaling can be given by conservative solute transport theory (Hunt *et al.*, 2014). It is critical to note that at scales greater than the pore scale, our theoretical construction does not yield equal solute and fluid velocities, because the solute travel time is proportional to the solute travel distance to a power greater than 1 (Lee *et al.*, 1999; Sheppard *et al.*, 1999). Such a result is also found in theoretical descriptions employing the continuous time random walk (CTRW) (Scher and Montroll, 1975; Scher *et al.*, 1991; Berkowitz and Scher, 1995; and Berkowitz *et al.*, 2006), or fractional advection-dispersion equations (FADE) (Pachepsky *et al.*, 2000; Meerscheart *et al.*, 1999, 2002; Sanchez *et al.*, 2005; and Park *et al.*, 2005). In percolation theory, CTRW and the FADE, the same non-linear scaling relationship between transport distance and transport time is obtained; what distinguishes percolation theory, though, is that it predicts the precise value of the power (Lee *et al.*, 1999; Sheppard *et al.*, 1999). Our interpretation of the distinction between solute and fluid velocities, which arises from the power-law scaling, is that it is a product of the fact that the fluid flow is incompressible, while the solute transport is not. In fact, if the shape of the solute distribution does not change with transport distance (Scher *et al.*, 1991), then the mean solute density must be inversely proportional to the mean solute transport distance.

Our treatment explains nearly all observations through a single functional dependence, uniting data from many different experiments and field studies in a single curve. The predicted functional dependence in time, *t*, for the solute velocity (and thus the reaction rate) is nearly a power-law, while the spatial dependence (on *x*) deviates further from power-law behavior. Functional dependences are referenced to a fundamental (pore-scale) fluid velocity *v*_{0} (which can often be determined or guessed) and to either a fundamental time scale *t*_{0} or a fundamental spatial scale *x*_{0}. When spatial heterogeneity is relevant at the scale of a single pore, then *x*_{0} is the pore length, the pore crossing time *t*_{0} = *x*_{0}/*v*_{0}, and from any two of these values (*x*_{0}, *t*_{0}, and *v*_{0}), the third can be determined. The assumed proportionality between the reaction rate and the solute velocity generates for reaction rates the same dependence on *x* and *x*_{0}, or *t* and *t*_{0} as for solute velocities. The interpretation of the scaling of surface reactions in porous media is thus greatly simplified.

Our conclusions do not depend on specifics such as mineralogy or tectonic regime. We have seen relevant variations in only two parameters, namely, temperature and velocity of the advecting fluid, both of which are affected by climate. In one experimental case (Liu *et al.*, 2008), the advective flow paths appeared to exhibit two-dimensional connectivity in invasion percolation rather than the expected three-dimensional connectivity in random percolation. The observed change in scaling was exactly as expected from the fractal dimension of the percolation backbone, *D*_{b} (Sheppard *et al.*, 1999), in different universality classes. The choice of universality class depends on conditions of saturation (the fraction of the pore space occupied by water) and the connectivity of the flow paths. Two-dimensional invasion percolation is consistent with unsaturated flow (the pore space only partially filled with water) along surfaces (e.g., fractures; Glass *et al.*, 2005), while 3D random percolation applies to saturated flow (the pore space completely filled with water). Weathering rind formation on surface clasts appears to follow scaling predicted for unsaturated flow conditions, as would be expected.

## THEORY

Our theory for conservative solute transport was first developed in Hunt and Skinner (2008), while the temporal evolution of the solute's spatial distribution was first discussed in Hunt and Skinner (2010a) and Hunt *et al.* (2014). Below we give a brief overview. New in this communication are the explicit calculation of the solute velocity from the temporal evolution of the spatial solute distribution, and the comparison of predicted solute velocity with experimental results from reaction rate scaling.

Many features of solute transport are poorly described by the traditional continuum mechanics approach. In particular, the Advection-Dispersion Equation (ADE), based on the assumption of Gaussian mixing (Bear, 1972), does not predict the observed long tailed distributions (Cortis and Berkowitz, 2004) nor does it yield the observed proportionality between the variance in the distribution and the square of the solute transport distance (Berkowitz and Scher, 1995; Berkowitz *et al.*, 2006). But the transport theory of Hunt and Skinner (2008), based on percolation theory applied to groundwater transport, both reproduces and predicts the experimental results. In particular, we see (Ghanbarian-Alavijeh *et al.*, 2012) that the fractal dimension of the percolation backbone, together with system-specific information about the pore size distribution, accurately predict the saturation dependence of the solute arrival time distribution. The predicted and experimentally observed distributions agreed almost exactly over three orders of magnitude, with only slight adjustment of the fundamental time scale, *t*_{0} (Ghanbarian-Alavijeh *et al.*, 2012).

Consider the flow paths through a complex, disordered medium described by a wide distribution of local flow resistances. Under a wide range of conditions, critical path analysis (Pollak, 1972; Hunt, 2001) can generate accurate results for the hydraulic conductivity *K* (Ghanbarian-Alavijeh and Hunt, 2012). The basis of critical path analysis is that *K* scales inversely with the smallest possible value of the largest resistance on a continuously interconnected path through the system. This conceptual approach addresses the two acknowledged barriers (Freeze, 1975; Gelhar, 1986; and Dagan and Neuman, 1997) to understanding flow and transport: heterogeneity and connectivity.

Our theoretical description of solute transport, while related to this picture of *K*, is based on enumerating all the transport paths through a system of a given length. We organize the enumeration using the structure of percolation theory (Stauffer and Aharony, 1994; Sahimi, 1994; Hunt, 2001; and Hunt and Skinner, 2010a). The organization has its conceptual basis in the tendency for flow to maximize dissipation by taking paths of lowest cumulative resistance. But avoiding large resistances tends to increase the length of a path. Critical path analysis uses the percolation threshold to define the smallest resistance value for which a continuous interconnected path across the system can be found. Exactly at the threshold, however, the length of interconnected paths tends to infinity. This competition, between minimizing the controlling resistance value, and minimizing the path length, allows determination of an optimal theoretical arrival time, which corresponds to the peak arrival time in experiments. The long tail in the distribution results from the increasingly tortuous pathways, with tortuosity diverging precisely at the threshold (Hunt and Skinner, 2008). An important distinction (Lee *et al.*, 1999) is that the fractal dimension of the percolation backbone governs solute arrival times, whereas the fractal dimension of the optimal (shortest; Cieplak *et al.*, 1996; Porto *et al.*, 1998) paths determines the geometrical tortuosity, relevant for path length calculations, and also useful for hydraulic and electrical tortuosity. The time result has been verified by comparison with transient photoelectric currents in amorphous semiconductors and polymers (Hunt *et al.*, 2011), while the distance has been verified through comparisons with measured and calculated hydraulic, electrical, and diffusive tortuosities (Ghanbarian *et al.*, 2013a; 2013b).

Consider first the probability distribution function, *W*(*g*), of local conductance values, *g*. These local conductances correspond to regions of pore space in the model, so they also contribute to the total pore volume. For a regular network model, *W*(*g*) can be normalized as

We are interested in the integral

which represents the cumulative distribution of conductance values greater than a given value, *g*. As *g* is reduced from its largest possible value, the probability *p* is increased. The value of *g* for which *p* = *p*_{c}—the critical percolation probability for the given network—is defined as *g*_{c}. The procedure to generate *g*_{c} thus satisfies our requirement that it represent the largest possible value of the smallest conductance on a path of continuously interconnected conductances that spans the system.

This problem was reformulated in terms of the total pore volume by allowing *g* to be a specific function of pore size. Many complicated methods for making this association are available in the literature, but we followed a relatively simple procedure. Initially, we used the Rieu and Sposito (1991; henceforth RS) discrete random truncated fractal model as a parsimonious means to represent pore-size heterogeneity, in other words for expedience. In order to be able to account for a wider range of media, we later generalized this model to a continuous form of the Bird *et al.* (2000) pore-solid fractal mode (e.g., in Ghanbarian-Alavijeh *et al.*, 2012). So far, this theoretical construction has been almost universally successful, indicating that either (1) it is well suited to describing many real systems, or (2) details of the medium are not critical to the final result.

To date our comparisons with *chemical reaction rates* have only utilized the original, and simpler, RS model. In the present treatment, we write the porosity *ϕ* as a fractional pore volume, integrated over all pore radii. The volume of a pore is proportional to *r*^{3}, while the probability distribution function (pdf) for the pore radius *r* is given by *W*(*r*) (between the limits *r*_{0} and *r*_{m}). The pore volume with radius between *r* and *r* + *dr* is *r*^{3}*W*(*r*). Use of *r*^{3} is consistent with any pore shape which preserves the relationships between the mean pore dimensions with increasing pore size (a constant aspect ratio). The porosity *ϕ* is thus

For a truncated power law pore size distribution, *W*(*r*) is proportional to *r*^{−D−1}, where *D* is the fractal dimensionality of the pore space. For example, a value *D* = 1.0 indicates a narrow pore size distribution; wide pore size distributions have *D* approaching 3.0. A suitable choice of normalization constant now retrieves the known RS result for the porosity

While not shown here, the result for the RS soil water retention curve also follows from this treatment. Assuming a constant aspect ratio across all pore sizes, combined with Poiseuille's law for fluid flow through a tube, generates a conductance (see Hunt, 2001) that is proportional to *r*^{3} as well. The percolation probability *p* and its critical value *p*_{c} now correspond to a volume fraction, *V*, and its critical value, *V*_{c}. Thus, one can write for an arbitrary smallest pore radius, *r* ≈ *g*^{1/3}, on a path of interconnected conductances

With these equations, we now find an expression for *V* – *V*_{c}:

The probability of finding an interconnected path of conductances greater than or equal to this particular *g*, that spans a system of length *x*, is then found using the cluster statistics of percolation (Hunt, 1998). This probability can also be regarded as the probability that a given value of a bottleneck resistance governs the flow through the cluster. The cluster statistics of percolation theory describe the frequency of occurrence of clusters for any bond (site or volume) fraction *p* with a given number *s* of interconnected bonds or sites, where *s* is interpreted as a cluster volume (Stauffer, 1979). The process of generating the desired distribution involves integrating the cluster statistics over all cluster sizes larger than the system in question, and also transforming the variables appearing in the cluster statistics (Hunt, 1998). One such transformation, mentioned above, is the substitution *p* – *p*_{c} → *V* – *V*_{c}, which converts the system from site or bond percolation to continuum percolation. Following this transformation, Eq. (6) allows us to generate a function of the minimum conductance value, *g* and its critical value, *g*_{c}. The second significant change of variables arises from the need to change cluster volumes to cluster lengths; in principle, this involves the mass fractal dimensionality of the percolation clusters. However, use of the dimensionally dependent scaling law (Toulouse, 1974) generates a length dependence that is independent of details of the cluster (Hunt, 1998), and even (to lowest order approximation) independent of dimension, provided one uses the dimensionally correct value of the critical exponent *ν* for the correlation length.

For a given rate-limiting hydraulic conductance *g*, at a given distance *x* from the source, the desired pdf for the probability that a specific value of a bottleneck resistance governs the flow through the cluster is given in terms of an exponential integral, *Ei* (Hunt and Skinner, 2008)

where *L* is a fundamental length scale (which for weathering rates on particle surfaces we take to be on the order of a particle size), *g*_{c} is the critical hydraulic conductance (proportional to *r*_{c}^{3}, where *r*_{c} is the critical pore radius in the medium), *ν* is the correlation length exponent from percolation theory (0.88 in a three dimensional system), and *g*_{min} and *g*_{max} are (Ghanbarian-Alavijeh *et al.*, 2012)

Here, *θ*_{t} is the percolation threshold of the medium, expressed as a volumetric moisture content.

The topology of percolation theory, and the distribution of controlling conductances for *g*_{min} < *g* < *g*_{max}, are then used to calculate the time that solutes take to travel across the system on such a path. This time, *t*, is a function of *g*

where *t*_{0} is the pore crossing time, and

where *D*_{b} is the fractal dimensionality of the percolation backbone. Note that *h*(*g*) diverges as a power law for *g* = *g*_{c}. The exponents *D*_{b} and ν, whose values are known from percolation theory, are critical to the form of the long-tailed distribution and thus the scaling of the solute velocities. Equations (7) and (9) (Hunt and Skinner, 2008, and subsequent publications) may be combined using *W*(*t* | *x*)d*t*= *gW*(*g* | *x*)d*g* to find the arrival time distribution, *W*(*t* | *x*) as a function of position, *x*, where *W*(*g* | *x*) is calculated using Eq. (7), and d*t*/d*g* using Eqs. (9) and (10). Note the factor *g* in the identity. This factor makes the solute transport proportional to the fluid flux on each particular pathway, since that flux is proportional to the smallest conductance included. The actual calculation of the arrival time distribution involves a numerical procedure which was described in detail by Hunt and Skinner (2008).

Following Hunt and Skinner (2008) and Hunt *et al.* (2011), we used the results of Lee *et al.* (1999) to assign *D*_{b} as the exponent for the typical system crossing time. Thus, the typical time taken for a solute particle to enter a system at one side and exit the system a distance *x* downstream is proportional to *x ^{D}*

^{b}, where

*D*

_{b}> 1. In three dimensions and for random percolation (appropriate for fully saturated conditions; see Sheppard

*et al.*, 1999),

*D*

_{b}= 1.87. Since

*D*

_{b}> 1, the solute velocity is a decreasing function of both time and of solute transport distance. If

*D*

_{b}= 2, and

*h*(

*g*) could be ignored, then solute velocity would scale exactly like diffusion. But because these conditions are not fulfilled, we have the result that the solute velocity scaling looks very much like diffusion at first, but changes to a steeper dependence at later times. At the length scale of a single pore, where the ADE is valid (Neuman, 1990; Hunt

*et al.*, 2011),

*t*

_{0}is proportional to a distance divided by a fluid velocity: the pore-scale solute velocity is identical to the pore-scale fluid velocity. Since

*t*

_{0}is actually obtained as an integral over the conductance distribution,

*t*

_{0}is also inversely proportional to the minimum conductance on the path, meaning that the solute velocity so obtained is equivalent to a solute flux.

When the known statistics of *g* are applied, Eq. (9) might be thought useful for finding a mean travel time. However, with *D*_{b} = 1.87, the tail of the distribution is so long that such a mean travel time does not exist, as discussed in Berkowitz and Scher (1995). But even without a mean travel time, one can still find the mean position, 〈*x*〉, of the solute at any given time, from an integral over the derived spatial distribution of the solute, using an equation analogous to Eq. (9) (Hunt and Skinner, 2008). Such a calculation was verified (Hunt *et al.*, 2011) to give a result for the typical (not mean) time for a particle to cross the system, *t*(*x*) ∝ *x ^{D}*

^{b}, in the case of electron transport in disordered semiconductors and polymers (Pfister, 1974; Pfister and Scher, 1977; Pfister and Griffiths, 1978; Tiedje, 1984; and Bos

*et al.*, 1989). Particularly important in this comparison is that

*D*

_{b}values from random percolation (Sheppard

*et al.*, 1999) generate the appropriate scaling of typical travel time as a function of system length, for both the three-dimensional transport in amorphous semiconductors and the two-dimensional transport in polymer systems.

Using 〈*x*(*t*)〉, one can obtain a mean solute velocity: *v =* d *⟨ x*(*t*) *⟩ /* d*t*. When the value of the solute distribution variance is also known, one can predict related quantities, such as the longitudinal dispersion coefficient and the dispersivity as functions of either transport distance or of transport time (the elapsed time over which transport occurs) (Hunt and Skinner, 2010a; 2010b; and Hunt *et al.*, 2011). The predictions given in those papers were made with the same calculations that are applied here. The mean solute velocity predicted depends on transport time and can be written as

where *v*_{0} ≡ *L*/*t*_{0} and the function *f* must be determined numerically via the procedure referenced above using the probabilistic transformation. However, evaluation of Eq. (11) at *t* = *t*_{0} gives *f*(1) = 1. Using the mean solute velocity, one can transform Eq. (11) so that *x* is the independent variable. Note that the approximate equality, which is all that can be obtained analytically, holds only for certain ranges of time, specifically at smaller time values. The full numerical solution predicts a change to a steeper reduction in reaction rate with time scale at larger length scales. The hypothesis here, that surface reaction rates *R*(*t*) are solute transport limited, makes them proportional to solute velocities and provides an unexpected opportunity to verify that the theoretical result in Eq. (11) holds also for groundwater transport. Therefore

where *R*_{0} is the initial reaction rate, i.e., the value of *R* when *t* = *t*_{0}, or *x* = *x*_{0}. Again, the approximate equality holds only for early time. We stress that we do not have a general theoretical approach for finding *R*_{0}, the value of which will depend on quantities such as fluid velocity, saturation conditions, temperature, pH, and potentially other inputs such as changes of the surface characteristics over time. But what has been considered most puzzling about reaction rates has been their time (or length) dependence, and that is what is addressed in Eq. (12). Equation (12) can be transformed so that *x* is the independent variable by substituting *t*(*x*), as obtained from Eq. (9). After such a transformation, the approximate power-law result for the reaction rate would be *R*_{0}(*x*/*L*)^{1−Db}, which can also be written as

When Eq. (13) is valid (at shorter length scales, as will be seen), reaction rates decay according to the power 1 − *D*_{b}.

The value of *t*_{0} provides a relevant scale to the time (horizontal) axis. Particle-sized reactant sources (e.g., for chemical weathering and Uranium dissolution) provide a particle-size scale for heterogeneity, meaning that *L* can often be assumed to be a characteristic particle size, as noted above. We typically do not know the value of *t*_{0} in Eqs. (11) and (12), although it can be calculated under well-defined laboratory conditions. For comparison with experimental data, it is also necessary to define a vertical scale. This can be either a fundamental rate constant (*R*_{0} from Eq. (12)), or a fundamental velocity (*v*_{0} from Eq. (11)).

Sometimes one is interested in how reaction rates depend on the scale of measurement, rather than the scale of transport. For example, Navarre-Sitchler and Brantley (2007) found a power-law increase of apparent reaction rates with increasing scale of measurement, and interpreted this result in terms of a reaction front with a fractal structure. Percolation theoretical treatments also generate such fractal surfaces, and in principle, it is possible to calculate the geometry of such a reaction front within the same theoretical framework that generates the solute velocities. In particular, the perimeter of a percolation cluster in three dimensions has *two* contributions (Kunz and Souillard, 1976; Hunt and Ewing, 2009): one proportional to the square of the radius, and one proportional to the volume of the cluster. The theoretical results from percolation theory apply only when the linear dimension of such a cluster is on the order of at least ten individual units (such as bond lengths, pore separations, or a surface roughness scale; Hunt, 2001). Let us define the radius of a large cluster as the correlation length, *χ*, from percolation theory. Consider first the term in the surface area proportional to the cluster volume. The volume of a three-dimensional percolation cluster of length *χ* is proportional to *χ* ^{2.5}: *M* $\u221d$ *χ ^{D}*

^{m}, where

*D*

_{m}≈ 2.5 is the universal mass fractal dimension of large clusters near the percolation threshold (Stauffer and Aharony, 1994; Hunt and Ewing, 2009). This means that the surface area

*A*has two terms: one proportional to

*χ*

^{2}, and the other proportional to

*χ*

^{2.5}. Thus, we have

Equation (14) contains an unknown numerical factor *C*, and a scale factor *χ*_{0}. The scale factor *χ*_{0} represents a length scale above which the fractal properties of the surface of percolation cluster begin to dominate, so *χ*_{0}^{2} is also a fundamental area. The value of *χ*_{0} cannot be smaller than the size of a molecule that probes the surface. In fact, however, the above argument on the minimum cluster size implies the necessity of using a value for *χ*_{0} of at least a factor 10 larger than a molecular diameter, since the cluster surface does not exhibit fractal characteristics at smaller length scales.

## COMMENT ON SCALE FACTORS

It has already been pointed out (Chao *et al.*, 2000) that two values of the length scale can be important for solute transport. One is the transverse extent of the initial solute distribution, and the second relates to the medium. If there is only a single scale of heterogeneity in the medium, then the larger length scale, whether characteristic of the medium or the plume, is most important. However, when the medium is heterogeneous at many (or all) length scales, then the only relevant length scale is that set by the experimenter. In a previous work (Hunt *et al.*, 2011), we suggested that the uniformity in the observed data for the dispersivity (1600 values), which were all consistent with the choice of a length scale of 1 m, was due to the relative consistency of the size of experimental apparatus, which in turn was related to the size of the experiment. We will see that for most experiments related to chemical reaction rates in porous media, the relevant length scale is close to the size of a single pore, meaning that the relevant heterogeneity scale is set by the variability in particle mineralogy. Interestingly, a quite different dominant length scale emerges when we apply our calculations to industrial agriculture.

## NATURAL FLOW RATES

In order to compare the theory with field weathering rates and fundamental hydrological scales, it is necessary to characterize flow rates. We contend that for the majority of field and laboratory data related to these topics, reaction rates are controlled primarily by transport. The following data summaries indicate that the predominant range of observed flow velocities is almost precisely the theoretical range of fluid flow velocities over which reaction rates are transport limited.

Consider a saturated hydraulic conductivity value *K*_{S} ≈ 10^{−4} cm/s (a median order of magnitude for crustal materials; Anderson and Anderson, 2010) under conditions of gravity flow. This combination, for porosity *ϕ* = 0.5, would be consistent with a single pore flow rate of about 2 *μ*m/s. This flow rate, i.e., 2 × 10^{−6} m/s, is the median value given by Blöschl and Sivapalan (1995) as dominating groundwater flows (unconsolidated media). The darkly shaded region in Fig. 2 of Blöschl and Sivapalan (1995), within which most measurements are alleged to fall, corresponds to a variation of about two orders of magnitude in flow velocity centered around this value. While *K*_{S} in crystalline rocks can be much smaller, similar *K* values are found in sedimentary rocks where, however, a smaller porosity generates larger pore-scale flow rates. Such a range of flow rates is nearly the same as the 0.5 *μ*m/s suggested by Maher (2010), and the 100 *μ*m/s suggested by Molins *et al.* (2012), as bounding the regime where chemical reactions are transport-limited. In particular, reaction rates are shown by Maher (2010; her Figures 4 and 5) to be proportional to flow rates in the transport-limited regime (flow rates greater than 16 m/yr in her analysis, which is equal to 0.5 *μ*m/s). As an upper limit, Molins *et al.* (2012) assert that transport control is expected only for flow velocities less than 100 *μ*m/s, about two orders of magnitude greater than the lower limit. A similar range of flow velocities under which chemical reaction rates are transport-controlled is given by Salehikhoo *et al.* (2013).

Fig. 1 is adapted from Loague and Corwin (2006). In addition to reproducing their placement of various processes, we also plot a fluid flow rate of 2 *μ*m/s and the associated 3D solute velocity function. In each case, we give the appropriate two order of magnitude envelope in the scaling, attributable to a two order of magnitude range in hydraulic conductivity. The solute and fluid velocities in the figure are constrained to be equal at the scale of 1 *μ*m. The choice of 1 *μ*m as a fundamental scale for heterogeneity is intended to be compatible with a pore scale, without contending that this value is exact or universal.

In this context, the noteworthy features of Fig. 1 are that (1) the vegetation time and length scales (relevant to industrial agriculture: plant canopy, tillage ridge, crop) tend to follow the line of constant fluid velocity, while (2) the chemical process length and time scales associated with chemical processes (relevant to weathering, chemically limited landscape evolution, and soil development) follow the line of constant solute velocity. Most vegetation features lie between the two lines. Our interpretation is that natural processes that follow the solute velocity line are chemical-limited, while those that follow the fluid flow velocity line are water-limited. The equality of solute and fluid velocities at the scale of a single pore suggests that the relevant heterogeneity is set by mineralogical variability among the individual particles. We will see below that the fundamental length scales associated with chemical weathering are indeed frequently in the range of microns to millimeters. However, as we previously noted (Hunt *et al.*, 2011), field experiments on solute dispersion tend to introduce a relevant length scale of approximately 1 m. This fundamental length scale is traceable to the experimentally imposed solute plume, and shows up in universal scaling behavior of the dispersivity at length scales substantially longer than 1 m. Interestingly, a fundamental length scale for solute transport that is orders of magnitude larger than a pore scale can also result from relatively uniform application of chemical fertilizers over agricultural length scales.

If the fundamental scale of heterogeneity employed for Figure 1 is increased from 1 *μ*m to 1 m, the solute transport line is displaced downwards by 6 orders of magnitude and diverges from the constant flow line at about the placement of tillage ridge. Such a change in the heterogeneity scale can be effected by applying nutrients homogeneously over scales of at least meters. While 1 m of vertical soil development, tied to natural heterogeneity at roughly the micron scale, shows up at about 1000 years, tillage ridge structures, also at nearly a meter, are displaced to a time scale of less than a month, more than four orders of magnitude shorter. Of course, it is possible to dig a structure at virtually any length scale in days or less, but only surface patterns like tillage ridges are suitable to agriculture. This is then a critical advantage of industrial agriculture, which can reduce time scales for vegetation development by orders of magnitude with the simple expedient of homogenizing the nutrient delivery.

## MATERIALS AND METHODS

In order to evaluate the proposed percolation theoretical approach, we used several datasets available in the literature. For each of these datasets, we discuss here our estimation of the relevant parameters. The datasets are divided into those for which we can estimate both *x*_{0} and *t*_{0}, and those for which only *x*_{0} is known or estimated.

### Datasets for which we can estimate both *t*_{0} and *x*_{0}

#### White and Brantley (2003)

The authors reported weathering rates versus time for silicate minerals, in particular, plagioclase, K-feldspar, hornblende, and biotite, for both laboratory and field conditions (Tables IV–VII in the original article). This database covers a time range from 10^{−2} to 10^{7} years. While accurate flow rates cannot be calculated for field conditions, the article presents data for an experiment intended to correspond to field conditions (fresh Panola plagioclase or granite), at least in its general characteristics. White and Brantley (2003) crushed 750 g of granite and placed it in a column 100 cm long and 2.4 cm in diameter, leading to a bulk density ρ_{b} of 750 g/[100 cm(π)(1.2 cm)^{2}]. Using a value of ρ = 2.7 g/cm^{3} for the particle density of granite, and applying the equation ρ_{b} = (1 − *ϕ*) ρ, one finds *ϕ* ≈ 0.4 for the porosity. At the given volume flux of 10 ml/h through this column, mass conservation gives a pore-scale velocity of about 15 *μ*m/s. Given that the particle sizes in the experiment range from 0.25 mm to 0.85 mm, we estimate a typical pore radius as 0.3 times the typical particle radius (Gvirtzman and Roberts, 1991), or about 0.15 mm. The inferred pore size distribution is actually so narrow that one should not necessarily expect percolation theory to be the best means to enumerate the solute transport paths. The purpose of this calculation is to check the reasonableness of the absolute time scales inferred from the theory. With a typical pore size of 150 *μ*m, a typical pore crossing time is about 150 *μ*m/15 *μ*m/s = 10 s. With a column length of 100 cm, one can find between 1000/0.85 = 1176 and 1000/0.25 = 4000 particles along the length of the core, and a similar number of pores. Core *solute transit times* of between (10 s) (1176)^{1.87} = 0.17 yr and (10 s) (4000)^{1.87} = 1.7 yr would be expected (from Eq. (9)), comparable to the time values in the experiment that range from 0.19 yr to 6.2 yr. For establishing the pore scale reaction rate, we use the scaling relationship *R* = *R*_{0} *x*^{1−Db} from Eq. (13) to translate the measured reaction rate at the shortest time interval (for the full column length) to its value at a single pore scale. For this single calculation, we used a geometric mean (1176 × 4000)^{0.5} value for the number of pores along the column. Note that one cannot then use inferences from the ADE to extract a spatial tortuosity factor as proportional to the temporal tortuosity. The spatial tortuosity exponent in three dimensions (for conditions of full saturation; Sheppard *et al.*, 1999) is 1.43, which would yield total solute path lengths of 21 m–35 m for a 1 m long column.

For comparison with the field data, we plot the tabulated data of White and Brantley on a horizontal axis in units of years, and normalize the vertical scale to the greatest observed rate.

#### Du *et al.* (2012)

The experiment was performed in a stirred flow cell 25 mm tall with 0.5 g of particles in 10 ml of liquid—in other words, it was a suspension, not a grain-supported medium. The data were collected from elution experiments for six different samples that were sieved into 5 distinct size fractions: <75 *μ*m, 75–500 *μ*m, 500–2000 *μ*m, 20–2000 *μ*m, and two replicates for the 2000–4000 *μ*m gravel fraction. A synthetic groundwater with a pH value in the range of 7.8–7.9 was applied to study Uranium release from these different size fractions. The output data were concentration as a function of pore volume, which for constant flow rate could be assumed to be proportional to travel time. In Table I, we use the published particle sizes and the particle numbers derived therefrom to calculate mean particle separations and the flow time (neglecting effects of stirring). The mean distance between the particles was calculated by assuming the particles were equidistantly dispersed. The total flow time was calculated from the product of an interparticle flow time and the quotient of the column length to interparticle distance raised to the power 1.87. Results from two different runs of the experiment on the largest particles were compared (Fig. 2). In order to generate the same regression coefficients in both experimental results, the time scale of the second run had to be chosen a factor 1.86 larger than for the first run (9.4 h rather than 5.2 h, as given in column 7), leading to the conclusion that the flow rates were probably not identical, but carried an uncertainty of about a factor 2. For later comparison with experiment, we used first the calculated column crossing time to normalize the time axis, and subsequently the values in parentheses.

Diameter, μm
. | Number of particles in 0.5 g . | Mean distance between particles, μm
. | Flow time between particles, s . | Column crossing time, h . | ||
---|---|---|---|---|---|---|

Lower bound . | Upper bound . | Geometric mean . | ||||

20 | 75 | 38.7 | 6 202 812 | 117 | 28 | 177 (93) |

75 | 500 | 193.6 | 49 622 | 586 | 141 | 43.8 (68) |

500 | 2000 | 1000.0 | 360 | 3028 | 727 | 10.5 (34) |

20 | 2000 | 200.0 | 45 044 | 606 | 145 | 42.3 |

2000 | 4000 | 2828.4 | 16 | 8563 | 2550 | 5.2^{a} |

Diameter, μm
. | Number of particles in 0.5 g . | Mean distance between particles, μm
. | Flow time between particles, s . | Column crossing time, h . | ||
---|---|---|---|---|---|---|

Lower bound . | Upper bound . | Geometric mean . | ||||

20 | 75 | 38.7 | 6 202 812 | 117 | 28 | 177 (93) |

75 | 500 | 193.6 | 49 622 | 586 | 141 | 43.8 (68) |

500 | 2000 | 1000.0 | 360 | 3028 | 727 | 10.5 (34) |

20 | 2000 | 200.0 | 45 044 | 606 | 145 | 42.3 |

2000 | 4000 | 2828.4 | 16 | 8563 | 2550 | 5.2^{a} |

^{a}

For comparison in Fig. 3, 5.2 h was used for both runs.

#### Liu *et al.* (2008)

Liu *et al.* (2008) measured the desorption of U(VI) from contaminated sediments collected at the Hanford 300 area of the Hanford site. In a small column (10.5 cm length × 2.4 cm diameter), sediments with size of 2 mm and less were packed with a porosity of 0.41. In a large column (80 cm length × 15 cm diameter), cobbles and gravels were also included to investigate the effect of such large particles on U(VI) desorption. Applying the medium descriptions with the given pore-scale velocity, we used the geometric mean pore size to calculate fluid advection times for single pores, *t*_{0}, to be 0.08 min for the short column and 0.19 min for the long. The geometric mean pore size was inferred from the geometric mean particle size, itself derived from the particle-size distribution. *R*_{0} was calculated similarly as for White and Brantley's (2003) Panola plagioclase, exploiting the known scaling of *R* with time at small times: *R*(*t*/*t*_{0})^{0.47} = *R*_{0}. The columns were saturated from the bottom with synthetic groundwater, and then a constant flow rate was applied to leach the sediments: 8.6 cm/h for the short column and 3.52 cm/h for the long. Particularly with such coarse Hanford sediments, wall flow may not be negligible (Cherrey *et al.*, 2003), and indeed our investigations (Ghanbarian-Alavijeh *et al.*, 2012) of those arrival time distributions showed evidence of 2D flow attributable to wall flow. For the present data, comparison with experiment will indicate that wall flow under unsaturated conditions dominated the solute delivery in the short columns.

#### Zhong *et al.* (2005)

The sediment was collected from the Oak Ridge site in eastern Tennessee to investigate the microbial immobilization of groundwater U(VI). Uranium(VI) sorption to the sediment and desorption were performed with and without Fe(II) at different pH values. For the purpose of this study, we considered only the experiments that showed U(VI) desorption at pH 7. We used the initial data points to generate both the fundamental reaction rate and the fundamental time scale in each case. For the case that the iron concentration was 0 mmol/g, the initial Uranium concentration was 158 *μ*g/l at 5.3 days. For the case that the iron concentration was 0.02 mmol/g, the initial Uranium concentration was 226 *μ*g/l at 5.3 days.

### Datasets for which only *x*_{0} is known

#### Maher (2010)

The data, collected from different published papers, are chemical weathering rates of granitic alluvial materials and sea-floor sediments as a function of fluid residence time in the range of 10^{−3}–10^{5} yr, weathering rates versus flow rates, and surface ages (10^{4}–10^{6} yr) versus flow rates. A fundamental “surface controlled” reaction rate of 10^{−2.3} yr^{−1} is assumed in Maher. For analyzing her data, we appropriate this value for use as *R*_{0}. In her Figure 7, the collected data for the spatial scaling of reaction rates extrapolate to 10^{−2.3} yr^{−1} at about 1 mm. Using the flow rate of 0.05 m/yr (from Maher's Figure 7) assumed relevant for these data, and a fundamental time scale of about 10^{−2} yr (from Maher's Figure 2) one finds a fundamental length scale about 0.5 mm. We used the mean of these two values, *x*_{0} = 0.75 mm, as the fundamental length scale for this dataset.

We also compare Maher's data for the temporal scaling of reaction rates directly with the data of White and Brantley (2003), and use the actual time scales for both. However, in that comparison we scale the reaction rates from Maher to 10^{−2.3} yr^{−1}, just as in the spatial comparison.

#### Salehikhoo *et al.* (2013)

Column experiments on magnesite dissolution were performed in columns of length from 5 to 22 cm, and at various flow velocities. Results were reported in terms of the Dammkohler number, Da_{I} = *x*/(*v*_{0} τ_{R}), where τ_{R} is a reaction time and *x* the column length. For constant fluid velocity *v*_{0}, the Dammkohler number is therefore a proxy for the column length. In their experiments, magnesite and quartz were both ground to a size range of 300 *μ*m–540 *μ*m, leaving a typical particle size (geometric mean) of about 400 *μ*m. Using the commonly assumed ratio of pore size to particle size of 0.3 gives a mean pore size of about 120 *μ*m, which we take as the fundamental length scale. Salehikhoo *et al.* (2013) demonstrated that data for different column lengths and flow velocities could be reasonably collapsed to a single curve when represented in terms of the Dammkohler number. We digitized their data and chose to represent the dependence of the reaction rate on system length, in accord with the observation above that for a constant flow velocity, the Dammkohler number is proportional to the system length. This means that a factor of 5 range in Dammkohler numbers corresponds to a range of system lengths of the same factor 5 and allows us to find approximate positions for 5 cm and 22 cm on the horizontal axis, at the beginning and end of their data range.

#### Peng *et al.* (2012)

Peng *et al.* measured Uranium and Copper concentrations as functions of depth from the surface of an individual basaltic clast using laser ablation. The depths ranged from 0.6 *μ*m to 29 *μ*m. They reported that the copper and uranium concentrations were highly correlated, so we plotted them simultaneously. Our strategy for choosing fundamental length and time scales was to normalize both the copper and the uranium depths and their rates to their initial values generated at a depth of slightly less than 1 *μ*m. Presented in this manner, the two data series are nearly indistinguishable.

For more information about the datasets used in this study, the interested reader is referred to the original published papers.

## RESULTS AND DISCUSSION

The above calculations or estimations of *R*_{0} and *t*_{0} allow graphical representations of different datasets on the same plot. In order to emphasize universality, we utilize this information wherever possible. Note that we use the same theoretical function in every case, meaning that the model morphology and the fluid flow characteristics are always assumed to be the same. These particular assumptions correspond to random percolation exponents (valid for saturated conditions) and a fractal dimensionality of the pore space of 2.95. For porosities of 0.2–0.4, *D* = 2.95 corresponds to pore size ranges of between 100 and 20 000, respectively, representing fairly heterogeneous materials.

The data from Du *et al.* are shown together in Fig. 3(a). To construct this figure, we used the raw data for the concentrations, and transformed the pore volume data to a column crossing time using the values calculated in Table I, column 7. We did not assume any variation in the fundamental reaction rate (or, therefore, the initial concentration). Thus, the comparison with experiment uses two adjustable parameters, *t*_{0} and *R*_{0}. Because there is uncertainty in the flow rate of about a factor 2 (see earlier discussion of Du *et al.*, 2012), we adjusted the times of three of the six experiments by a factor of roughly 2 (the values in parentheses in Table I, column 7) to generate Fig. 1(b). This increases to 5 (three experimental, two theoretical) our total number of adjustable parameters for comparing one function with six different experimental conditions. With these 5 total parameters we can fit four values of the slope, and three times at which the slope changes for each of six different experiments, for a total of 42 different constants.

Next consider the experimental data from Zhong *et al.* (2005), White and Brantley (2003), and Liu *et al.* (2008) (Fig. 4). For these three datasets, knowing the actual flow velocity allowed us to calculate *t*_{0} quite accurately. In Zhong *et al.* (2005) and Liu *et al.* (2008), we used the characteristics of the column to generate *R*_{0} at the time scale of *t*_{0} (and scale of a single pore) by applying the simple power-law scaling of reaction rates at short times (Eq. (12)). Note that this power-law scaling would represent a slight underestimation of the initial *R*_{0} for the White and Brantley data (fresh Panola plagioclase), and hence overestimation of *R*/*R*_{0}, since the values of *t*/*t*_{0} exceed somewhat the maximum value for which the simple power law relationship is valid. Consequently, we use the full numerical solution to generate *R*_{0} from the values of *R* at the smallest times measured. For the first data set (Zhong *et al.*, 2005), we used the value of the time at which a maximum concentration was reached for *t*_{0}, and that value of the concentration for *R*_{0}. Consequently, Fig. 4 compares a prediction with zero adjustable parameters to three separate experiments.

In Fig. 5, we show together the field data for weathering rates compiled by White and Brantley (2003) and Maher (2010). *R*_{0} for Maher was given as 10^{−2.3}, while for White and Brantley we simply used the largest value measured, which for most silicate minerals was 3.6 × 10^{−11}. Note therefore that for both datasets, the largest individual values of log(*R*/*R*_{0}) are 0. The time scales for both are in years; we have not attempted to find the fundamental time scales for the hundreds of individual reaction rates tabulated. This particular comparison uses a two-parameter fit, with both *R*_{0} and *t*_{0} chosen to optimize the agreement between experiment and theory. Recall that our prediction for the reaction rates would exhibit a two order of magnitude variability at any time scale due to the uncertainty in flow rate, and this two order of magnitude variation is in exact accord with the variability in flow rates for which reactions may be considered to be transport controlled (Maher, 2010; Molins *et al.*, 2012). Under the circumstances, the agreement obtained between the theory and experiment is quite good. Navarre-Sitchler and Brantley (2007) argue that when temperature variability is accounted for, the remaining variability in reaction rates is about two orders of magnitude. Thus, one could reasonably conjecture that the variability beyond that predicted is attributable to the lack of control of the temperature variable.

In Fig. 6, we show field data and laboratory data together for the scaling of chemical reactions with distance or depth. The data that were given directly in terms of spatial variables are taken from Maher (2010) and Peng *et al.* (2012). The data from Salehikhoo *et al.* (2013), originally given in terms of the Dammkohler number, Da_{I}, are represented here as a function of length. The values of *R*_{0} and *x*_{0} for the data of Maher (2010) were discussed above in Materials and Methods. *R*_{0} and *x*_{0} for the data of Peng *et al.* (2012) were simply chosen as the largest *R* and corresponding *x*_{0} value, respectively. The choice of *x*_{0} for Salehikhoo *et al.* (2013) was also described in Materials and Methods. Note that use of the value of *R*_{0} reported by Salehikhoo *et al.* (2013) would not allow such outstanding correspondence between theory and experiment. The value we use is obtained by extrapolating the dependence *R* = *R*_{0}(*x*/*x*_{0})^{−0.87} from the smallest system sizes (5 cm) back to the dimensions of a single pore (120 *μ*m); this implies that, in our model, the length dependence of the reaction rate sets on already at 120 *μ*m, rather than at about 5 mm, as assumed by Salehikhoo *et al.* (2013). Our value for *R*_{0} is thus about 1.5 orders of magnitude higher than the value given by the authors. Consequently, Fig. 6 may be assumed to have a single adjustable parameter. But our use of this parameter implies that the reaction is transport-limited over a larger range of length scales than was assumed by Salehikhoo *et al.* (2013), a discrepancy that should be further investigated. Our hypothesis that reaction rates are proportional to the velocity, together with our result that *R* is proportional to *x*^{−0.87} for smaller distances, and a rather larger negative power at larger distances, generates *R* ≈ Da_{I}^{−1} as a very good approximation. Thus, our general results are compatible with the statement of Salehikhoo *et al.* (2013) that reaction rates can be scaled by the Dammkohler number.

While we have here shown that a large amount of data are consistent with a single prediction, it should be emphasized that not all data can be explained by 3D random percolation exponents. The exceptions that we have seen so far are all consistent with 2D invasion percolation theory, which is applicable for unsaturated flow along surfaces. The appropriate value of *D*_{b} in this case is 1.22 (Sheppard *et al.*, 1999), while ν takes the value 4/3 (Stauffer and Aharony, 1994). This yields a scaling of *R* with *t*^{(Db−1)/}^{Db} ≈ *t* ^{−0.18}. Of all the data we have investigated in detail, we have only one specific case for which this complication is observed: the short column experiment reported in Liu *et al.* (2008). In Fig. 7, we show the clear distinction between the results of the short column and those of the long column. The exponent extracted from the short column experiments is −0.17, and from the long column experiments, −0.47, almost precisely the predicted values, as was already apparent for the latter case in the comparison from Fig. 4. An analogous result was also found in the arrival time distribution in one Hanford site sediment. In that case, the interpretation was that the two-dimensional unsaturated conditions were a result of the relevance of wall flow in the case that the individual particles were large compared with the dimensions of the experimental apparatus. Note, however, that unsaturated flow and flow along fracture surfaces can also be dominant in arid environments (Glass *et al.*, 2005). Thus, it may not be surprising that a collection of field data (Sak *et al.*, 2003) for reaction rate scaling in alpine areas of New Zealand (Chinn, 1981; Knuepfer, 1994) show a weathering rind thickness dependence *x* ≈ *t*^{0.8}, which our theory would predict to be *x* ≈ ∫*v*(*t*)d*t* = ∫*t*^{−0.18}d*t* = *t*^{0.82}. The difference between the observed exponent (0.8) and its theoretical value (0.82) is probably insignificant.

In Figure 8, we compare the predicted and observed dependences of reaction rates with the scale of the measurement. We also compare these two results with those of a theory that predicts a fractal reaction front, with a fractal dimensionality of 2.33 as a fitted parameter. In percolation theory, there are two contributions to the perimeter of a cluster: one is proportional to the square of the cluster size and the other is proportional to the volume of the cluster. The first term provides the horizontal asymptote for the reaction rate at small cluster sizes. But percolation theory yields the mass fractal dimensionality (≈2.5) of clusters near the percolation threshold, so that the fractal dimensionality is not a fitting parameter; the percolation model therefore has one fewer adjustable parameter than the fractal model.

We now consider why our particular results, derived for solute transport, unite in understanding the various phenomena discussed, and why this could not have been accomplished otherwise. The most important factors are (1) that the solute velocity is proportional to the flow velocity, (2) the time dependence of the solute velocity resembles the dependence from diffusion closely (under saturated conditions of 3D flow connectivity), and (3) is referred to a fundamental length scale of heterogeneity. The resemblance to diffusion coupled with the experimental results triggers association with the ADE. But the proportionality of reaction rates also to the flow velocity, as revealed in analysis of the experiments of Salehikhoo *et al.*, 2013 in terms of the Dammkohler number, leads to confusion if a diffusion mechanism is invoked. Further, an interpretation of the time scales in Figure 1 in terms of a flow velocity for water-limited processes, and diffusion for chemical-limited processes, would not have united these processes at the length scale of a pore. Such an interpretation would also never have revealed the role of homogenization of solute delivery relevant to industrial agriculture. It is the fractal structure of the critical (and related) flow paths near the percolation threshold that govern the solute velocity. Homogenization of the solute delivery to any particular length scale larger than the natural heterogeneity leads to the elimination of the relevance of fluctuations in solute concentration at smaller length scales. Reduction of the size range of relevant fractal structures impeding transport diminishes the tortuosity of the dominant paths, increasing the speed of the solute transport. Finally, diffusion, as a means to overcome the limitations of solute gradients, would never, in the same theoretical framework, generate a change of solute velocities at larger times to dependences much slower than diffusion.

Thus, in order to simultaneously explain all the results reported here, a theory must have all the features that our derivation supplied. Note that all these features were present in the original derivation; since the original publication (Hunt and Skinner, 2008) neither the theory nor the results have changed. The same solute distributions have been used to predict the arrival time distributions of solutes as a function of saturation (Ghanbarian-Alavijeh *et al.*, 2012), the dispersivity as a function of length scale (Hunt *et al.*, 2011), the variance as a function of time (Hunt and Skinner, 2010b), and the scaling of the typical arrival time with system length (Hunt *et al.*, 2011).

We note that weathering rates in other circumstances have also been studied frequently, particularly for rocks lying on the earth's surface (surface clasts). Summaries of the development of the weathering rinds on these rocks, all from Alpine New Zealand were given in Sak *et al.*, 2008. For those cases (Chinn *et al.*, 1981; Knuepfer *et al.*, 1994), which were assumed to be power law in nature, i.e., rind thickness *x* ≈ *t ^{δ}*, the extracted values of

*δ*were 0.76, 0.82, and 0.86, with mean 0.81. The simplest prediction, from Eq. (9) is that

*x*≈

*t*

^{1/Db}. Use of

*D*

_{b}= 1.217 (Sheppard

*et al.*, 1999), appropriate for 2-D flow paths through unsaturated media, generates

*δ*= 1/

*D*

_{b}= 1/1.217 = 0.82, extremely close to the mean measured value. For further confirmation, we took data from four experiments on the weathering rinds of basalt clasts in tropical, or otherwise highly vegetated regions (Sak

*et al.*, 2003; Pelt

*et al.*, 2008; Ma

*et al.*, 2012; and Oguchi and Matsukura, 1999), and plotted the logarithm of the weathering rind thickness as a function of the logarithm of time (Fig. 9). The extracted power in the power-law relationship was 0.688, which compares well with the theoretical value of 1/

*D*

_{b}= 1/1.458 = 0.686 for 3-D invasion percolation (Sheppard

*et al.*, 1999).

## CONCLUSIONS

We have shown that the available data for the scaling of surface reactions in porous media conform closely to the predictions of our transport theory. The comparisons used a minimum of adjustable parameters, and sometimes none. The implication is that the observed reaction-rate scaling is a consequence of transport limitations, whether on the availability of reagents or the removal of products. In the context of percolation theory, the evidence from experiment implies that the solute transport is chiefly through the bulk medium under saturated conditions. In one experimental case and in some field summaries, evidence for unsaturated two-dimensional flow conditions controlling reaction rates was found. In an experiment, such flow conditions could be caused by incompletely saturated wall flow through coarse sediments. Field results indicating similarly unsaturated flow conditions, obtained in Alpine regions, suggest flow along fracture surfaces: flow conditions within fractures are much more frequently unsaturated than in the surrounding matrix.

We emphasize that we have avoided discussion of the details of the reactions themselves. We have sought the guidance of unifying, rather than distinguishing, characteristics. As an example, we have not addressed the distinct influences of chemical gradients in the mineralogy and those in the fluid. Ignoring such differences means, for example, that our upscaling technique as presently developed does not address complexities that can develop from secondary deposition of reaction products.

Our conclusions for the media and experiments considered are:

Weathering reaction rates are fundamentally transport limited.

The relevant type of transport limitation is advective.

An increase of reaction rate with measurement scale is generated by the known increase in surface area of percolation clusters; the fractal dimensionality of the surface is equal to the mass fractal dimensionality of percolation theory.

Reaction rate scaling (in time and space) can be predicted by scaling the velocity of conservative solute transport.

The solute transport scaling cannot be predicted from conventional continuum mechanics modeling (i.e., the Advection Dispersion Equation, ADE), but must be predicted using the percolation theoretical framework.

The observed scaling exponents in space and time relate simply and directly to the fractal dimensionality of the percolation backbone.

Known variability of fluid velocities and previously published limitations on the range of fluid flow velocities for transport-controlled reactions reproduce very well the observed range of data for field weathering rates over 11 orders of magnitude in the time scale.

The same dependence of solute transport on fluid flow rates and fundamental length scales of heterogeneity organizes existing graphical summaries of length and time relationships in subsurface processes (Blöschl and Sivapalan, 1995; Loague and Corwin, 2006).

Massive use of fertilizer sets a much larger fundamental length scale for solute transport in agriculture than is derived from natural heterogeneity, accelerating agricultural time scales by many orders of magnitude.

It is interesting that similar effects of homogenization on nutrient transport are found in streams (Helton *et al.*, 2011; Petersen *et al.*, 2001), whereas the homogeneous nutrient delivery to crops reflects their nearly homogeneous delivery by fluvial processes to the near-shore environment of the ocean. Increasing productivity of the land surface through industrial agriculture leads to greater nutrient loading of the streams. Sufficient overloading of the streams produces an overloading of the ocean in any case. But homogenization of, in particular, river headwaters, contributes to a more rapid delivery of the increased nutrient load to the oceans, and increased oceanic productivity. The transport of nutrients from river headwaters to trunk streams involves repeated nutrient recycling, as is also possible in reactive contaminant transport in soils. But we have seen that the temporal dependence of the chemical reactions in porous media is not primarily a consequence of their adsorption and desorption; rather it is a consequence of the relevance of the fractal paths that the solute follows. The fractal dimensions of percolation theory are relevant to these paths because of the tendency for the advecting fluid, i.e., water, to find the paths of minimum resistance. But a similar argument can be applied for nutrient transport along streams, which is likely to be controlled by the tendency for water to find the path of minimum resistance through the heterogeneity of the substrate. Whatever delays are possible due to interactions with biota, the number of such relevant interactions will increase dramatically if the path of least resistance is fractal. But the range of length scales over which percolation governs the transport is reduced with homogenization of the streambed. Thus similar changes, drastic reduction in nutrient delivery times, in solute transport in both soils and streams may result from analogous physical constraints.

## ACKNOWLEDGMENTS

This research was supported in part by NSF grant EAR 142806.

## References

_{2}Se

_{3}

_{2}Se

_{3}