Interactive comment on “ Modeling the impact of soil aggregate size on selenium immobilization ” by

1. p12054, ln 8-9: You state that the individual soil aggregates were placed in the center of the reaction cell and then a constant flow rate was applied. Was anything else packed in the cell or was the soil aggregate allowed to move around during the experiment? Based on the larger cell dimensions and the spherical nature of the aggregate, I would imagine that the aggregate would not remain stationary while under flow conditions. Do you think this could affect your modeling results? C8088


Introduction
Selenium is a trace element of great ecological significance due to its dual role as both an essential nutrient and environmental contaminant.The intake range between deficiency and toxicity of selenium for most animals is exceedingly small (Lenz and Lens, 2009).For example, humans need around 40 µg day −1 for optimal health; however, an intake of more than 400 µg day −1 can cause selenosis (WHO, 1996).Selenium is heterogeneously distributed across terrestrial landscapes, with seleniferous (> 0.5 mg kg −1 ) and selenium-deficient soils (< 0.1 mg kg −1 ) sometimes occurring as close as 20 km from one another (Lenz and Lens, 2009).The irrigation of seleniferous soils can lead to the production of sub-surface drainage water containing high selenium concentrations that cause ecological damage in downstream aquatic ecosystems.This process was first recognized in the late 1980s after the occurrence of numerous deaths and embryonic deformities in waterfowl populations near the agricultural evaporation ponds of the Kesterson reservoir, a former wildlife refuge in California (Presser, 1994).In the Western United States alone, nearly 400 000 km 2 may be at risk of irrigation-induced selenium contamination (Seiler et al., 1999) and similar problems have also been observed in Canada, Egypt, Israel, and Mexico (Lemly, 2004).Enhanced understanding of selenium transport and retention in surface soils will improve the management of seleniferous soils to diminish the risk of contamination for downstream ecosystems.
The known microorganisms capable of selenium reduction are scattered across a very diverse set of prokaryotic taxa (including Crenarchaeota (Domain: Archaea), low-and high-GC gram-positive bacteria, β-, γ -, and ε-Proteobacteria) (Stolz et al., 2006).We chose Enterobacter cloacae SLD1a-1 as the benchmark organism for the reactive transport model presented in this paper.Since its original isolation from seleniferous agricultural drainage water, this γ -Proteobacterium has received a particularly large amount of research attention (Losi and Frankenberger, 1998).E. cloacae SLD1a-1 is known to grow both aerobically and anaerobically using a variety of electron donors, including pyruvate (Losi and Frankenberger, 1997b).It can reduce selenate via selenite to solid elemental Se (Losi and Frankenberger, 1997a) and its membrane-bound selenate reductase has been identified (Watts et al., 2003) and purified (Ridley et al., 2006).Furthermore, it has been used as a model organism in several studies investigating the molecular genetics of selenium reduction (Ma et al., 2007(Ma et al., , 2009;;Yee and Kobayashi, 2008;Yee et al., 2007) and has been considered as a candidate for bioremediation schemes (Losi and Frankenberger, 1998).E. cloacae SLD1a-1 can carry out selenium reduction under oxic conditions, though at rates that are an order of magnitude lower than under anoxic conditions (Losi and Frankenberger, 1997b).This inhibition is in accordance with thermodynamic expectations, since at neutral pH reduction of selenate to selenite is not energetically favored above a pE of 7.5, placing it generally between nitrate and manganese oxide on the redox ladder (Sposito et al., 1991).Unsurpris-ingly, oxygen inhibition of microbial selenium reduction has also been observed in environmental samples (Oremland et al., 1989).Selenium reduction rates are thus sensitive to local redox potential.
Redox conditions in soils are known to display a large amount of spatial heterogeneity, owing to the heterogeneous size and distribution of pore spaces in a soil's physical structure (Arah and Vinten, 1995;Sexstone et al., 1985).A common approach to model this spatial heterogeneity in mechanistic biogeochemical models has been to represent the soil matrix as a series of spherical microporous structures or aggregates (Arah and Vinten, 1995).Soil aggregates are mmto cm-sized structural units of mineral particles bound by roots, fungal hyphae, and organic matter that naturally occur in most soils (Brady and Weil, 2002).Whereas a variety of different shapes exist in nature, granular spherical aggregates (granular structure) are particularly common in surface soils (A horizons) (Brady and Weil, 2002).Aggregates represent well-defined natural systems in which to study the impact of mm-to cm-scale spatial heterogeneity on biogeochemical reactions (Tokunaga et al., 2003).
Dissolved chemical concentrations within aggregates depend on the local coupling of reactions and transport.Whereas fast advective solute transport is prevalent in the inter-aggregate macropores, transport inside aggregates is mostly diffusive due to the low permeability of the microporous structures.In conjunction with local microbial metabolic activity this often leads to the formation of steep chemical gradients within aggregates (Arah and Vinten, 1995;Sexstone et al., 1985;Kaurichev and Tararina, 1972), as the consumption (or production) of a chemical species inside an aggregate outpaces the diffusive supply from (or export to) the surrounding macropore (Tokunaga et al., 2003).Full transitions from oxic to anoxic have been observed in aggregates as small as 4 mm in diameter (Sexstone et al., 1985).Tokunaga et al. (1994) demonstrated that such anoxic zones can influence the mm-scale distribution and mobility of selenium based on synthetic, flat-aggregate no-flow systems.We previously developed a system of saturated-flow reactors containing single idealized soil aggregates, for the investigation of aggregate scale gradients in the microbial reduction and secondary mineral formation of iron (Pallud et al., 2010b) and a mathematical model describing the coupled transport and reactions of iron within the system (Pallud et al., 2010a).This unique experimental approach, emulating field-like transport conditions at the aggregate scale while allowing for controlled experiments, was later applied to investigations of arsenic mobility (Masue-Slowey et al., 2011) and heterogeneity in microbial selenium reduction (Kausch et al., 2012).In the latter study, we found increasing concentrations of reduced selenium towards the core of aggregates containing E. cloacae SLD1a-1 under both oxic and anoxic conditions independent of selenate or carbon concentrations in solution.This was hypothesized to be due to limitations in the transport rate of reduced selenium products out of the aggregate, which may imply a general role of aggregate size in selenium reduction and retention (Kausch et al., 2012).
In this study we formally tested this hypothesis by constructing a mechanistic mathematical model of the coupled transport and transformations of selenium species inside an aggregate containing selenium-reducing bacteria (E.cloacae) and in the surrounding solution.The objectives were to evaluate the theoretical impact of variations in aggregate size on selenium retention within a comprehensive mechanistic framework and to assess implications for the management of seleniferous soils that can be a source of water contamination.Furthermore, we wanted to disentangle the roles of heterogeneity in reaction rates and chemical fluxes in producing the experimentally observed concentration gradients and resolve the production and distribution of elemental selenium in aggregates (which could not be directly measured in experiments).To this end we reproduced our experimental results (Kausch et al., 2012) to validate the model, and then simulated selenate and selenite reduction in aggregates ranging from 1-2.5 cm in diameter.

Overview of experimental data used
In developing and validating our reactive transport model we used data obtained from flow-through experiments involving novel artificial soil aggregates made of sand and containing Enterobacter cloacae SLD1a-1 (Kausch et al., 2012).The experiments were carried out under a set of eight different chemical conditions (0.4 or 0.8 mM selenate and 0.3 or 1.2 mM pyruvate in the input solution under oxic and anoxic conditions).We used data on bacterial cell densities and concentrations of reduced selenium within three concentric sections of the aggregates, as well as data on temporal evolution of selenite and flow tracer (bromide) concentrations in the reactor outflow.The flow tracer data were used to establish that the physical flow model was representative of experimental flow and transport.The data on solution and solid-phase selenium concentrations were used to parameterize the global expressions for chemical reaction rates and later to assess that the reactive transport model represented the experimentally observed spatial and temporal dynamics of the system under the full set of chemical conditions.The data on cell density were used to constrain cell density in the reactive transport model by checking whether it varied between concentric sections of the aggregates and between the beginning and end of experiments.

Aggregate construction
Spherical (2.5 cm diameter) aggregates were constructed in the laboratory using the protocol developed by Pallud et al. (2010b).Aggregates were composed of pure quartz sand (IOTA 4 pure quartz powder, grain size of 150-250 µm, Unimin Corporation) homogenously inoculated with a suspension of E. cloacae to a density of approximately 10 8 cells g −1 dry aggregate mass.Right after the sand was inoculated with E. cloacae, it was mixed with hydrogel agarose to promote particle aggregation (150 g of 0.5 wt % agarose, 30 mL of bacterial suspension and 300 g of sand) and then molded into aggregates using a spherical press.The average pore diameter and porosity of aggregates thus constructed are 39 µm and 0.58, respectively (Pallud et al., 2010b).No carbon source or selenium was present in the aggregates at the time of synthesis.Pure quartz sand was chosen since adsorption of selenium oxyanions is negligible under the experimental conditions, thus allowing the investigation of transport mechanisms independent of sorption.

Flow-through experiments
Flow-through systems consisted of individual artificial soil aggregates placed on hollow plastic stands, which supported them in a stable position without inhibiting flow, at the center of a polycarbonate reactor cell (5.1 cm diameter, length 3.7 cm).A constant flow rate of 1.0 ± 0.1 mL h −1 was imposed at the lower boundary of the reactors using a peristaltic pump and effluent samples were collected in 8 h increments for 192 h.The input solution was a sterile, artificial groundwater medium (prepared in 10 mM PIPES and containing the following (in mg L −1 ): NaCl, 30; NH 4 Cl, 0.95; KCl, 5; MgSO 4 , 50; KH 2 PO 4 , 0.950 and 1 mL L −1 mineral Wolfe's solution with pH adjusted to 7.2), with selenate concentrations of either 0.4 or 0.8 mM, and pyruvate concentrations of either 0.3 or 1.2 mM (added from filtered sterile stock solutions).Bromide (2 mM, as KBr) was included in inflow solutions as an unreactive flow tracer, but not in the solution used to fill the reactors initially.For each experiment the composition of the input solution remained constant throughout.
Anoxic experiments were fully contained within an anaerobic chamber (Coy Laboratories, Inc., Grass Lake, MI), and the input solutions (single batches of 0.5 L) were made anoxic by sparging with O 2 -free N 2 for 15 min.The efficacy of the chosen sparging procedure and time was ascertained through tests with the redox indicator die resazurin (0.5 mg L −1 ).For oxic experiments, thorough aeration was insured by using oxic input solutions and continuously bubbling air into the solution surrounding the aggregates throughout the experiment.

Chemical analysis
Selenium extraction from the solid phase was performed with 12 M HNO 3 (Kausch et al., 2012) and analysis of solid extraction and elution samples was performed via hydride generation coupled to optical emission spectrometry (PerkinElmer 5300 DV) as previously described by Kausch et al. (2012).While the procedure allows for direct analysis of selenite concentrations in liquid samples, the extraction and analysis of solid-phase samples yields a single experimental fraction of reduced selenium including both selenite and elemental selenium (Kausch et al., 2012).Outflow samples were analyzed for selenite and total selenium, with selenate concentrations computed as the difference between the two measurements.Selenate outflow concentrations did not differ significantly from input concentrations and are thus not reported.Additionally, for outflow samples from two oxic and two anoxic experiments each, bromide concentrations were determined on a Dionex ion chromatograph (ED40 electrochemical detector) with an AS23 column (Thermo Fisher Scientific Inc.) using 4.5 mM Na 2 CO 3 , 0.8 mM NaHCO 3 eluent.Bromide and selenite standards were prepared from certified stock (Spex CentriPrep Group L.L.C. and VHG Labs, respectively).

Bacterial cell counts
Counts of bacteria per gram of sand were obtained prior to and at the conclusion of experiments.Prior to experiments, the sand-agarose-cell mixture from which aggregates were shaped was sampled, whereas at the conclusion of experiments separate samples from the exterior and mid-sections of each individual aggregate were used.The volume of sand in core sections was too small to allow bacterial counts for that aggregate section.The 1 g (wet weight) samples were each transferred into 9 mL of sterile saline solution (NaCl, 8 g L −1 ) and shaken thoroughly.A 10-fold dilution series in sterile NaCl (8 g L −1 ) was prepared from each aggregate sample solution.Then 100 µL aliquots of appropriate dilutions were plated in triplicate on solid medium (30 g L −1 tryptic soy broth, 15 g L −1 agar).Plates were incubated for 2-5 days at 25 • C before colony forming units (CFU) were counted.

Reactive transport model
The transport model and geometry used here are based on previous work modeling iron reduction in artificial aggregates (Pallud et al., 2010a).The experimental setup was emulated in a 2-D axi-symmetric reactive transport model (Fig. 2a) describing physical solute transport by advection and molecular diffusion and the effect of kinetically driven redox reactions on mass balance expressions for pyruvate, oxygen, selenate, selenite, and elemental selenium.
The flow and pressure fields were computed using the stationary, incompressible Navier-Stokes equations in the free fluid surrounding the aggregate and the Brinkman equation inside the aggregate (see Appendix A for an explicit account of the governing equations).The aggregate permeability was computed based on grain size and porosity measurements and a uniform flow velocity was imposed at the lower boundary of the reactor (inflow), with no-slip conditions at the reactor walls (Pallud et al., 2010a).
Reactions of selenium species, pyruvate, and (in oxic simulations) oxygen were formulated by taking into account aerobic respiration and selenium reduction coupled with pyruvate oxidation by E. cloacae (Losi and Frankenberger, 1998).All reactions were implemented using kinetic rate expressions with a Monod-type dependency on both electron acceptor (i = oxygen, selenate, or selenite) and electron donor (pyruvate) (Rittman and VanBriesen, 1996) as where R i is the reaction rate and k i the cell-specific rate constant for the respective electron acceptor (i), ρ cells denotes the cell density of E. cloacae, and K pyruvate mi and K i m are the half-saturation constants for electron donor and acceptor, respectively.I i ([oxygen]) is the oxygen concentration dependent inhibition term (Van Capellen and Gaillard, 1996) specific for selenate or selenite reduction (I i ([oxygen]) = 1 for aerobic respiration).I i ([oxygen]) was parameterized by fitting the data collected by Losi and Frankenberger (1997b) on the inhibition of selenate and selenite reduction by E. cloacae at various oxygen concentrations (Supplement: Fig. S1).
A cell death term was introduced that led to a linear decrease in cell density over time (see Appendix A for details) based on our observations that the average cell density measured at the conclusion of anoxic experiments was an order of magnitude lower than the initial cell density (see below).This term was only used in fitting the model to experimental results and not in the predictive simulations.Other than those relating to cell density, model parameters were kept constant among all simulations (i.e., different chemical conditions).For the predictive simulations, the cell density was also held constant at 1 × 10 8 cells g −1 for all chemical conditions.The values of reaction parameters used in the model are listed in Table 1.A detailed description of the model is given in Appendix A.

Model fitting
For all eight chemical scenarios in parallel, the reactive transport model was fitted to both temporal data of selenite breakthrough curves and spatial data of reduced selenium solid-phase concentrations, by manually adjusting reaction rate and cell density parameters.Among the reaction rate parameters (Table 1), only the maximum cell-specific rate constants, the electron donor half-saturation constant for c Oxygen inhibition was modeled based on best-fit to data from Losi and Frankenberger (1997b), with an inverse-Monod relationship for selenate reduction (half-inhibition constant in units of mM) and an exponential inhibition term for selenite reduction (unitless inhibition constant).See Appendix A for details.selenate reduction, and the oxygen half-saturation constant were adjusted; all other parameters were taken from the literature.Whereas reaction rate parameters were kept constant between chemical scenarios, the initial cell density values were varied between 1-3 × 10 8 cells g −1 .See Appendix A for a discussion of the cell death term.

Predictive simulations
Predictive simulations were performed by varying the aggregate diameter for chemical scenarios equivalent to those run experimentally (0.3 and 1.2 mM pyruvate, with 0.4 and 0.8 mM selenate input solutions under both oxic and anoxic conditions).Each predictive scenario consisted of a set of 9 simulations in which the aggregate diameter varied between 1 cm (which corresponds to the size of core sections in experiments) and 2.5 cm (which corresponds to the full size of experimental aggregates).Results were linearly interpolated between the 9 simulated diameters.Simulation time (192 h) and kinetic parameters (Table 1) were left unchanged and the cell density was constant at 1 × 10 8 cells g −1 .3 Results

Experimental selenite breakthrough curves
Experimental selenite breakthrough curves from oxic and anoxic experiments are shown in Fig. 1a and b, respectively.It can be seen that under both oxic and anoxic conditions, selenite outflow concentrations reached their highest value between 80 and 150 h.Selenite concentrations in the reactor outflow were 6-20 times higher under anoxic than oxic conditions, with peak selenite concentrations of respectively 40-70 µM and 2-12 µM.The initially steeply increasing selenite concentrations plateaued as they approached their peak value; however, a clear steady state was not maintained under all experimental conditions.Selenite concentrations decreased after reaching the peak value in the oxic experiment with high-pyruvate (1.2 mM) and selenate (0.8 mM) concentrations in the input solution and in all anoxic experiments, except the one with low-pyruvate (0. the outflow under both oxic and anoxic conditions (Fig. 1).Under oxic conditions, peak selenite concentrations approximately tripled with an increase of the pyruvate input concentration from 0.3 mM to 1.2 mM, whereas they approximately doubled with an increase of selenate input concentrations from 0.4 mM to 0.8 mM (Fig. 1a).In contrast, under anoxic conditions peak selenite concentrations increased by a factor of less than 1.5 with both an increase of the pyruvate input concentration from 0.3 mM to 1.2 mM, and an increase of selenate input concentrations from 0.4 mM to 0.8 mM (Fig. 1b).

Experimental solid-phase distribution of reduced selenium
Reduced selenium, which corresponds to the sum of selenite and elemental selenium, was detected in core, mid-, and exterior sections of all aggregates, with positive gradients in reduced selenium concentrations observed from the aggregate exterior towards the cores of aggregates under all investigated chemical conditions (Fig. 2).These gradients were steeper under oxic than under anoxic conditions.Specifically, the ratio in solid-phase concentrations of reduced selenium between core and exterior sections of aggregates was between 2.2-3.9 and 1.8-2.0under oxic and anoxic conditions, respectively.Absolute concentration values of solid-phase reduced selenium varied between 0.01-0.05µmol g −1 and 0.16-0.34µmol g −1 under oxic and anoxic conditions, respectively (Fig. 2).The differences within the oxic and anoxic sets of solid-phase concentrations were not statistically significant (Kausch et al., 2012), but are here presented separately for the purpose of showing the fits of individual model scenarios.For example, under oxic conditions the concentrations of reduced solid-phase selenium appear systematically higher for pyruvate input concentrations of 1.2 mM as compared to pyruvate input concentrations of 0.3 mM, and while this difference in experimental results was not large enough to be statistically significant it was nevertheless reproduced by the simulated results.

Flow tracer (bromide) breakthrough curves and transport model
Steady-state concentrations of bromide (within 10 % of the 2 mM input concentration) in the outflow from the four reactors tested were reached between 80 and 120 h (Supplement: Fig. S2).This compares well with results from the transport model, which predicts that 90 % of the input concentration of an unreactive tracer should be reached after 116 h in the outflow.Experimental data points scatter around the simulated curve with a maximum distance of 0.5 mM (Fig. S1).

E. cloacae cell counts prior to and at termination of experiments
The density of E. cloacae (average ± standard deviation) determined in aggregate material prior to experiments was (1.2 ± 0.5) × 10 8 CFU g −1 (N = 3).Under anoxic conditions, cell densities decreased by an order of magnitude over the course of the experiments.In fact, the average value from replicates from all aggregate sections measured at the experiment conclusion was (8 ± 7) × 10 6 CFU g −1 (N = 25).However, under oxic conditions, cell densities remained relatively constant over time with an average cell density from replicates of all aggregate sections at the end of experiments of (1.1 ± 0.6) × 10 8 CFU g −1 (N = 26), one outlier aggregate section excluded.The only exception was the exterior section of the oxic aggregate run with high-pyruvate (1.2 mM) and high-selenate (0.8 mM) input concentrations, where E. cloacae density reached a value of (1.4 ± 0.9) × 10 9 CFU g −1 at the conclusion of experiments (N = 3), which is an order of magnitude higher than the initial value.Aside from this no systematic differences in the cell density based on aggregate section or input solution composition were discerned.

Model fitting and validation
The temporal and spatial concentration data produced by the reactive transport model are compared to experimental observations in Figs. 1 and 2, respectively.The model accurately fitted both the breakthrough curves and solid-phase data with a single kinetic parameterization (Table 1) for the full set of experiments.Similarly to what was observed in experiments, simulations showed quasi-steady state in selenite outflow concentrations being reached after around 100 h regardless of the chemical conditions (Fig. 1).The decrease in selenite outflow concentrations after 100 h in three of the anoxic and one of the oxic experiments was also reproduced as a result of a decrease in modeled cell densities (active cell death term).The reactive transport model also accurately reproduced the differences in outflow selenite concentrations observed experimentally for different chemical conditions.In fact, peak selenite concentrations in the modeled outflow ranged between 2-11 µM and 40-70 µM for oxic and anoxic conditions, respectively, compared to experimental values of 2-12 µM and 40-70 µM.The experimentally observed spatial trend of increasing concentrations of reduced selenium towards the core of aggregates was also reproduced by the model (Fig. 2).Comparing modeled and experimentally determined solid-phase concentrations of reduced selenium individually based on aggregate section and chemical conditions (Fig. 2), it can be seen that the modeled results lie within one standard deviation of the experimental measurements in 8 out of 12 and 5 out of 12 sections for oxic and anoxic conditions, respectively.Furthermore, it is obvious that while absolute concentrations were not always reproduced, the spatial trend was reproduced for all chemical scenarios.Specifically, the modeled concentration ratio between the simulated core and exterior sections was 1.8 for anoxic conditions, and 2.8-4.9 for oxic conditions, compared to experimental values of 1.8-2.0 and 2.2-3.9.

Dynamics of concentration fields and reaction rates
The spatial and temporal behavior displayed by chemical species in our reactive transport model differed qualitatively based on whether they were reactants or products of reactions, whether they were solids or solutes, and also based on the bulk oxidation conditions.Here we use concentration and reaction rate maps from the oxic simulation with 0.3 mM pyruvate and 0.4 mM selenate as an example to illustrate spatial and temporal patterns that were consistent among the different chemical scenarios explored.Figure 3 shows the final concentration fields for pyruvate, oxygen, selenate, selenite, and elemental selenium at the termination of this simulation.Movies showing the dynamics of pyruvate, oxygen, selenate, selenite, and elemental selenium concentration fields (Movies S1-S5) and of reaction rate fields (Movies S6-S7) are included in the Supplement.
Whereas the concentrations of pyruvate, selenate, and selenite approached steady state between 64 and 192 h of simulation time (Movies S1, S3, S4), those of oxygen reach steady state within 32 h (Movie S2).Concentrations of elemental selenium, on the other hand, kept increasing inside the aggregate domain throughout the simulation time (Movie S5).Over the course of simulations, all compounds developed smooth gradients of approximately radial symmetry within the aggregate domain (the symmetry is slightly shifted in the direction of flow).However, whereas concentrations of re- actants (pyruvate, selenate, and when present, oxygen) decreased towards the core of the aggregate domain, those of products (selenite and elemental selenium) increased towards the core (Fig. 3).The gradients of selenite, elemental selenium and, except in the case of selenate, also of reactants, were steeper under oxic than under anoxic conditions (data not shown).The maximum concentrations of the reactants pyruvate, selenate, and oxygen found in the free fluid surrounding the aggregate at steady state correspond to the input concentrations (Fig. 3, note: 0.26 mM correspond to atmospheric equilibrium oxygen concentration).
Diffusive and advective transport, respectively, dominated inside the aggregate and in the surrounding solution.This is illustrated for the case of selenite in Fig. 4. It can be seen that while the diffusive selenite flux was near zero at the core of the aggregate, it increased towards the aggregate exterior and is 0.4 nmol m −2 s −1 at the interface between the aggregate and the surrounding solution (Fig. 4a).In the surrounding solution the diffusive flux decreased with increasing distance from the aggregate.This was, however, compensated by an increasing advective flux component downstream of the aggregate, where it lead to a maximal total selenite flux of 0.5 nmol m −2 s −1 (Fig. 4b).
Similar to solute concentrations, rates of selenate and selenite reduction reached a steady state after 60 and 150 h of simulation time, respectively (Movies S6-S7).The rates of selenate reduction approached a steady state that is relatively homogenous throughout the aggregate (Fig. 4c).In contrast, the rates of selenite reduction approached a steady state that mirrored the concentration field of selenite (with increasing values towards the core of the aggregate domain, Fig. 4d).In the example simulation the steady-state rates of selenate reduction were 130-190 nM s −1 everywhere in the aggregate domain (Fig. 4c), whereas those of selenite reduction smoothly increased from 0 at the aggregate exterior to 12 nM s −1 at the core (Fig. 4d).

Predictive simulations of aggregate size impact on selenium retention
Simulations in which the diameter of the spherical aggregate was varied between 1 and 2.5 cm (corresponding to volumes of 0.5 cm 3 and 8.2 cm 3 , respectively) showed that the aggregate solid-phase concentrations of reduced selenium scale with aggregate size under all chemical conditions investigated (Fig. 5).This effect was more pronounced under oxic conditions than under anoxic conditions.In fact, under oxic conditions the average solid-phase concentrations of total reduced selenium were 0.6-1.6 nmol g −1 in aggregates of 1 cm diameter and 4-17 nmol g −1 in aggregates of 2.5 cm diameter, which corresponds to a difference by a factor of 6-11 between the two size extremes (Fig. 5a-d).On the other hand, under anoxic conditions the average solid-phase concentrations of total reduced selenium were only about 4 times higher in 2.5 cm diameter aggregates than in 1 cm diameter aggregates, with values of 57-140 nmol g −1 and 14-36 nmol g −1 , respectively (Fig. 5e-h).Under oxic conditions, solid-phase concentrations of elemental selenium in aggregates were particularly sensitive to aggregate size, scaling from 0.1-0.2 to 0.8-4.5 nmol g −1 over the investigated size range (factor 8-23).Thus, the fraction of elemental selenium in aggregates varied greatly with aggregate size under oxic conditions (11-26 % of reduced selenium).In contrast, under anoxic conditions elemental selenium makes up around 60 % (55-61 %) of reduced selenium in aggregates regardless of aggregate size and input solution concentrations.
Under anoxic conditions, the curves describing reduced selenium solid-phase concentrations as a function of aggregate size displayed negative curvature throughout the different reactant concentration scenarios investigated (Fig. 5e-h).However, under oxic conditions the curves corresponding to high-pyruvate (1.2 mM) concentration scenarios (Fig. 5c-d) were markedly steeper towards the larger end of the aggregate size spectrum (positive curvature) than those for lowpyruvate (0.3 mM) concentration scenarios (negative curvature) (Fig. 5a-b).

Model validity
The reactive transport model developed in this work reproduces the general spatial and dynamic behavior observed in experiments across the full spectrum of chemical conditions investigated.More specifically, our model reproduces the time to reach quasi-steady state, the impact of aeration conditions and input solution composition on outflow selenite concentrations (which are a proxy for experimental reac-  and (d and h) 1.2 mM pyruvate/0.8 mM selenate.The amount of selenite in the solid phase is shown in grey, while the amount of elemental selenium in the solid phase is shown in red.Each panel shows the results of 9 simulations carried out at different aggregates sizes; the space between simulations was linearly interpolated.tion rates; Kausch et al., 2012), and the spatial distribution of reduced selenium inside the aggregates (i.e., concentrations increase towards the core).For the solid-phase data, however, when comparing individually for specific aggregate sections and chemical conditions, the modeled results lie within one standard deviation of experimental measurements for only a little over half of all aggregate sections.Overall, the solidphase model fits are convincing when viewed in the context of the large uncertainty introduced by uncontrollable biological factors.Deviations between experimental and model results are likely due to the inability to perfectly control experimental conditions rather than a misinterpretation of the biogeochemistry and physics that occur in the idealized soil aggregates.
The greatest error source in simulating aggregate experiments is the variability in cell densities and microbial activity since they factor into every rate law, are difficult to control experimentally, and their temporal evolution is hard to measure.In fact, the random error in the solid-phase determination of cell density, which was already large at the beginning of experiments (standard deviation was 42 % of the average measurement), increased in the measurements conducted at the conclusion of experiments.Whereas some of the error in measuring cell densities is likely to have been corrected by using cell density as a fitting parameter in simulations, this does not include potential errors introduced by temporal variations in cell-specific reaction rates that could not be measured or controlled.Additionally, cell densities were modeled as spatially homogenous within aggregates, whereas spatial variation is likely to have occurred, though eclipsed by the large uncertainty range of experimental observations.Based on the data collected for this work with E. cloacae and previous work with Shewanella putrefaciens CN-32 (Pallud et al., 2010a, b), we are confident that our procedure of aggregate construction leads to a distribution of cell density that is initially homogenous.However, spatial heterogeneity may occur over the course of experiments as cell populations grow or decay differentially in different sections of aggregates.For example, in our work with S. putrefaciens we found a trend at the conclusion of experiments indicating that cell densities may have been slightly higher (factor 2) at the exterior of aggregates (Pallud et al., 2010b), where elevated concentrations of electron donors may have led to more growth.Similarly, in the data presented here, the elevated CFU count value for the exterior section of the oxic aggregate run with high pyruvate (1.2 mM) and high selenate (0.8 mM) provides an indication of preferential growth in certain sections of the aggregate under certain chemical conditions.For example, our model shows that under oxic conditions, concentrations of both pyruvate and oxygen are expected to decrease steeply towards the core of aggregates.Since these are the substrates for aerobic growth, it is plausible that E. cloacae may have grown preferentially at the aggregate-solution boundary where concentrations of pyruvate and oxygen were the highest.
The cell density values (including the cell death term) used in the model are similar to the experimentally obtained values from bacterial cell counts (see above).The maximum cell specific rate constants which were the results of our model fitting (13 and 0.9 × 10 −17 mmol s −1 cell −1 for selenate and selenite reduction respectively) align well with what has previously been published for E. cloacae.Losi and Frankenberger (1997) found that the bacterium reduces 5 µmol g −1 of selenate over 30 min in a medium with 127 µM selenate and 2.8 mM pyruvate.Assuming a bacterial cell mass of 10 −12 g (Davis et al., 1973) and the rate law (Eq.1) and half-saturation constants (Table 1) given in the methods section, this would amount to a maximum cell specific rate of about 10 −17 mmol s −1 cell −1 , which is within an order of magnitude from our value.Similarly, Ma et al. (2009) found that E. cloacae reduces 1 mM concentrations of selenate and selenite at rates of 1.6 and 0.17 mmol s −1 g −1 , respectively, which following the same conversion as above amounts to maximum cell-specific rates of 10 and 1 × 10 −17 mmol s −1 cell −1 , respectively.These values are nearly identical to the ones that resulted from our model fitting (13 and 0.9 × 10 −17 mmol s −1 cell −1 for selenate and selenite reduction, respectively).
In summary, the simulations appear to capture the dynamics of the model aggregate system well and the kinetic parameters resulting from the model fitting represent the model organism (E.cloacae) in congruence with the literature.Thus, the reactive transport model presented here is a suitable framework for a theoretical exploration of aggregate scale impacts on selenium reduction.

Controls of aggregate scale gradients in reduced selenium
Previous experimental work showed substantial mm-scale radial intra-aggregate variations in reduced selenium distribution, with concentrations that increase from the advection boundary towards the core of artificial, 2.5 cm diameter aggregates (Kausch et al., 2012).With the reactive transport model presented here, we can now provide a mechanistic explanation of how this accumulation of reduced selenium at the core of aggregates occurs.
The rates of selenate reduction to selenite are relatively homogenous throughout the aggregate, yet selenite accumulates at the core of aggregates because increasing distance from the advective boundary surrounding the aggregate increases mass transfer limitations (lower selenite flux).The diffusive selenite flux is directed away from the aggregate's core throughout the aggregate (since the aggregate is a selenite source); the local magnitude of the flux, however, depends on the steepness of the selenite concentration gradient.Consequently, the diffusive flux is fastest at the aggregate solution boundary, since it is the boundary between the source domain (aggregate) on the one hand and a domain of zero production (free solution) on the other.At the boundary, this gradient is maintained by fast advective transport of selenite away from the aggregate-solution interface.Inside the aggregate, the concentration gradient is smoothed out by the diffusive flux, which decreases with increasing distance from the boundary along with the steepness of the gradient.The absolute steady-state selenite concentrations depend on the reduction rates since faster reduction shifts the equilibrium between production and diffusive export of selenite to higher concentrations.However, the decrease in concentrations towards the advective boundary occurs across the full experimental spectrum of reduction rates.Additionally, in the presence of oxygen, reduction is increasingly inhibited as oxygen concentrations increase towards the exterior of the aggregate.This explains why the experimentally observed ratio between reduced selenium concentrations in the core and exterior aggregate sections was greater under oxic than under anoxic conditions.The effect of differential oxygen inhibition in the exterior of aggregates is however secondary to the effect of diffusive mass transfer limitation in creating the observed gradients in reduced selenium.This is evidenced by the similarity in reduced selenium gradients observed under anoxic and oxic conditions.Another effect of oxygen is that by enabling aerobic respiration it increases the limitations imposed by the electron donor (pyruvate).Oxygen can thus indirectly inhibit selenium reduction at the aggregate core by limiting electron donor concentrations if the supply from the surrounding solution is small.In contrast to selenite concentrations, those of pyruvate decrease with increasing distance from the advective boundary.This explains why higher concentrations of pyruvate in the input solution lead to much higher reduced selenium concentrations at the core of aggregates under oxic conditions, but not under anoxic conditions.
The concentrations of solid-phase elemental selenium inside aggregates mirror those of selenite because selenite is the primary limiting factor in the production of elemental selenium (selenite reduction).The half-saturation constant for selenite in selenite reduction by E. cloacae is 0.7 mM (Ma et al., 2007); thus, at selenite concentrations representative of our experiments (∼ 0.01 and 0.05 mM for oxic an anoxic conditions, respectively), selenite reduction rates depend almost linearly on selenite concentrations (Eq.1).Consequently, the rates of selenite reduction are very heterogeneous inside aggregates with faster rates at the core where selenite accumulates.Once produced, elemental selenium will not move, and is thus not affected by transport, as selenite is.This also explains why elemental selenium concentrations in aggregates respond more strongly than selenite concentrations to conditions that promote the production of both forms of reduced selenium (i.e., for example, anoxic condition).Increasing selenite concentrations inside the aggregate will also lead to an increase of diffusive flux out of the aggregate.This negative feedback equilibrates concentrations of selenite in aggregates, which is why selenite distributions reached a steady state in all simulations.This is not the case for elemental selenium, which is not affected by diffusive transport and can thus grow boundlessly while selenite reduction occurs.The rates of elemental selenium production (selenite reduction) reach a selenite concentration dependent steady state, but the concentrations of elemental selenium keep increasing (at a steady rate) for as long as chemical input and microbial activity remain stable.

Impact of soil aggregate size on selenium retention
The same process that leads to the accumulation of reduced selenium at the core of aggregates causes the retention of reduced selenium to scale with aggregate size.Specifically, selenite (and consequently also elemental selenium) accumulates inside aggregates by virtue of diffusive mass transfer limitations that are more pronounced in larger aggregates.Smaller aggregates have a shorter average diffusion path length to the surrounding advective boundary (selenite sink) than larger aggregates.Consequently, selenite export from smaller aggregates to the surrounding solution is less limited than in the case of larger aggregates, and selenite will accumulate to a lesser extend at their cores.The impact of aggregate size is enhanced in the presence of oxygen since oxygen inhibition of selenium reduction will be stronger in smaller aggregates (shorter average diffusion path from the surrounding oxic solution).It is for this reason that, under bulk oxic conditions, the concentrations of reduced solid-phase selenium scaled more strongly with aggregate size than under anoxic conditions.Additionally, under oxic conditions, the interaction between electron donor concentrations and aggregate size becomes more important.The electron donor is needed both for aerobic respiration, which reduces oxygen concentrations in aggregates, and for selenium reduction.Thus, rates of electron donor consumption are faster under oxic conditions, and the electron donor concentrations inside the aggregate are lower than under anoxic conditions.At the core of aggregates, the electron donor can become limiting for selenium reduction, but less so in aggregates of smaller size where the average diffusion path to the advective boundary (electron donor supply) is shorter than for larger aggregates.As a result there is an interactive effect of electron donor concentrations and aggregate size under conditions that are oxic in bulk: the concentrations of reduced selenium retained in the solid-phase will scale more strongly with aggregate size if the supply of electron donor is larger.In summary, all other things equal, the amount of reduced selenium retained by a soil is expected to scale with increased aggregate size or electron donor concentrations, and with decreased aeration.Additionally, there are important interactive effects between these variables since, in the presence of oxygen, organic carbon (electron donor) concentrations increase the impact of aggregate size on reduced selenium concentrations.
With regards to selenium retention, the impact of aggregate size on elemental selenium concentrations is particularly relevant.Elemental selenium is insoluble and environmentally observed rates of oxidation are 3-4 orders of magnitude below those of reduction (Stolz et al., 2002(Stolz et al., , 2006)).Thus, soil elemental selenium is a good pool for long-term selenium retention, and the transformation of selenium oxyanions to solid elemental selenium has been suggested as a remedial strategy for contaminated sites (Darcheville et al., 2008;Dungan and Frankenberger, 1999;Oremland et al., 1991).Average concentrations of both solid-phase selenite and elemental selenium scale with aggregate size, but the effect is particularly pronounced for the elemental selenium fraction.Furthermore, in contrast to selenite, elemental selenium concentrations are not expected to reach a steady state according to our reactive transport model.It can therefore be expected that, as long as conditions are favorable to selenium reduction persist, the impact of aggregate size on elemental selenium content will keep increasing over time as the concentrations of elemental selenium build up.Increased soil aggregation may thus lead to a significant increase in long-term selenium retention.
Further studies are needed to establish whether the impact of aggregate size on selenium retention described here is significant at the field scale.A single, artificial aggregate surrounded by saturated flow cannot capture the full complexity of a structured soil, and relevant processes are likely to have been excluded.For example, selenite diffusing from aggregates may build up in macropores along the diffusion path, thereby decreasing the concentration gradient at the aggregate-macropore interface and reducing the impact of aggregate size on intra-aggregate selenite concentrations.Furthermore, the trends here discussed were obtained with selenate supplied in excess, while diffusive limitations in selenate supply may reduce selenium reduction in larger aggregates under field conditions.On the other hand, the dynamic saturation conditions in a natural surface soil may serve to increase the impact of aggregate size on selenium reduction through steeper redox gradients, as the macropores surrounding aggregates are expected to be filled with air rather than water for part of the time.The advantage of simplified models in studying complex natural systems is that they allow for the isolation of individual processes that may be inextricable in nature.The results here described point to general reactive transport mechanisms that may lead a soil with a larger mean aggregate size to retain more selenium.This would have implications for the management of irrigated seleniferous soils, since the fraction of macro-aggregates (> 0.25 mm in diameter) in a soil is sensitive to agricultural management on time scales of 2 yr or less (Degens, 1997).Notably, conservation tillage has been found to lead to significant increases in both mean soil aggregate size and organic matter content of surface soils (Angers et al., 1993;Carter, 1992).Similarly, the addition of manure can result in a rapid increase of mean aggregate size within the first two years of application (Whalen et al., 2003).It is particularly interesting to note that organic carbon input can increase aggregate size since the concomitant increase in electron donor availability and aggregate size is expected to synergistically enhance selenium retention according to our model.The beneficial effects of soil aggregation for plant growth have long been established (Dexter, 1988).If there were to be additional external benefits to enhancing soil structure for seleniferous soils, in the form of reduced selenium exports, conservation tillage and organic amendments may prove to be an attractive management technique in regions prone to irrigation-induced selenium contamination.Viewed in this context, our results suggest that field studies investigating the effect of conservation tillage and organic matter applications on selenium retention may be warranted in the search for management strategies to reduce selenium export.
In conclusion, our reactive transport model suggests that enhanced soil aggregation may increase the fraction of elemental selenium in aggregates and consequently increase the degree to which selenium is retained.This effect could be useful in the prevention of irrigation-induced selenium contamination.Managing seleniferous agricultural soils in a way that optimizes aggregation (i.e., reduced tillage and organic amendments) may reduce the amount of exported selenium and thus alleviate the impact of selenium pollution for aquatic ecosystems downstream of such soils.

Fig. 1 .
Fig. 1.Measured (diamonds) and simulated (solid lines) selenite concentrations in the reactor outflow under oxic (a) and anoxic (b) conditions.Colors correspond to specific combinations of input selenate and pyruvate concentrations (see legend).

Fig. 2 .
Fig. 2. Measured (solid bars) and simulated (hollow bars) concentrations of reduced selenium (selenite + elemental Se) in the solid phase of concentric aggregate sections (core, mid-section, exterior) under oxic (a) and anoxic (b) conditions.Error bars correspond to the standard deviation of replicate measurements (N = 2 for the core section and 3 otherwise).Colors correspond to specific combinations of input selenate and pyruvate concentrations (see legend).

Fig. 3 .
Fig. 3. Modeled concentration fields of pyruvate (a), O 2 (b), selenite (c), and Se(0) (d) within the aggregate (circle) and the surrounding solution after 192 h simulation time (steady state) in an oxic simulation with 0.3 mM pyruvate and 0.8 mM selenate input concentrations.Cold colors denote the absence of the species (minimum concentration is zero in all panels), whereas warm colors represent high concentrations, with maximum values as indicated on the panels.Dimensions L and R correspond to the length and radius of experimental reactors (3.7 and 2.55 cm respectively) and the dotted line represents the axis of cylindrical symmetry.The black rings correspond to the boundaries between the exterior, mid-, and core sections that were experimentally sampled.Molarities are expressed with respect to the pore/solution volume.The white arrows represent the flow field.

Fig. 4 .
Fig. 4. Modeled fields of diffusive selenite flux (a), total selenite flux (b), selenate reduction rates (c), and selenite reduction rates (d) within the aggregate (circle) and the surrounding solution after 192 h simulation time (steady state) in an oxic simulation with 0.3 mM pyruvate and 0.8 mM selenate input concentrations.Cold colors denote low values (minimum concentration is zero in all panels), whereas warm colors represent high concentrations, with maximum values as indicated on the panels.The black rings correspond to the boundaries between the exterior, mid-, and core sections that were experimentally sampled.Molarities are expressed with respect to the pore/solution volume.The white arrows in (a) and (b) represent flux vectors.

Fig. 5 .
Fig. 5. Results of oxic (a-d) and anoxic (e-h) simulation series evaluating the amount of reduced selenium per gram of the solid phase as a function of aggregate size for conditions with (a and e) 0.3 mM pyruvate/0.4 mM selenate, (b and f) 0.3 mM pyruvate/0.8 mM selenate, (c and g) 1.2 mM pyruvate/0.4 mM selenate,and (d and h) 1.2 mM pyruvate/0.8 mM selenate.The amount of selenite in the solid phase is shown in grey, while the amount of elemental selenium in the solid phase is shown in red.Each panel shows the results of 9 simulations carried out at different aggregates sizes; the space between simulations was linearly interpolated.

Table 1 .
Reaction rate parameters used in simulations.Cell-specific rate constants were established in a global fit to the experimental data shown here.