Solid phase evolution in the Biosphere 2 hillslope experiment as predicted by modeling of hydrologic and geochemical fluxes

Abstract. A reactive transport geochemical modeling study was conducted to help predict the mineral transformations occurring over a ten year time-scale that are expected to impact soil hydraulic properties in the Biosphere 2 (B2) synthetic hillslope experiment. The modeling sought to predict the rate and extent of weathering of a granular basalt (selected for hillslope construction) as a function of climatic drivers, and to assess the feedback effects of such weathering processes on the hydraulic properties of the hillslope. Flow vectors were imported from HYDRUS into a reactive transport code, CrunchFlow2007, which was then used to model mineral weathering coupled to reactive solute transport. Associated particle size evolution was translated into changes in saturated hydraulic conductivity using Rosetta software. We found that flow characteristics, including velocity and saturation, strongly influenced the predicted extent of incongruent mineral weathering and neo-phase precipitation on the hillslope. Results were also highly sensitive to specific surface areas of the soil media, consistent with surface reaction controls on dissolution. Effects of fluid flow on weathering resulted in significant differences in the prediction of soil particle size distributions, which should feedback to alter hillslope hydraulic conductivities.


INTRODUCTION
The goal of the Biosphere 2 (B2) institutional hillslope experiment is to help develop an understanding of how climate, rock and biota interact over time to influence hillslope structure formation, and its impact on the movement of water in the landscape (Huxman et al., 2009). In order to address the coupled processes involved, three replicated, zero-order watersheds (hillslopes) were designed (Hopp et al., 2009) and are now under construction. The hillslopes will comprise a loamy-sand textured granular basalt as the starting material. Basalt lithology was chosen because of (i) the global significance of basalt to CO 2 consumption in silicate weathering (Suchet et al., 2003) (ii) the relatively high rate of chemical weathering for basalt mineral assemblages (Nesbitt and Wilson, 1992), and (iii) the large mass concentrations of nutrient cations whose dissolution can effectively support microbial and plant life (Hinsinger et al., 2001).
One process coupling of particular interest is that between hydrologic flow path evolution and mineral weathering. Since rock-forming primary silicates are thermodynamically unstable at Earth surface conditions, they are transformed to more stable secondary minerals, typically of smaller particle size and lesser crystallinity (Nordstrom and Munoz, 1994). The rates of these reactions and the nature of solid phase products are dictated by the through-flux of fresh water and dissolved or gaseous dioxygen (O 2 ). The rate of mineral weathering depends on the types of minerals present initially, in addition to temperature and water flux rate (e.g., White and Brantley, 2003;White and Blum, 1995).
An important question is whether significant structural change will occur in the hillslope porous media over the assumed 10 year lifespan of the experiment. A key component of structural change is the incipient dissolution of the initial primary mineral assemblage, and the resulting precipitation of neoformed, "secondary" solids. These primary to secondary mineral transformations (termed "incongruent" weathering reactions when they result in the formation of solid-phase mineral products) tend to shift the particle size distribution to smaller values because the newly-formed, clay-sized particles accumulate at the expense of primary sand and silt sized materials (Sposito, 2008). Increasing the proportion of fine particles generally increases water holding capacity and decreases saturated hydraulic conductivity. While prior chronosequence studies have indicated the effects of long-term pedogenesis on the hydrologic behavior of soils (Lohse and Dietrich, 2005), the B2 hillslope experiment provides an opportunity to observe these coupled processes in real time as affected by well-constrained climatic and vegetative forcings.
Since dissolution-precipitation reactions are rate-limited, they tend to occur nonuniformly in the hillslope, with kinetics being dependent on the local relative saturation of the aqueous phase with respect to both dissolving and precipitating mineral solids, i.e., the "chemical affinity" (Oelkers, 2001). As chemical weathering drives the particle size distribution to smaller values, this will also affect the hydraulic conductivity and surface chemistry of the soil medium, with potentially important feedbacks on subsurface flow path evolution, nutrient and organic matter retention, etc. For example, incongruent dissolution of primary silicates to form layer silicate clays and oxyhydroxides can facilitate plant growth by (i) solubilizing many lithogenic nutrients required for microbial and plant growth, as well as by (ii) producing additional reactive surface area for the adsorptive retention of nutrient cations (e.g., Ca 2+ , Mg 2+ , and K + ) and organic molecules (humic substances and biomolecules). In addition, this incongruent dissolution can (iii) affect the development of subsurface heterogeneity (e.g., soil profile horizonation) which, in turn, (iv) feeds back to affect water transport paths and plant establishment along the hillslope.
Such biogeochemical-hydrologic couplings are considered fundamental to the evolution of structure in the Earth's "Critical Zone" (Chorover et al., 2007;Brantley et al., 2006), defined as the region of the Earth's surface extending from the outer periphery of the vegetation canopy to the lower limits of groundwater (NRC, 2001). The evolution of Critical Zone structure and function in response to climatic forcing is a principal focus of Earth science research at Biosphere2.
Our present objective is to predict the time scale over which significant chemicalweathering-induced changes will occur in the B2 hillslopes. In this contribution, we use the reactive transport code CrunchFlow2007 (Steefel, 2008) to evaluate the effects of geochemical weathering on the evolution of the soil mineral assemblage and the associated particle size distributions. A number of reactive transport models can be used to predict chemical weathering in soils and bedrock, e.g., SAFE , WITCH (Goddéris et al., 2006), and CrunchFlow (Steefel and Van Cappellen, 1990;Steefel and Maher, 2009;Maher et al., 2006;Maher et al., 2009;Soler and Lasaga, 1998). CrunchFlow2007, which is the result of 20 years of development in multicomponent reactive transport software, was chosen for the simulations because of its successful use in the past for modeling of chemical weathering and water-rock interaction (Steefel and Van Cappellen, 1990;Steefel and Lasaga, 1994;Steefel and Lichtner, 1998;Soler and Lasaga, 1998;Steefel and MacQuarrie, 1996;Giambalvo et al., 2002;Maher et al., 2006;Maher et al., 2009;Steefel and Maher, 2009).
The code CrunchFlow2007 solves the multicomponent reactive transport equation given by where φ is the porosity, C is the concentration of the ith primary species, D is the molecular diffusion coefficient in porous media, u is the velocity vector, and R r and R m are the homogeneous (aqueous phase) and heterogeneous (mineral) reaction rates respectively. The symbol ν is the stoichiometric reaction coefficient  We explore the impacts of initial particle size (specific surface area of primary minerals) and climatic regime (translated to flow rate and saturation degree) and assess the impact of chemical weathering on the time and space evolution of hydraulic conductivity. The study consists of two distinct parts: (1) an analysis of model sensitivity to selected parameters that are expected to influence weathering in the hillslope experiment and (2) predictions of hillslope geochemical development under two proposed climatic regimes. The sensitivity analysis includes one-dimensional simulations evaluating effects of mineralogy, particle size (surface area), flow velocity, and dioxygen concentration on dissolution and precipitation of minerals.
Two-dimensional simulations of a hillslope cross-section were then used to predict changes in the constructed B2 hillslope soils.

