Coupled modeling of hydrologic and geochemical fluxes for prediction of solid phase evolution in the Biosphere 2 hillslope experiment

Coupled modeling of hydrologic and geochemical fluxes for prediction of solid phase evolution in the Biosphere 2 hillslope experiment K. Dontsova, C. I. Steefel, S. Desilets, A. Thompson, and J. Chorover B2 Earthscience, The University of Arizona, Tucson, AZ, USA Department of Soil, Water & Environmental Science, The University of Arizona, Tucson, USA Earth Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA, USA Department of Hydrology and Water Resources, The University of Arizona, Tucson, AZ, USA Department of Crop and Soil Sciences, University of Georgia, Athens, GA, USA Received: 24 May 2009 – Accepted: 26 May 2009 – Published: 18 June 2009 Correspondence to: K. Dontsova (dontsova@email.arizona.edu) Published by Copernicus Publications on behalf of the European Geosciences Union.


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, firstorder 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)  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 Biosphere 2. 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.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 onedimensional simulations evaluating effects of mineralogy, particle size (surface area), Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion 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 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 (Fig. 1).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 Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc

Printer-friendly Version
Interactive Discussion 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).
The latter material was composed by sieving and remixing different size fractions of the original basalt to give a larger fraction of fine particles.Particle size of the basaltic material as received was 97.58% sand, 2.20% silt, and 0.22% clay, whereas that of the reconstituted loamy sand was 77.74% sand, 21.04% silt and 1.23% clay.Specific surface areas were 1.438 m 2 g −1 as received and 3.238 m 2 g −1 for loamy sand.Measured particle density was 2.81 g cm −3 .
Total elemental composition of basalt is characterized by the following mass percent of oxides: 47.82% 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 (Sverdrup and Warfvinge, 1995).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 298 K were scaled using the Arrhenius relation (Lasaga, 1984): where E a is activation energy, k m,298 is the rate constant at 298 K, 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 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).Primary aqueous species (H + , Na + , Mg 2+ , Ca 2+ , NO − 3 , SO 2− 4 , Fe 2+ , K + , Al 3+ , SiO 2(aq) , HCO − 3 , O 2(aq) , Cl − , Ti(OH) 4(aq) , HPO 2− 4 , 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 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 1-D simulations, to provide a comparison, whereas diffusion-controlled gas flux was used for 2-D 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 1-D 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 2-D 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 1-D simulations and 30 1-m segments in the Introduction

Conclusions References
Tables Figures

Back Close
Full Dissolution rates in CrunchFlow2007 (Steefel, 2008) are described using: where , Q is the ion activity product for the mineral-water reaction, K SO is the corresponding solubility product constant, and a n i 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).
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 rate at conditions closer to equilibrium; m 1 =1.17, m 2 =0.000078, and m 3 =3.81.For other minerals where dual rate constants are not available, a single rate constant (with m 1 =0.5, m 2 =1, and m 3 =2) was employed.
Equation 3 describes transition state theory (TST) if parameters m 1 , m 2 , and m 3 are all equal to 1 and there is no inhibition or catalysis of the reaction by ions in solution.Linear transition state theory has been used to describe precipitation of secondary minerals, although there is some evidence that kaolinite dissolution may be governed Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc

Printer-friendly Version
Interactive Discussion 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).

