In the latter study, we found increasing concentrations of reduced selenium towards the core of aggregates containing E. cloacae 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 . 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 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 . To this end, we reproduced our experimental results to validate the model, and then simulated selenate and selenite reduction in aggregates ranging from 1-2.5 cm in diameter.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 .We used data on bacterial cell densities and concentrations of reduced selenium within three concentric sections of the aggregates,flower bucket as well as data on temporal evolution of selenite and flow-tracer 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 constrain this biological factors in the reactive transport model by checking whether cell density varied between concentric sections of the aggregates and between the beginning and end of experiments. Whereas the concentrations of pyruvate, selenate, and selenite approached steady state between 64 and 192 h of simulation time , those of oxygen reach steady-state within 32 h . Concentrations of elemental selenium on the other hand kept increasing inside the aggregate domain throughout the simulation time . Over the course of simulations, all compounds developed smooth gradients of approximately radial symmetry within the aggregate domain . However, whereas concentrations of reactants decreased towards the core of the aggregate domain, those of products increased towards the core . The gradients of selenite, elemental selenium and, except in the case of selenate, also of reactants, were steeper under oxic than under anoxic conditions . 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 . Diffusive and advective transport, respectively, dominated inside the aggregate and in the surrounding solution. This is illustrated for the case of selenite in Fig. 5. 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/ at the interface between the aggregate and the surrounding solution .
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/ . Similar to solute concentrations, reaction rates of selenate and selenite reduction reached a steady state after 60 and 150 h of simulation time, respectively . The rates of selenate reduction approached a steady-state that is relatively homogenous throughout the aggregate . In contrast, the rates of selenite reduction approached a steady-state that mirrored the concentration field of selenite . In the example simulation the steady state rates of selenate reduction were 130-190 nM/s everywhere in the aggregate domain , whereas those of selenite reduction smoothly increased from 0 at the aggregate exterior to 12 nM/s at the core .Simulations in which the diameter of the spherical aggregate was varied between 1 and 2.5 cm showed that the aggregate solid phase concentrations of reduced selenium scale with aggregates size under all chemical conditions investigated . 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 in aggregates of 1 cm diameter and 4-17 nmol/g in aggregates of 2.5 cm diameter, which corresponds to a difference by a factor of 6-11 between the two size extremes . 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 and 14-36 nmol/g, respectively . 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 over the investigated size range . Thus the fraction of elemental selenium in aggregates varied greatly with aggregate size under oxic conditions . In contrast, under anoxic conditions elemental selenium makes up around 60% 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 . However, under oxic conditions the curves corresponding to high pyruvate concentration scenarios were markedly steeper towards the larger end of the aggregate size spectrum than those for low pyruvate concentration scenarios .The reactive transport model we 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, and the spatial distribution of reduced selenium inside the aggregates . 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,plastic flower bucket the solid phase 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 bio-geochemistry 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 , 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 which 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 , 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 at the exterior of aggregates , 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 and high selenate 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 used in the model are similar to the experimentally obtained values from bacterial cell counts . The maximum cell specific rate constants which were the result of our model fitting align well with what has previously been published for E. cloacae. Losi and Frankenberger found that the bacterium reduces 5 µmol/g of selenate over 30 minutes in a medium with 127 µM selenate and 2.8 mM pyruvate. Assuming a bacterial cell mass of 10-12 g and the rate law and half-saturation constants given in the methods section, this would amount to a maximum cell specific rate of about 10-17 mmol/s/cell, which is within an order of magnitude from our value. Similarly, Ma et al. found that E. cloacae reduces 1 mM concentrations of selenate and selenite at rates of 1.6 and 0.17 µmol/s/g respectively, which following the same conversion as above amounts to maximum cell specific rates of 10 and 1×10-17 mmol/s/cell respectively. These values are nearly identical to the ones that resulted from our model fitting . 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 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.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 . 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 . The diffusive selenite flux is directed away from the aggregate’s core throughout the aggregate , 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 on the one hand and a domain of zero production 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.