Abstract The weathering of sedimentary rocks in high‐elevation catchments influences freshwater quality and the global carbon cycle. While individual biogeochemical mechanisms involved in this process are relatively well understood, quantifying their contributions to solute export and carbon fluxes under natural, transient conditions remains challenging. Here, we implement a numerical multidimensional and multiphase model to simulate coupled hydrological and biogeochemical processes in a shale‐underlain, snow‐dominated hillslope in the Rocky Mountains, Colorado. The model captures the dynamic interplay between soil respiration, mineral weathering, and climate‐driven hydrological forcing, reproducing observed soil CO2 dynamics, groundwater chemistry, and subsurface flow. Our results reveal that seasonal snowmelt enhances carbonate weathering by promoting the infiltration of CO2‐rich water to depth, while pyrite oxidation is primarily sensitive to low water saturation that facilitates O2 diffusion through the regolith. Topography modulates the spatial distribution of shale weathering, as steeper slopes enhance lateral drainage, favoring the delivery of reactants to greater depths. While shale weathering at our site acts as a transient carbon sink, with silicates and carbonates buffering acidity and promoting atmospheric CO2 consumption (1% of soil‐derived CO2), the exported dissolved inorganic carbon is predominantly geogenic (∼73%). Consequently, when accounting for long‐term marine carbonate precipitation. The current weathering regime represents a net source of carbon to the atmosphere. The oxidation of pyrite and petrogenic organic carbon together release approximately 0.9 mol·m−2·yr−1 of CO2. Our findings highlight the role of topography, hydroclimate, and the coupling between acid‐base reactions in shaping the carbon balance and the solute exports in mountainous critical zones.
The growing demand for hydrogen necessitates innovative production strategies beyond conventional fossil fuelderived methods. Stimulated geologic hydrogen generation through serpentinization of ultramafic rock is a promising alternative. However, the feedback mechanisms between fluid transport, mineral reactions, and rock deformation that govern the kinetics and evolution of serpentinization remain poorly constrained. This study investigates hydrogen yields and mineral phase transformations during laboratory serpentinization of peridotite. Experiments conducted at 90 degrees C and 0.1 MPa, and 200 degrees C and 6 MPa, show hydrogen yields of 0.5 to 0.9 & micro;mol/g rock and 20.8 to 31.2 & micro;mol/g rock after 30 days, respectively. High-resolution X-ray computed tomography captures the spatiotemporal (4D) evolution of serpentine vein networks during serpentinization, revealing serpentine growth primarily along pre-existing veins and sustaining fluid-mineral interactions through submicron vein permeability. Quantitative volume analyses show that serpentine volume fraction increases by 6.7 +/- 0.93 to 10.2 +/- 1.42 vol% from an initial value of 38 to 39 vol%, resulting in a volumetric strain of + 0.005 +/- 0.00025 to + 0.006 +/- 0.00012 (positive for volume increase) under zero effective stress. These results indicate that serpentine vein growth during serpentinization is primarily driven by mineral replacement at olivineserpentine interfaces rather than by stress-induced fracturing followed by mineral infill, with fluid transport playing an important role in reaction progression. Overall, the findings highlight serpentinization as a coupled hydro-thermo-chemo-mechanical process that can occur over short timescales and at relatively low temperatures in peridotite, features that are desirable for scalable geologic hydrogen stimulation and production.
Supercritical (sc) CO2 in geologic carbon sequestration (GCS) can chemically and mechanically deteriorate wellbore cement, raising concerns for long-term operations. In contrast to the conventional view of "sulfate attack" on cement, we found that adding 0.15 M sulfate to the acidic brine can significantly reduce the impact of scCO2 attack on Portland cement, resulting in stronger cement than that found in a sulfate-free system. Scanning electron microscopy revealed a decreased total attack depth in reacted cement in the presence of sulfate. With a newly defined minimum porosity term in reactive transport modeling, our model suggests that sulfate caused CaCO3 to fill more nanopore spaces in the cement. Small angle X-ray scattering experiments also showed that sulfate can decrease the pore sizes of the carbonate layer. The results suggest that the interactions between sulfate and cement can generate a less porous CaCO3 layer, which better resists acidic brine. Using this mechanism as a proof-of-concept, we tested the incorporation of sodium sulfate into Portland cement and synthesized new cement composites that show stronger resistance against scCO2 attacks. These newly discovered interfacial interactions between CaCO3 and sulfate provide new insights into engineering mechanically strong and green materials for safer GCS.
Land surface models consider the exchange of water, energy, and carbon along the soil‐canopy‐atmosphere continuum, which is challenging to model due to their complex interdependency and associated challenges in representing and parameterizing them. Differentiable modeling provides a new opportunity to capture these complex interactions by seamlessly hybridizing process‐based models with deep neural networks (DNNs), benefiting both worlds, that is, the physical interpretation of process‐based models and the learning power of DNNs. Here, we developed a differentiable land model, JAX‐CanVeg. The new model builds on the legacy CanVeg by incorporating advanced functionalities through JAX in the graphic processing unit support, automatic differentiation, and integration with DNNs. We demonstrated JAX‐CanVeg's hybrid modeling capability by applying the model at four flux tower sites with varying aridity. To this end, we developed a hybrid version of the Ball‐Berry equation that emulates the water stress impact on stomatal closure to explore the capability of the hybrid model in (a) improving the simulations of latent heat fluxes and net ecosystem exchange , (b) improving the optimization trade‐off when learning observations of both and , and (c) benefiting a multi‐layer canopy model setup. Our results show that the proposed hybrid model improved the simulations of and at all sites, with an improved optimization trade‐off over the process‐based model. Additionally, the multi‐layer canopy set benefited hybrid modeling at some sites. Anchored in differentiable modeling, our study provides a new avenue for modeling land‐atmosphere interactions by leveraging the benefits of both data‐driven learning and process‐based modeling.
The modeling and simulation of the Cement-clay Interaction-Diffusion field (CI-D) experiment at the Mont Terri site in Switzerland presented here demonstrates that it is possible to capture the multiscale physical and chemical features of natural and engineered barrier systems for radionuclides. The simulations are successfully carried out with the newly developed CrunchODiTi high-performance computing software that accounts for multiple continua, including a continuum representing the electrical double layer (EDL) developed along negatively charged clay particles in clay rock. The simulation also accounts for both the complex three-dimensional (3D) geometry, expected as the norm in a geological waste repository, and the anisotropy of the geological formation. In addition, the high resolution of the model makes it possible to include "skin effects" developed at the interface between highly reactive materials, in this case between the high pH cement and the circumneutral but electrostatic Opalinus Clay. The successful history matching with the field experiment demonstrates that the distinct geochemical and physical properties of the cement and the Opalinus Clay in the CI-D experiment can be accounted for. Such analyses are essential for developing a defensible safety case for the underground storage of radioactive waste.
Radionuclide transport in smectite clay barrier systems used for nuclear waste disposal is controlled by diffusion, with adsorption significantly retarding transport rates. While a relatively minor component of spent nuclear fuel, Se-79 is a major driver of the safety case for spent fuel disposal due to its long half-life (3.3 x 10(5) yr) and its low adsorption to clay (K-D < 10 L/kg), thus a thorough understanding of Se diffusion through clay is critical for understanding the long-term safety of spent fuel disposal systems. Through-diffusion experiments with tritiated water (HTO, conservative tracer) and Se(VI) were conducted with a well-characterized, purified montmorillonite source clay (SWy-2) under a constant ionic strength (0.1 M) and three different electrolyte compositions: Na+, Ca2+, and a Na (+) -Ca2+ mixture at pH 6.5 in order to probe the effects of electrolyte composition and interlayer cation composition on clay microstructure, Se(VI) aqueous speciation, and ultimately diffusion. The results were modeled using a reactive transport modeling approach to determine values of porosity (epsilon), D-e (effective diffusion coefficient), and K-D (distribution coefficient for adsorption). HTO diffusive flux was higher in Ca-montmorillonite (D-e = 1.68 x 10(-10) m(2) s(-1)) compared to Na-montmorillonite (D-e = 7.83 x 10(-11) m(2) s(-1)). This increase in flux is likely due to a greater degree of clay layer stacking in the presence of Ca2+ compared to Na+, which leads to larger inter-particle pores. Overall, the Se(VI) flux was much lower than the HTO flux due to anion exclusion, with Se(VI) flux following the order Ca (D-e = 1.03 x 10(-11) m(2) s(-1)) > Na-Ca (D-e = 2.12 x 10(-12) m(2) s(-1)) > Na (D-e = 1.28 x 10(-12) m(2) s(-1)). These differences in Se(VI) flux are due to a combination of factors, including (1) larger accessible porosity in Ca-montmorillonite due to clay layer stacking and smaller electrostatic effects compared to Na-montmorillonite, (2) larger accessible porosity for neutral-charge CaSeO4 species which makes up 32% of aqueous Se(VI) in the pure Ca system, and (3) possibly higher Se(VI) adsorption for Ca-montmorillonite. Through a combination of experimental and modeling work, this study highlights the compounding effects that electrolyte and counterion compositions can have on radionuclide transport through clay. Diffusion models that neglect these effects are not transferable from laboratory experimental conditions to in situ repository conditions.
We present a new pore-scale model for multicomponent advective-diffusive transport with coupled mineral dissolution and precipitation. Both dissolution and precipitation are captured simultaneously by introducing a phase transformation vector field representing the direction and magnitude of the overall phase change. An effective viscosity model is adopted in simulating fluid flow during mineral dissolution-precipitation that can accurately capture the velocity field without introducing any empirical parameters. The proposed approach is validated against analytical solutions and interface tracking simulations in simplified structures. After validation, the proposed approach is employed in modeling realistic rocks where mineral dissolution and precipitation are dominant at different locations. We have identified three regimes for mineral dissolution-precipitation coupling: (a) compact dissolution-precipitation where dissolution is dominant near the inlet and precipitation is dominant near the outlet, (b) wormhole dissolution with clustered precipitation where dissolution generates wormholes in the main flow paths and precipitation clogs the secondary flow paths, and (c) dissolution dominant where all solid grains are gradually dissolved. In the three regimes, the proposed approach provides reliable porosity-permeability relationships that cannot be described well by traditional macroscale models. We find that the permeability can increase while the overall porosity decreases when the main flow paths are expanded by dissolution and adjacent pore spaces are clogged by precipitation. A new pore-scale model is developed for multicomponent reactive flow with coupled mineral dissolution and precipitation The proposed model is validated against analytical solutions and reference simulations, and then applied to realistic rocks Three dissolution-precipitation coupling regimes with different flow patterns and porosity-permeability relationships have been identified
Models have historically represented fractured porous media with continuum descriptions that characterize the media using bulk parameters. The impact of small-scale features is not captured in these models, although they may be controlling the performance of subsurface applications. Pore-scale models can simulate processes in small-scale features by representing the pore space geometry explicitly but are computationally expensive for large domains. The alternative multiscale approach entails the combination of pore-scale and continuum-scale descriptions in a single framework. We use Chombo-Crunch, a computational capability that discretizes complex geometries with an adaptive, embedded boundary method to contrast these two approaches. Chombo-Crunch takes advantage of recent computational performance and memory bandwidth improvements resulting from the emergence of exascale computing resources. These combined improvements enable the efficient simulation of reactive transport in fractured media with a high degree of fidelity and the ability to capture the control small-scale processes exert on the overall medium evolution.
Climate influences near-surface biogeochemical processes and thereby determines the partitioning of carbon dioxide (CO 2 ) in shale, and yet the controls on carbon (C) weathering fluxes remain poorly constrained. Using a dataset that characterizes biogeochemical responses to climate forcing in shale regolith, we implement a numerical model that describes the effects of water infiltration events, gas exchange, and temperature fluctuations on soil respiration and mineral weathering at a seasonal timescale. Our modeling approach allows us to quantitatively disentangle the controls of transient climate forcing and biogeochemical mechanisms on C partitioning. We find that ~3% of soil CO 2 (1.02 mol C/m 2 /y) is exported to the subsurface during large infiltration events. Here, net atmospheric CO 2 drawdown primarily occurs during spring snowmelt, governs the aqueous C exports (61%), and exceeds the CO 2 flux generated by pyrite and petrogenic organic matter oxidation (~0.2 mol C/m 2 /y). We show that shale CO 2 consumption results from the temporal coupling between soil microbial respiration and carbonate weathering. This coupling is driven by the impacts of hydrologic fluctuations on fresh organic matter availability and CO 2 transport to the weathering front. Diffusion-limited transport of gases under transient hydrological conditions exerts an important control on CO 2(g) egress patterns and thus must be considered when inferring soil CO 2 drawdown from the gas phase composition. Our findings emphasize the importance of seasonal climate forcing in shaping the net contribution of shale weathering to terrestrial C fluxes and suggest that warmer conditions could reduce the potential for shale weathering to act as a CO 2 sink.
The constitutive relations of the Richardson‐Richards equation encode the macroscopic properties of soil water retention and conductivity. These soil hydraulic functions are commonly represented by models with a handful of parameters. The limited degrees of freedom of such soil hydraulic models constrain our ability to extract soil hydraulic properties from soil moisture data via inverse modeling. We present a new free‐form approach to learning the constitutive relations using physically constrained neural networks. We implemented the inverse modeling framework in a differentiable modeling framework, JAX, to ensure scalability and extensibility. For efficient gradient computations, we implemented implicit differentiation through a nonlinear solver for the Richardson‐Richards equation. We tested the framework against synthetic noisy data and demonstrated its robustness against varying magnitudes of noise and degrees of freedom of the neural networks. We applied the framework to soil moisture data from an upward infiltration experiment and demonstrated that the neural network‐based approach was better fitted to the experimental data than a parametric model and that the framework can learn the constitutive relations.
Most of the available data on diffusion in natural clayey rocks consider tracer diffusion in the absence of a salinity gradient despite the fact that such gradients are frequently found in natural and engineered subsurface environments. To assess the role of such gradients on the diffusion properties of clayey materials, throughdiffusion experiments were carried out in the presence and absence of a salinity gradient using salt-diffusion and radioisotope tracer techniques. The experiments were carried out with vermiculite samples that contained equal proportions of interparticle and interlayer porosities so as to assess also the role played by the two types of porosities on the diffusion of water and ions. Data were interpreted using both a classical Fickian diffusion model and with a reactive transport code, CrunchClay that can handle multi-porosity diffusion processes in the presence of charged surfaces. By combining experimental and simulated data, we demonstrated that (i) the flux of water diffusing through vermiculite interlayer porosity was minor compared to that diffusing through the interparticle porosity, and (ii) a model considering at least three types of porous volumes (interlayer, interparticle diffuse layer, and bulk interparticle) was necessary to reproduce consistently the variations of neutral and charged species diffusion as a function of salinity gradient conditions.
Residual trapping is an important process that affects the efficiency of cyclic storage and withdrawal and in-situ production of hydrogen in geological media. In this study, we have conducted pore-scale modeling to investigate the effects of pore geometry and injection rate on the occurrence and efficiency of residual trapping via dead-end bypassing. We begin our theoretical and numerical analyses using a single rectangular pore to understand the key controls in bypassing. We further investigated two factors affecting bypassing: (a) a continuous cycle of injection-extraction of H2, and (b) variable pore geometry. Based on our pore-scale simulations, we found that: (a) a higher pore height/width ratio (h/w) and a higher injection rate cause more residual trapping, which is unfavorable for withdrawal of H2; (b) the trapping percentage increases with the h/w first and then decreases after h/w reaches 0.5; (c) and a converging-shaped pore can result in less trapping volume. Based on a theoretical comparison of the residual trapping behavior of H2 and CO2, we discuss the mechanisms that are applicable to CO2 residual trapping and the possibility of developing engineering controls of H2 storage and production.
We computationally explore the relationship between surface–subsurface exchange and hydrological response in a headwater-dominated high elevation, mountainous catchment in East River Watershed, Colorado, USA. In order to isolate the effect of surface–subsurface exchange on the hydrological response, we compare three model variations that differ only in soil permeability. Traditional methods of hydrograph analysis that have been developed for headwater catchments may fail to properly characterize catchments, where catchment response is tightly coupled to headwater inflow. Analyzing the spatially distributed hydrological response of such catchments gives additional information on the catchment functioning. Thus, we compute hydrographs, hydrological indices, and spatio-temporal distributions of hydrological variables. The indices and distributions are then linked to the hydrograph at the outlet of the catchment. Our results show that changes in the surface–subsurface exchange fluxes trigger different flow regimes, connectivity dynamics, and runoff generation mechanisms inside the catchment, and hence, affect the distributed hydrological response. Further, changes in surface–subsurface exchange rates lead to a nonlinear change in the degree of connectivity—quantified through the number of disconnected clusters of ponding water—in the catchment. Although the runoff formation in the catchment changes significantly, these changes do not significantly alter the aggregated streamflow hydrograph. This hints at a crucial gap in our ability to infer catchment function from aggregated signatures. We show that while these changes in distributed hydrological response may not always be observable through aggregated hydrological signatures, they can be quantified through the use of indices of connectivity.
The reactive transport code CrunchClay was used to derive effective diffusion coefficients ( D e ), clay porosities ( ε ), and adsorption distribution coefficients ( K D ) from through-diffusion data while considering accurately the influence of unavoidable experimental biases on the estimation of these diffusion parameters. These effects include the presence of filters holding the solid sample in place, the variations in concentration gradients across the diffusion cell due to sampling events, the impact of tubing/dead volumes on the estimation of diffusive fluxes and sample porosity, and the effects of O-ring-filter setups on the delivery of solutions to the clay packing. Doing so, the direct modeling of the measurements of (radio)tracer concentrations in reservoirs is more accurate than that of data converted directly into diffusive fluxes. While the above-mentioned effects have already been described individually in the literature, a consistent modeling approach addressing all these issues at the same time has never been described nor made easily available to the community. A graphical user interface, CrunchEase, was created, which supports the user by automating the creation of input files, the running of simulations, and the extraction and comparison of data and simulation results. While a classical model considering an effective diffusion coefficient, a porosity and a solid/solution distribution coefficient ( D e – ε – K D ) may be implemented in any reactive transport code, the development of CrunchEase makes it easy to apply by experimentalists without a background in reactive transport modeling. CrunchEase makes it also possible to transition more easily from a D e – ε – K D modeling approach to a state-of-the-art process-based understanding modeling approach using the full capabilities of CrunchClay, which include surface complexation modeling and a multi-porosity description of the clay packing with charged diffuse layers.
The weathering of shale exerts an important control on the hydrochemical fluxes to river systems, thus influencing the global carbon, nutrient, and geochemical cycles. However, the quantitative understanding of shale weathering and its impact on global biogeochemical cycles remains inadequate due to the com-plex interplay between hydrological, biogeochemical, and physical processes. In this study, we develop a novel modeling approach to quantitatively interpret the long-term chemical weathering of shale occur-ring since the last glaciation period (15,000 years) leading to the present geochemical conditions. The model explicitly considers processes occurring across multiple phases and involved in the weathering, including: (i) the infiltration of meteoric water, (ii) the interactions between the water and the mineral assemblage via dissolution/precipitation reactions, (iii) the microbially-mediated oxidation of organic matter, (iv) the evolution of porosity induced by mineral reactions, and (v) the exchange of gases between the subsurface and the atmosphere. To implement and test our model, we conduct this study at a well -instrumented hillslope underlain by Mancos Shale and located at the East River study site, Western Colorado. Consistent with field observations, the model successfully reproduces the stratified weathering front, the complex spatial distribution of organic carbon, and the gaseous emissions of carbon dioxide from the subsurface. Model simulations show that aerobic respiration exerts a fundamental control on the weathering of shale. While previous studies have highlighted the diffusion of oxygen from the atmo-sphere as the primary mechanism for shale weathering, our model simulations demonstrate that aerobic respiration limits the propagation of oxygen in the shallow subsurface, and thereby inhibits the dissolu-tion of pyrite at depth. Aerobic respiration is particularly favored in the top soil horizon due to the con-stant flux of oxygen from the atmosphere, the replenishment of fresh litter/plant-derived organic matter, and to a lesser extent the presence of fossil shale-associated organic matter. The acidic pore water gen-erated through aerobic respiration within the shallow subsurface is transported to greater depths, where it sustains the dissolution of carbonate (dolomite in this example). Overall, our results demonstrate that the evolution of pyrite and carbonate depletion fronts are significantly different, and primarily depend on the ability of microorganisms to carry out microbial respiration and the various transport pathways of reactants controlling the mineral reactions under partially saturated conditions.& COPY; 2022 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).