One-dimensional simulations
Initial simulations were of a 1-D, 30 m length soil column (to model maximum hillslope flow path length).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.238 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 Brantley, 2007).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 2-D 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).Climate data were taken from (i) the Lucky Hills site within the USDA-ARS Walnut Gulch Experimental Watershed, Arizona, USA (Scott et al., 2000), and (ii) a scenario representative of a semi-arid "Sky Island" Ponderosa Pine forest (Brown-Mitic et al., 2007), a subset of the US western sub-alpine forest (a semi-annual total rainfall of 367 mm distributed over 90 days, followed by a dry period of 90 days).For Lucky Hills conditions (total of 356 mm yr −1 , 0.2655 mm h −1 average daily rainfall rate for 56 days per year) steady state flow conditions were reached after 1044 h (43.5 days) of simulated rainfall.For simulated Sky Island forest climate (734.2 mm yr −1 , 0.0041 m h −1 average daily rainfall over 74 days per year) steady state was reached after 980 h (40.8 days).These HYDRUS-3D-derived steady-state saturation and velocity profiles were used to run 2-D CrunchFlow2007 simulations for 0, 168, 560, and 1000 days for the "Lucky Hills climate" (simulating 0, 3, 10, and 18 years of the actual hillslope experiment under Lucky Hills climate conditions) and for 0, 222, 740, and 1000 days for the "Sky Island climate" (simulating 0, 3, 10, and 14 years of the actual hillslope experiment under Sky Island climate conditions).For CrunchFlow2007, only days with rain were modeled and days without rain were omitted.Therefore, since the Sky Island climate had more days with rain per year (74 vs. 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 Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion 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 (2001).Rosetta is a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions.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.

One dimensional CrunchFlow2007 simulations
1-D simulations indicate primary mineral dissolution and secondary mineral precipitation in the porous basalt column representing a 30-m flowpath on the hillslope (Fig. 3).
Basalt glass shows the most rapid dissolution, followed by forsterite.In contrast, anorthite and diopside did not undergo significant dissolution during the time-scale of modeling.Irrespective of whether O 2 was limited by diffusion into the surface of the column, or whether it was present in excess throughout the column, smectite and goethite were the predominant newly-precipitating phases in the top (30 cm) of the column (Fig. 3).
Kaolinite and pyrolusite were also formed, but to a lesser extent.Farther along the Introduction

Conclusions References
Tables Figures

Back Close
Full column length, at 45 cm and deeper, under O 2 -limited conditions, a small amount of Fe(OH) 2 precipitation was observed.The total mineral transformation was slightly greater in the presence of excess O 2 , indicating slight limitations on weathering imposed by oxidant availability.
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 (Fig. 4).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 2-D 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 1-D simulations, 2-D secondary mineral phase distribution and Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion their respective saturation indices (log Q/K SO ) show smectite as the principal neophase formed in both climatic regimes, particularly in the vadose zone (Fig. 5).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 Fig. 6, we have plotted the volumetric depletion of primary minerals (left) and accretion of secondary minerals (right) following a steady-state water flux equivalent to that water flux deriving from (A) 3, (B) 10, and (C) 18 years of weathering in the Lucky Hills climate, or (A) 3, (B) 10, and (C) 14 years of weathering in the Sky Island climate.In both left and right sides of Fig. 6, large and small changes are indicated by hot and cool colors, respectively.The simulations show significant differentiation of the hillslope profile over the time-scale of the B2 experiment, largely in response to the non-uniformity in flow.Substantial primary mineral dissolution is observed in the surface layers of the porous basalt medium that are exposed to fresh infiltration water.Primary mineral dissolution is diminished deeper in the domain where pore-waters are enriched in aqueous phase reaction products, and the solution is closer to equilibrium saturation with respect to basalt-derived minerals.Distinct stratification in mineral precipitation patterns is also consistent with hillslope flow patterns (Fig. 2).The volumetric sum of secondary mineral precipitation was greatest near the hillslope surface where slow, unsaturated flow was predominant.In addition, a Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion band of precipitation is observed between saturated and unsaturated domains, where vertical flow from the surface intersects the water table (which is closer to the surface under Sky Island climate).The patterns for both dissolution and precipitation followed the flow patterns for both climates (Fig. 2), with greater mineral transformation in the Lucky Hills climate due to a vadose zone extent wherein most of the weathering occured (i.e.13.2% dissolution, 8.98% precipitation in Lucky Hills climate versus 10.5% dissolution, 7.02% precipitation in Sky Island climate).Loss of lithogenic elements from the hillslope (Fig. 7) showed a similar pattern, reflecting greater weathering in the Lucky Hills climate.

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 (Fig. 8).Throughout the hillslope, K sat was decreased from 0.044 m h −1 G, 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 Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion 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 (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 Brantley, 2007).Another proposed reason for discrepancies between laboratory-scale and fieldscale 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 2-D simulations with the understanding that predicted mineral transformations may be higher than those observed in the experiment.Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version

Interactive Discussion
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 time-scale 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 1-D simulations indicated that weathering kinetics were controlled by mineral dissolution rate under the modeled conditions.
Conversely, the 2-D 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 ad-Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion vance (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 1-D and 2-D 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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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 low-lying 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 toeslope positions (Fig. 6).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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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)  Full "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 physico-chemical parameters on mineral transformation rates in 1-D basalt soil columns, and (ii) to assess the influence on mineral weathering of 2-D fluid flow in the vadose and saturated zones, including feedbacks to hydrologic flow.The 2-D simulation represents the central axis of the first order hillslope studied byHopp et al. (2009).
direction and 11 0.091-m segments in surface to bottom direction in 2-D simulations.