In this study we use numerical experiments to analyze reactive flow and transport behavior in discrete intersecting fracture networks, focusing on (a) how reaction-induced changes in physical and chemical properties affect flow connectivity and (b) how fracture networks developed in the Earth's critical zone contribute to carbon sequestration via mineral weathering reactions. In the first part of the study, we used two-dimensional reactive flow and transport simulations to analyze the impacts of mixing in a natural discrete fracture network. We concluded that reaction-induced changes can substantially alter the flow connectivity, especially at fracture intersections. The second set of simulations considered the problem of natural weathering of fractured mafic and ultramafic rocks in the partially saturated Earth's critical zone as a function of infiltration rates, fracture permeability, and partially saturated flow parameters. As a model system, we considered an incongruent reaction network with dissolution of forsterite and precipitation of magnesite. The behavior is complex in terms of the rate-controlling processes because of the multicomponent nature of the system as shown by the grid Peclet number: the CO2 behavior is gas diffusion-controlled in the partially saturated zone, while the rate of water flow via the Damkӧhler number controls Mg2+ transport through the fracture network. The amounts of carbon that can be trapped are modest, but the naturally fractured domain considered here provides a useful “base case” against which various engineered solutions can be compared.

Flow-only and coupled flow-diffusion experiments at 95 °C and 100 bars pCO2 carried out in micro-capillary tubes packed with forsterite mineral grains were used to constrain a coupled classical heterogeneous nucleation and crystal growth reactive transport model describing the formation of secondary magnesite. The study made use of a novel experimental setup in which one capillary tube for flow is connected via a three-way tee to a perpendicular capillary tube sealed at the distal end in which only molecular diffusion is allowed to occur—an experimental analogue of a single fracture-rock matrix system. While the high flux of CO2 bearing fluids and their low pH did not result in the formation of secondary carbonates in the flow-dominated channels, as much as 2.7% magnesite formed in a diffusion-controlled capillary tube sample after 300 hours of reaction. The precipitation of secondary magnesite was not uniformly distributed, however, but showed a distinctive peak shaped pattern along the sample reacted as quantified by RAMAN spectroscopic analysis. About 50% of the total magnesite precipitation formed within a narrow interval of a few millimeters close to the middle of the 3 cm diffusion sample. To simulate the behavior of the diffusion–reaction column, and in particular the pronounced mm scale peak at approximately 1.8 cm, an interfacial free energy of approximately 70 mJ·m−2 combined with a relatively high crystal growth rate for secondary magnesite was required. In agreement with the observations based on RAMAN spectroscopy, the simulations suggest that more than 50% of Mg2+ dissolved from the primary forsterite precipitated as secondary magnesite. The magnesite precipitation at the central peak position in the diffusion sample showed such rapid growth that it created a local minimum in pore fluids Mg2+ concentration close to the peak, creating a “nucleation shadow” that focused continued growth there as nearby regions had no nucleation seeds available for crystal growth. Continued crystal growth on the initial magnesite band acted to further increase the reactive surface area at this point, thus enhancing the spatially focused crystal growth rate and creating a positive feedback leading to pattern formation. This highlights the potentially critical role of an initial nucleation event in controlling mineral precipitation patterns in subsurface porous media, patterns that may determine how the pore structure subsequently evolves physically and chemically over time due to reactive flow and transport.

A model based on the code CrunchClay is presented for a fracture-clay matrix system that takes electrostatic effects on transport into account. The electrostatic effects on transport include those associated with the development of a diffusion potential as captured by the Nernst-Planck equation, and the formation of a diffuse layer bordering negatively charged clay particles within which partial anion exclusion occurs. The model is based on a dual continuum formulation that accounts for diffuse layer and bulk water pore space, providing a more flexible framework than is found in the classical mean electrostatic potential models. The diffuse layer model is obtained by volume averaging ion concentrations in the Poisson-Boltzmann equation, but also includes the treatment of longitudinal transport within this continuum. The calculation of transport within the bulk and diffuse layer porosity is based on a new formulation for the Nernst-Planck equation that considers averaging of diffusion coefficients and accumulation factors at grid cell interfaces. Equations for function residuals and the associated Jacobian matrix are presented such that the system of nonlinear differential-algebraic equations can be solved with Newton’s method. As an example, we consider a 2D system with a single discrete fracture within which flow and advective transport occurs that is coupled to diffusion in the clay-rich matrix. The simulation results demonstrate the lack of retardation for anions (e.g., 36Cl−) of the contaminant plume within the fracture flow system because they are largely excluded from the charged clay rock, while the migration of cations (e.g., 90Sr++) is more strongly attenuated. The diffusive loss of divalent cations in particular from the fracture is accentuated by their accumulation in the diffuse layer within the clay-rich matrix.