Mineralogical and Geochemical Characterization
A black basaltic cinder material, obtained from the CEMEX Corp. Flagstaff, AZ, has been identified as a candidate material for hillslope construction. In order to provide a solid geochemical foundation for the reactive transport modeling, this material was characterized for chemical and mineral composition and physical properties using a range of analytical methods.
Total elemental composition was determined using lithium metaborate/tetraborate fusion followed by inductively coupled plasma mass spectrometry analysis (ICP-MS, Activation Laboratories Ltd., Ancaster, Ontario). Elemental composition of individual minerals was determined using microprobe analysis (Cameca SX50 electron microprobe, Gennevillier, France). From these data, mass fractions of the minerals were calculated using ModAn, a computer program for estimating mineral quantities based on bulk composition (Paktunc, 2001).
These concentrations were verified by quantitative X-ray diffraction analysis (XRD) using a PANalytical X'Pert Pro MPD X-ray Diffractometer (PANalytical B.V., Almelo, the Netherlands) with Cu-Kα radiation (R 2 = 0.88) that gave bulk values for the composition of the crystalline minerals and by image analysis (ArcView GIS 3.2, Environmental Systems Research Institute, Inc., Redlands, California) of element distribution maps obtained by microprobe ( Figure 1) (R 2 =0.90).
Porous media particle size distribution was measured using a Beckman Coulter LS 13 320 Laser Diffraction Particle Size Analyzer (Beckman Coulter, Inc., Fullerton, CA). Particle density , was determined using the pycnometer method (Flint and Flint, 2002). The solid density is a direct input parameter in CrunchFlow2007, while the particle size distribution is indirectly used via the reactive or physical surface area. Surface area was measured by N 2 adsorption (Brunauer et al., 1938) using a Beckman Coulter SA 3100 Gas Adsorption Surface Area Analyzer (Beckman Coulter, Inc., Fullerton, CA).
Particle size distribution and surface area were measured on the basaltic material as received and on the same material following reconstitution for a loamy sand texture that emerged as preferred texture from hydrologic simulations (Hopp et al., 2009 As discussed in the Results section, model output indicated precipitation of the secondary minerals: Fe(OH) 2 , kaolinite, smectite, gibbsite, pyrolusite, and goethite, depending on conditions. Hence, also listed also in Table 1 are the data for log K SO , E a , and log k m for the following dissolution/precipitation reactions (Wolery et al., 1990): Dissolution rate constants (k m ) were obtained from the literature for the minerals measured in the basalt, and they are listed in Table 1. For data selection, preference was given to laboratory studies on the specific minerals found in the Flagstaff basalt. While laboratory rates tend to be greater than those derived from field data, weathering rates in the B2 hillslopes are expected to fall somewhere in between . When installed as a poorly consolidated material devoid of preferential flow paths (at t = 0), the B2 porous media will present higher interfacial area and more reactive (fresh unweathered) surfaces to percolating water than is observed in typical field systems. In addition, field rates are thought to be lower than laboratory rates, in part, because of diminished interfacial area and weathering-passivated surfaces that contact mobile water in structured soils.
Literature k m values obtained at T other than 298K were scaled using the Arrhenius relation (Lasaga, 1984): where E a is activation energy, k m,298 is the rate constant at 298K, R is the gas constant, and T is temperature in Kelvin scale. Mineral densities used for mineral mass conversion to volume fraction were 2.69 g cm -3 for labradorite, 3.27 g cm -3 for forsterite, 3.4 g cm -3 for diopside, and 2.4 g cm -3 basaltic glass (http://webmineral.com/).

CrunchFlow2007 simulations
We used fluid flow output data from 3D HYDRUS simulations (Hopp et al., 2009) of the "Lucky Hills" (Scott et al., 2000) and "Sky Island" (Brown-Mitic et al., 2007) climatic conditions as input hydrologic parameters for our CrunchFlow2007 modeling. The CrunchFlow2007 modeling objectives were (i) to determine the influence of selected physicochemical parameters on mineral transformation rates in 1D basalt soil columns, and (ii) to assess the influence on mineral weathering of 2D fluid flow in the vadose and saturated zones, including feedbacks to hydrologic flow. The 2D simulation represents the central axis of the zero-order hillslope studied by Hopp et al. (Hopp et al., 2009).
Simulation duration was 1000 days, with a steady state water infiltration for that period that gives a total rainfall equivalent to that expected over the course of the 14 to 18 year (depending on selected climate) B2 hillslope experiment. A porosity value of 0.41 was used based on an average value for the loamy sand textural class (Carsel and Parish, 1988 , Mn 2+ ) were selected to account for all major elements in dissolving minerals. A list of auxiliary or secondary complexes to be included in the calculation of total concentrations (a total of 105 species, Attachment 2) was obtained by "sweeping" the EQ3 database for all relevant species. For both experiments, inflow solution was pure water in equilibrium with ambient concentrations of atmospheric gases, including CO 2 and O 2 (to simulate the reverse osmosis purified, atmosphere pre-equilibrated water that will be used in the experiment).
Boundary conditions were Dirichlet at the soil surface (to allow for gas diffusion) and flux at the bottom. CrunchFlow2007 allows one to either define constant gas flux across the domain or to impose diffusion limitations. Both options were used in 1D simulations, to provide a comparison, whereas diffusion-controlled gas flux was used for 2D simulations. When used, constant gas flux was set at 20 m h -1 to ensure O 2 was present in excess. The gas diffusion coefficient in basalt, D b (0.0089 cm 2 s -1 ), in 1D simulations was calculated using the formula suggested for soils by Moldrup et al. (2000): where D air is the gas diffusion coefficient in air, 0.219 cm 2 s -1 for O 2 (Welty et al., 1984), and θ air is air-filled porosity (0.0615 for 0.41 total porosity and 85% saturation). For 2D simulations, the value of 0.01 cm 2 s -1 was used, since this value lies between the calculated maximum gas diffusion coefficient in dry basalt (0.059 cm 2 s -1 ) and the zero gas diffusion coefficient operative under fully saturated conditions.
The initial condition contained pore water equilibrated geochemically with basalt.
Discretization was 30 1-m segments in 1D simulations and 30 1-m segments in the down-slope direction and 11 0.091-m segments in surface to bottom direction in 2D simulations.
Dissolution rates in CrunchFlow2007 (Steefel, 2008) are described using: where R is dissolution rate in mol m -3 s -1 , A m is mineral bulk surface area, m 2 m -3 , k m is the dissolution intrinsic rate constant in units of mol m -2 s -1 , E a is the activation energy (kJ mol -1 ), Q is the ion activity product for the mineral-water reaction, K SO is the corresponding solubility product constant, and ∏a i n is a product representing the inhibition or catalysis of the reaction by various ions in solution raised to the power n. Rate dependence on reaction affinity, g, (or Gibbs energy) is defined by the parameters m 1 , m 2 , and m 3 following relationships observed by Burch et al. (1993) and Hellmann and Tisserand (2006).
Reactive surface area for primary minerals was calculated from the measured specific surface area for the basalt and mass fraction of each identified mineral phase. In the model, therefore, bulk surface area decreases linearly with the volume fraction of the dissolving mineral phase. For secondary minerals an initial bulk surface area of 100 m 2 mineral/m 3 porous medium was assumed for lack of any specific data on how nucleation and precipitation in its early stages occurs. Once precipitation occurs, surface area is recalculated so that it increases with a 2/3 dependence on the secondary mineral volume fraction.
For the feldspar dissolution reactions in our study, a kinetic model with two parallel rates was used (Hellmann and Tisserand, (2006)): where k 1 represents a faster rate that dominates far from equilibrium and k 2 results from a slower there is some evidence that kaolinite dissolution may be governed by a different rate law (Yang and Steefel, 2008). In the present work, the same rate law was used for all secondary minerals.
Generally speaking, while there is evidence that dissolution and precipitation of minerals is well described by the TST model far from equilibrium (large Gibbs free energy), dissolution close to equilibrium is not predicted as well. However, few papers describe and quantify this phenomenon (Burch et al., 1993;Hellmann and Tisserand, 2006;Maher et al., 2009). Including Al-inhibition of dissolution rates has been shown to improve predictions compared to TST model for several minerals (Gautier et al., 1994;Wolff-Boenisch et al., 2004;Gérard et al., 1998).
However,  showed that a kinetic model with two parallel rates (one nearly linear, one strongly nonlinear) following the studies of Hellmann and Tisserand (2006) and Burch et al., 1993) gives similar results to the Al-inhibition model. The Hellmann-Tisserand kinetic model was used in this study.

One-Dimensional Simulations
Initial simulations were of a 1D, 30 m length soil column. A 30 m column represents the longest flow path possible on the modeled hillslope. This allowed us to follow changes to the material occurring both close to the surface and further down-gradient at greater depth.
Orientation of the column did not affect results of the simulations as flow direction was defined independently. Five different constant flow velocities (0.00005, 0.00026, 0.00041, and 0.005 m h -1 ) were evaluated, all at 85% saturation. These velocity values were within the range observed in HYDRUS-3D simulations of water transport in the B2 hillslopes for the loamy sand texture class (Hopp et al., 2009). Specific surface areas of the porous medium in the model included that measured using N 2 BET in our laboratory for the original loamy sand (3.24 m 2 g -1 ), but several additional lower values (10, 100 and 1000-fold lower than the measured N 2 BET) were also used, since prior studies have shown that N 2 BET may overestimate reactive mineral-water interfacial area (Navarre-Sitchler and . Although same specific surface area values were applied to all constituent minerals in a given run, volumetric abundances of each mineral were based on measurement (see Table 1). Running the model across a range in pore water velocity and specific surface area allowed us to evaluate sensitivity to these parameters.

Two-Dimensional Simulations
For the 2D simulations, a rectangular model was constructed computationally to simulate the hillslope central cross-section parallel to lateral water flow (i.e., at the channel). The flow velocity and saturation profile were generated by using steady-state conditions resulting from HYDRUS-3D simulations. As the HYDRUS-3D simulation domain consisted of a grid parallel to the channel bottom, conversion to the rectangular domain only required correction of the velocity vectors. Results of the CrunchFlow2007 simulations were then plotted in the HYDRUS-3D domain.
A constant precipitation rate was imposed that was set equal to the mean value for all days that received rainfall each year for two different climatic locations (Hopp et al., 2009 56), 1000 days of rain required fewer years. Figure 2 presents the velocity distribution and saturation in the profile for the two climate regimes.
The simulations contain several simplifications relative to the actual hillslope experiments being constructed at B2. In particular, we have simulated steady-state conditions, rather than the wet-dry cycles characteristic of actual rainfall-runoff events, in order to better enable calculation of geochemical transformations. The influence of wet-dry cycles will be explored in future modeling studies. Contribution of evapotranspiration to water flow was also ignored. Further, we have assumed here that only abiotic water-rock interaction is occurring, and we have not included the potential contribution of plants and microorganisms, or the indirect effects of bio-produced organic acids, on the weathering process (Berner et al., 2003;Hinsinger et al., 2001;Neaman et al., 2005) . Furthermore, we have assumed that the hillslopes would be initially homogeneous and isotropic, although we recognize that deposition of the granular material will likely result in some degree of particle sorting, which was not considered here.
To assess the feedback effects of mineral transformation on hillslope hydrology, we estimated hydraulic parameters from soil texture evolution using Rosetta software (Schaap et al., 2001). Rosetta is a collection of artificial neural network models that estimate saturated hydraulic conductivity and parameters in the van Genuchten equation from easily obtained soil properties such as soil texture and bulk density. Rosetta was calibrated on a database of 2134 samples containing texture, bulk density, water retention and saturated conductivity data.
Correlations and their performance metrics can be found in Schaap et al. (Schaap et al., 2001).
In our Rosetta simulations sand, silt and clay content were used as inputs. It was assumed that all secondary mineral precipitation contributed to an accretion of the clay fraction, consistent with distribution of secondary minerals in the soils. Figure 3 presents a schematic of the relations between modeling components used in this study to predict solid phase evolution.

One-Dimensional CrunchFlow2007 Simulations
1D For the soil column configuration, the predicted extent of weathering in the top 30 cm was more sensitive to particle surface area than to flow velocity ( Figure 5). There was a large increase in dissolution as specific surface area (SSA) increased from 0.3 to 3 m 2 g -1 , the highest value being that which was measured for the loamy-sand texture class of material. In general, diffusion limitations on O 2(g) availability appear to diminish mineral transformation rates, presumably because of their effect on oxidation of Fe(II), which would tend to increase the magnitude of the chemical affinity of term (Eq. 4) for Fe(II)-bearing solids. The strong effect of SSA and the weak effect of water velocity indicate that chemical dissolution kinetics rather than solubility or physical transport were rate limiting (Steefel, 2007).
The fact that dissolution-precipitation curves are not mirror images of each other reflects the time-dependent depletion of primary solid phase mass by solute efflux. As a result, the negative volumetric change in primary minerals is greater than the positive volumetric change in secondary minerals. Trends shown here are specific for basalt mineralogy that is high in weatherable glass. Results from preliminary simulations of different primary mineral assemblages with lower log K SO (not shown) exhibited significant effects of flow rate on dissolution, indicating solubility, rather than a kinetic or surface reaction rate limitation on dissolution, in contrast to the dissolution rate control noted above. In addition, these latter simulations also showed a much stronger impact on mineral transformation of O 2(g) diffusion limitations.

Two-Dimensional CrunchFlow2007 Simulations
Results from the 2D simulations facilitate visualization of mineral transformations and their distribution in response to the more complex hillslope cross-section flow system (Fig. 2).
Similar to the 1D simulations, 2D secondary mineral phase distribution and their respective saturation indices (log Q/K SO ) show smectite as the principal neophase formed in both climatic regimes, particularly in the vadose zone ( Figure 6). Other secondary minerals shown are those for which CrunchFlow2007 identified the solution phase as being close to equilibrium saturation, based on aqueous phase speciation. Local water filling of pores and diminished access to oxygen appeared to strongly impact the saturation index and mineral distribution results through its impact on Fe and Mn oxidation states. For example, solutions in the groundwater zone were strongly undersaturated with respect to goethite (α-Fe III OOH) and pyrolusite (Mn IV O 2 ), whereas vadose zone solutions were at equilibrium or supersaturated with respect to these same phases.
Conversely, solutions were near equilibrium with respect to ferrous hydroxide (Fe II (OH) 2 ) in the saturated zone.
Since the net transformation of primary to secondary minerals is a principal driver for soil textural evolution, the total volume fraction of primary minerals dissolved and secondary minerals precipitated provides insight into the hillslope reaction fronts as a function of space, time and climate. In When hillslope weathering rates were compared for several surface area scenarios ( Figure 9), it was observed that total dissolution (i.e., hillslope chemical denudation) increased with surface area in a trend similar to that observed in 1-D simulations. The one exception was that Mg release decreased at the highest specific surface area, consistent with enhanced precipitation of Mg-containing smectite driven by greater supersaturation. Calculated denudation rates for the measured mineral surface area (3.2 m 2 g -1 ) were in line with field observations (White and Blum, 1995).

Rosetta Predictions
Saturated hydraulic conductivity (K sat ) was the hydraulic parameter most sensitive to changes in the basaltic media particle size distribution. All areas of the hillslope are predicted to show a change in K sat because of incongruent weathering processes ( Figure 10). Throughout the hillslope, K sat was decreased from 0.044 m h -1 , its original calculated value. Consistent with the weathering zonation, a larger decrease in K sat is observed in the upper layers of the domain, with the saturated zone retaining higher K sat values.
Differences between the two climate regimes were also observed. The Sky Island climate, with higher rainfall showed a greater reduction in the saturated hydraulic conductivity, particularly in the backslope and footslope positions, but the Lucky Hills climate showed a reduction throughout a larger areal extent, consistent with a more extensive vadose zone.

DISCUSSION
A number of assumptions were made in order to facilitate transfer of hydrologic information from HYDRUS into the CrunchFlow2007 reactive transport code. Perhaps most significant was the assumption of a steady-state flow regime. Although this better enabled stable and computationally efficient geochemical calculations for reactive transport simulation, it is not a condition likely to develop on the natural hillslope. Nonetheless, we contend that this model configuration permits a straightforward assessment of model response to input parameterization.
The effects of transient seasonal flow will be an interesting subject for future investigations.
Results of the reactive transport model are clearly sensitive to input parameterization, as indicated by the effects of specific surface area (SSA) and climatic forcing on mineral transformation rates and feedbacks to hydraulic conductivity parameters. In particular, the calculation of dissolution rates (Eq. 3) is dependent on quantifying mineral-specific SSA values, but the estimation of reactive interface in soils is one of the greatest sources of uncertainty (White and Blum, 1995). Given this uncertainty, we used laboratory measured values to provide an upper limit for geochemical reaction modeling, but also explored the influence of decreasing SSA by several orders of magnitude. Justification for this latter approach derives from the widespread observation that laboratory dissolution rates considerably exceed those observed in the field, and this difference has been attributed, in part, to overestimates of mineral surface area in contact with percolating water in field soils (Steefel, 2007;Drever and Clow, 1995). For this reason, surface area has been used as an adjustable parameter to reconcile differences between lab and field rates of weathering (Navarre-Sitchler and . Another proposed reason for discrepancies between laboratory-scale and field-scale weathering rates is the aging of mineral surfaces and the passivation that results from formation of protective secondary layers (Drever and Clow, 1995). However, as discussed in the Materials and Methods section, the B2 hillslope system, comprising fresh primary mineral surfaces and a relatively homogeneous, isotropic porous medium, more closely represents a laboratory experiment than a field-weathered soil. For this reason (absent experimental data on hillslope scale weathering rates for adjustment), the measured N 2 BET SSA values of the basalt were used for the 2D simulations with the understanding that predicted mineral transformations may be higher than those observed in the experiment.
The quantitative model assessments of geochemical weathering and K sat change in the B2 hillslopes must therefore be seen as hypothetical; they provide a means for postulating the timescale of structural change, the locations of reaction fronts in relation to water table dynamics, etc.
The B2 hillslope geochemical model will be increasingly constrained through data acquisition over the course of the experiment and this is expected to lead to an iterative and continuous evolution of the quantitative assessments. Model results are intended to facilitate prediction of where and why to expect accumulation of secondary minerals and to aid in the design of sampling strategies to capture and quantify changes within the hillslope. To the best of our knowledge this is also the first study that predicts spatially-resolved changes in hydraulic properties of a catchment as a result of geochemical weathering.
Dissolution can be either surface-reaction or transport-controlled depending which rate is slower (Steefel, 2007). When weathering is controlled by chemical kinetics, the solution may be undersaturated with respect to dissolving mineral(s) but slow surface reaction rate impedes weathering nonetheless. Under transport control, the aqueous phase quickly equilibrates with the dissolving solid, decreasing the magnitude of the chemical affinity term (log Q/K SO ), and thereby diminishing the dissolution rate until solutes can be transported away from the dissolving solid by advection or diffusion. Lack of flow velocity effects on dissolution of primary minerals in the 1D simulations indicated that weathering kinetics were controlled by mineral dissolution rate under the modeled conditions.
Conversely, the 2D simulations revealed indications of both transport and dissolution control. A 10-fold decrease in particle SSA resulted in a decrease in the extent of dissolution and precipitation on the hillslope cross-section, consistent with surface reaction rate limitations (not shown). However, the distinct pattern of changes in the amount of primary and secondary minerals could be clearly related to the flow patterns on the hillslope cross-sections, and was different between two climate regimes indicating transport control. These results support previous studies indicating that within a watershed, both transport and dissolution extend control over weathering front advance (Steefel et al., 2005;Steefel, 2007). While on the large scale, overall progress of the weathering front is controlled by transport, as can be observed by less weathering under dry climates, local control can derive from the chemical dissolution rate of constituent minerals, which depends on the reactive surface area and site density of the particles.
In agreement with previous experimental studies (Oelkers and Gislason, 2001;Gislason and Oelkers, 2003;Pokrovsky and Schott, 2000;Nesbitt and Wilson, 1992), both in 1D and 2D simulations, basaltic glass was the main phase undergoing dissolution, along with a smaller contribution from forsterite (olivine). It is known that the relatively high solubility of basaltic glass leads to stabilization of more crystalline igneous minerals (Gislason and Arnorsson, 1993).
For example, release of Si, Al, Ca, Mg, Fe from the glass reduced the extent of undersaturation with respect to olivine and induced supersaturation of solution with respect to pyroxene.
Consequently, the dissolution of basaltic rock constituents is typically sequential, depending on the relative amount and composition of the igneous minerals and glass present, with glass dissolving first, then olivine, and finally pyroxene and plagioclase.
Of the common secondary minerals that form from incongruent weathering of basalt that were selected for inclusion in CrunchFlow2007 simulations, smectite precipitated in the largest amounts. Initial precipitation of smectite can be explained by the high Si to Al ratio of the dissolving glass. The 2:1 layer silicate mineral smectite has a higher Si/Al ratio than the 1:1 layer silicate kaolinite and is, therefore, more stable in the initial stages of basalt weathering (Chadwick and Chorover, 2001). Smectite is commonly observed as a result of weathering of basaltic material, particularly in the presence of glass (Tosca, 2009;Singer, 2007). Deeper in the profile, particularly for the Sky Island climate with thicker saturated layer, under oxygen-limited conditions, iron (II) hydroxide was found to precipitate. Development of deeply anoxic conditions resulting in Fe(OH) 2 precipitation is unlikely in semi-arid climatic regimes, but it may occur as a metastable phase during the rainfall season, which occurs in both climates studied.
The weathering-induced evolution of heterogeneity in hydrologic properties of the initially homogeneous hillslope is consistent with our current knowledge of spatial variability in the natural hillslopes. On the scale of the soil profile or pedon, a finer particle size and smaller hydraulic conductivity are usually observed in association with the weathering front. The decrease in particle size promotes increased water holding capacity, ion exchange capacity, nutrient retention, and, therefore, an improved environment for plant growth. The model simulation indicates that the surficial layer of basalt that resides in the unsaturated zone is subjected to more primary mineral dissolution and secondary mineral precipitation than the lower part of the profile that resides in the saturated zone. This results in lower saturated hydraulic conductivities in the near surface than at depth.
At the hillslope scale, it is also observed that neoformed fine particles accumulate in lowlying areas of the landscape. While this can result from particle translocation by physical erosion, it can also result from the accumulation of solutes along water flow paths that result in supersaturation and precipitative deposition. The latter is observed in the model simulations that predict an accumulation of secondary minerals in the toe-slope positions (Figure 7). Areas of concentrated flow, below the water table, became equilibrium saturated with respect to primary minerals undergoing dissolution and, as a result, showed less dissolution and precipitation.
Dissolution in this portion of the hillslope can be characterized as transport limited while secondary phase accumulations were apparently limited by precipitation kinetics and low oxygen concentrations (inhibiting Fe(II) oxidation). The precipitation band that is observed to form on the border between the saturated and unsaturated zones is particularly interesting, and is due to the intersection of geochemical conditions that promoted the precipitation of smectite (above) and ferrous hydroxide (below).
While not captured in the current simulations, newly formed heterogeneity in hydraulic conductivity is expected to feedback to influence the further the formation of soil structure, including preferential flow paths, thereby resulting in additional changes to water transport, solute or sediment flux and mineral weathering (Uchida et al., 2005). For example, a decrease in K sat in the near surface may diminish infiltration relative to runoff generation during intense rainfall events, thereby promoting overland flow via the Horton runoff generation mechanism (Maxwell and Kollet, 2008). However, it should also be pointed out that this study assumes bulk changes to the particle size of the soil resulting from mineral transformation processes. While this may be the case to some extent, it is also expected that neoformed clay particles will be concentrated in close association with existing particles and in small pores, thereby possibly having less influence on large pores and saturated hydraulic conductivity values (relative to small pores and unsaturated hydraulic conductivity values). Precipitation of secondary clays in small pores will also promote aggregate formation, with impacts on broadly relevant processes in the B2 experiment such as carbon sequestration and evapotranspiration.
In the present study, results of CrunchFlow2007 simulations were transferred into the Rosetta code to calculate hydrologic parameters based on change in particle size distribution of the domain. To enable future model coupling, this hydrologic parameterization would need to be translated back into HYDRUS for a redefinition of the hillslope hydrologic properties. Since each such translation between models (three in this case) requires multiple participants and the potential associated introduction of error, the B2 hillslope experiment has the goal of developing a fully-coupled hydrologic and geochemical model that will be capable of evaluating feedbacks between geochemical reactions and hydrologic response in single model runs. Ideally, a fully coupled model would also describe the interaction between climatic, hydrologic, geochemical, and biological processes in the evolution of fluid flow paths. This paper and the companion papers dealing with hydrologic and vegetation modeling of the B2 hillslopes (Hopp et al., 2009) are a first approach toward this goal.  (Wolery et al., 1990), with Gibbs energy contribution from entropy calculated assuming ideal mixing b (Lelli et al., 2008); c EQ3/EQ6 database (Wolery et al., 1990); d (Yang and Steefel, 2008), rate recalculated from T = 22 °C; e (Carroll and Knauss, 2005); f (Berner et al., 1980); g (Grandstaff, 1980); h (Gislason and Oelkers, 2003), rate recalculated from T = 7 °C; i (Ganor et al., 1995); j (Cama et al., 2000), rate recalculated from T = 80 °C; k (Ganor et al., 1999); l (Cornell et al., 1975);m (van Hees et al., 2002); m (Hellmann and Tisserand, 2006), rate recalculated from T = 150 °C; o (Schott et al., 1981) ; p (Wogelius and Walther, 1991) .