Models test our understanding of processes and can reach beyond the spatial and temporal scales of measurements. Multi-component Reactive Transport Models (RTMs), initially developed more than three decades ago, have been used extensively to explore the interactions of geothermal, hydrologic, geochemical, and geobiological processes in subsurface systems. Driven by extensive data sets now available from intensive measurements efforts, there is a pressing need to couple RTMs with other community models to explore non-linear interactions among the atmosphere, hydrosphere, biosphere, and geosphere. Here we briefly review the history of RTM development, summarize the current state of RTM approaches, and identify new research directions, opportunities, and infrastructure needs to broaden the use of RTMs. In particular, we envision the expanded use of RTMs in advancing process understanding in the Critical Zone, the veneer of the Earth that extends from the top of vegetation to the bottom of groundwater. We argue that, although parsimonious models are essential at larger scales, process-based models are the only way to explore the highly nonlinear coupling that characterizes natural systems. We present seven testable hypotheses that emphasize the unique capabilities of process-based RTMs for (1) elucidating chemical weathering and its physical and biogeochemical drivers; (2) understanding the interactions among roots, microorganisms , carbon, water, and minerals in the rhizosphere; (3) assessing the effects of heterogeneity across spatial and temporal scales; and (4) integrating the vast quantity of novel data, including-omics‖ data (genomics, transcriptomics, proteomics, metabolomics), elemental concentration and speciation data, and isotope data into our understanding of complex earth systems. With strong support from data-driven sciences, we are now in an exciting era where integration of RTM framework into other community models will facilitate process understanding across disciplines and across scales.

Why do solute concentrations in streams remain largely constant while discharge varies by orders of magnitude? We used a new hydrological land surface and reactive transport code, RT‐Flux‐PIHM, to understand this long‐standing puzzle. We focus on the nonreactive chloride (Cl) and reactive magnesium (Mg) in the Susquehanna Shale Hills Critical Zone Observatory (SSHCZO). Simulation results show that stream discharge comes from surface runoff (Qs), soil lateral flow (QL), and deeper groundwater (QG), with QL contributing >70%. In the summer, when high evapotranspiration dries up and disconnects most of the watershed from the stream, Cl is trapped along planar hillslopes. Successive rainfalls connect the watershed and mobilize trapped Cl, which counteracts dilution effects brought about by high water storage (Vw) and maintains chemostasis. Similarly, the synchronous response of clay dissolution rates (Mg source) to hydrological conditions, maintained largely by a relatively constant ratio between “wetted” mineral surface area Aw and Vw, controls Mg chemostatic behavior. Sensitivity analysis indicates that cation exchange plays a secondary role in determining chemostasis compared to clay dissolution, although it does store an order‐of‐magnitude more Mg on exchange sites than soil water. Model simulations indicate that dilution (concentration decrease with increasing discharge) occurs only when mass influxes from soil lateral flow are negligible (e.g., via having low clay surface area) so that stream discharge is dominated by relatively constant mass fluxes from deep groundwater that are unresponsive to surface hydrological conditions.

Model development in hydrology and geochemistry has been advancing separately with limited integration. We developed a watershed hydrogeochemical code RT‐Flux‐PIHM to understand complex interactions between hydrological processes (PIHM), land‐surface processes (FLUX—Noah Land Surface Model), and multicomponent subsurface reactive transport (RT). The RT module simulates geochemical processes including aqueous complexation, surface complexation, mineral dissolution and precipitation, and cation exchange. The RT module is verified against the widely used reactive transport code CrunchFlow. The code uses semidiscrete finite volume method and irregular gridding and offers data harvesting capabilities from national databases. The application of RT‐Flux‐PIHM is demonstrated in the Susquehanna Shale Hills Critical Zone Observatory (SSHCZO). We aim to understand key processes that govern hydrogeochemical dynamics of the nonreactive chloride and reactive magnesium. Simulation results indicate that watershed characteristics, in particular topography, dictate the spatial distributions of water content and soil dissolution rates. Ion exchange provides buffering capacities and leads to a hysteresis loop of concentration and discharge relationship of magnesium, which differs from the open hysteresis of chloride. RT‐Flux‐PIHM offers physics‐based modeling capabilities to integrate the vast amount of water and chemistry data that have now become available, to differentiate the relative importance of competing processes, and to test hypotheses at the interface of hydrology and geochemistry.

We aim to understand the scale-dependent evolution of uranium bioreduction during a field experiment at a former uranium mill site near Rifle, Colorado. Acetate was injected to stimulate Fe-reducing bacteria (FeRB) and to immobilize aqueous U(VI) to insoluble U(IV). Bicarbonate was coinjected in half of the domain to mobilize sorbed U(VI). We used reactive transport modeling to integrate hydraulic and geochemical data and to quantify rates at the grid block (0.25 m) and experimental field scale (tens of meters). Although local rates varied by orders of magnitude in conjunction with biostimulation fronts propagating downstream, field-scale rates were dominated by those orders of magnitude higher rates at a few selected hot spots where Fe(III), U(VI), and FeRB were at their maxima in the vicinity of the injection wells. At particular locations, the hot moments with maximum rates negatively corresponded to their distance from the injection wells. Although bicarbonate injection enhanced local rates near the injection wells by a maximum of 39.4%, its effect at the field scale was limited to a maximum of 10.0%. We propose a rate-versus-measurement-length relationship (log R' = -0.63 log L - 2.20, with R' in μmol/mg cell protein/day and L in meters) for orders-of-magnitude estimation of uranium bioreduction rates across scales.

Geologic carbon sequestration represents a promising option for carbon mitigation. Injected CO(2), however, can potentially leak into water systems, increase water acidity, and mobilize metals. This study used column experiments to quantify the effects of environmental controls on cadmium desorption during CO(2) leakage in subsurface systems without ambient flow. Results show that fast leakage rates are responsible for earlier and larger amounts of Cd desorption. Long weathering time of Cd laden clay leads to low Cd desorption. Calcite content as low as 10% can mitigate the effect of pH reduction and result in zero Cd desorption. Increasing the salinity of the leaking fluid has a relatively minor effect, primarily due to the offsetting impacts of an increased extent of ion exchange and the decrease in CO(2) solubility (and therefore acidity). This work systematically quantifies, for the first time, the effects of environmental controls on Cd desorption and points to key parameters for risk assessment associated with metal mobilization during CO(2) leakage.