Molecular diffusion of dissolved species is a fundamental mass transport process affecting many environmental and technical processes. Whereas diffusive transport of single tracers can be described by Fick's law, a multicomponent approach based on the Nernst-Planck equation is required for charge-coupled transport of ions. The numerical solution of the Nernst-Planck equation requires special attention with regard to properties that are required at interfaces of numerical cells when using a finite difference or finite volume method. Weighted arithmetic and harmonic averages are used in most codes that can solve the Nernst-Planck equation. This way of averaging is correct for diffusion coefficients but inappropriate for solute concentrations at interfaces. This averaging approach leads to charge balance problems and thus to numerical instabilities near interfaces separating grid volumes with contrasting properties. We argue that a logarithmic-differential average should be used. Here this result is generalized, and it is demonstrated that it generally leads to improved numerical stability and accuracy of concentrations computed near material interfaces. It is particularly relevant when modeling semipermeable clay membranes or membranes used in water treatment processes.

Microscale numerical modeling is potentially an effective approach for understanding salt creep, but it is quite challenging because natural salt grain boundaries can have complex and evolving geometry, and because contacts between these deformable salt mineral grains are dynamic as a result of compaction, chemical reaction, fluid flow, and heat transfer. In this study, we have overcome these challenges and developed a new microscale mechanical-chemical (MC) model to analyze creep of salt at the microscale, accounting for coupled deformation, dynamic contacts, and chemical reaction in granular systems with realistic geometric representations. The MC model was realized by linking a new microscale mechanical code based on the numerical manifold method (NMM) to a reactive transport code named Crunch. We simulated the processes of reorganization of the salt grains, microfracturing, and pressure solution that contribute to creep at larger scales. Based on this first quantitative microscale model, we found that sharp corners of mineral grains can dominate the contact dynamics, microfracturing, and pressure solution when salt is compacted, thus governing the structural changes and porosity loss of the system. We found that pressure solution, which preferentially dissolves sharp corners and edges, can lead to relatively high porosity loss in the system, thus playing an important role in the creep of salt. Our analysis shows that the dynamic changes of salt granular systems involving grain relocation and pressure solution can occur repeatedly and continuously, thus contributing to longer-term creep of salt at larger scales.

The complex geometric features of subsurface fractures at different scales makes mesh generation challenging and/or expensive. In this paper, we make use of neural style transfer (NST), a machine learning technique, to generate mesh from rock fracture images. In this new approach, we use digital rock fractures at multiple scales that represent ‘content’ and define uniformly shaped and sized triangles to represent ‘style’. The 19-layer convolutional neural network (CNN) learns the content from the rock image, including lower-level features (such as edges and corners) and higher-level features (such as rock, fractures, or other mineral fillings), and learns the style from the triangular grids. By optimizing the cost function to achieve approximation to represent both the content and the style, numerical meshes can be generated and optimized. We utilize the NST to generate meshes for rough fractures with asperities formed in rock, a network of fractures embedded in rock, and a sand aggregate with multiple grains. Based on the examples, we show that this new NST technique can make mesh generation and optimization much more efficient by achieving a good balance between the density of the mesh and the presentation of the geometric features. Finally, we discuss future applications of this approach and perspectives of applying machine learning to bridge the gaps between numerical modeling and experiments.

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.

Abstract. High resolution gridded datasets of meteorological variables are needed in order to resolve fine-scale hydrological gradients in complex mountainous terrain. Across the United States, the highest available spatial resolution of gridded datasets of daily meteorological records is approximately 800 m. This work presents gridded datasets of daily precipitation and mean temperature for the East-Taylor subbasin (in western United States) covering a 12-year period (2008–2019) at a high spatial resolution (400 m). The datasets are generated using a downscaling framework that uses data-driven models to learn relationships between climate variables and topography. We observe that downscaled datasets of precipitation and mean temperature exhibit smoother spatial gradients compared to their coarser counterparts. Additionally, we also observe that when downscaled datasets are reaggregated to the original resolution (800 m), the mean residual error is almost zero, ensuring spatial consistency with the original data. Furthermore, the downscaled datasets are observed to be linearly related to elevation, which is consistent with the methodology underlying the original 800 m product. Finally, we validate the spatial patterns exhibited by downscaled datasets via an example use case that models lidar-derived estimates of snowpack. The presented dataset constitutes a valuable resource to resolve fine-sale hydrological gradients in the mountainous terrain of the East-Taylor subbasin, which is an important study area in the context of water security for southwestern United States and Mexico. The dataset is publicly available at https://doi.org/10.15485/1822259 (Mital et al., 2021).