Reactive Mechanics and Transport Modeling
Mengsu Hu & Carl Steefel
Overview: A new coupled reactive mechanics and transport (RMT) simulation code will be developed that simulates the pore-discontinuum scale evolution of porous rock and fractures via individual and coupled pathways. The simulation capability will be achieved by linking two LBNL codes: the numerical manifold method for continuum-discontinuum mechanics with contact alteration and the Crunch family of codes for chemical reaction and transport, including transport in nanoporous media. The codes utilize different meshes and interpolations in space and time and they are coupled with optimized routines.
The RMT code will contribute to the studies of the mechanical compaction of sedimentary carbonates by simulating settling, reorganization and dewatering of biogenic carbonate grains of realistic morphology. The RMT code will contribute to studies of chemical–mechanical compaction of sedimentary carbonates by incorporating single-interface strain-rate laws developed in the Pressure Dissolution project and testing those laws by simulating grain-scale compaction studies analyzed by X-ray nanotomography and other experimental approaches. The code will enable numerical investigation of the effects of slip on rough surfaces and of mechanical instabilities that will be used to guide subsequent experimental investigation.
Background
In this proposal, we seek to understand and predict the fundamental pore- and grain-scale chemical–mechanical processes of rock deformation, fracturing, healing and lithification. A predictive model capable of represented all fundamental processes at this length scale is essential for testing the role(s) of individual processes, their coupling, and their contribution to continuum-scale properties and emergent behavior. Although numerous numerical codes have been to developed to simulate the flow of fluids in the subsurface and the associated mechanical, chemical and thermal effects, none are fully suitable for this project. Common approaches for modeling large-scale fractured or porous media use either implicit representation of fractures or pores (Rutqvist et al. 2013; Gan and Elsworth 2016) or highly simplified explicit representations of fractures such as line segments or solid elements with different material properties (Ulven et al. 2014; Yasuhara et al 2016; Ogata et al. 2019; Evans et al. 2018; Zhang et al. 2019). Recent developments in the modeling of single fractures use explicit representations of rough rock surfaces in contact, but without consideration of the mechanics and contact alteration of the fracture (McDermott et al. 2015; Pan et al. 2016; Yasuhara and Elsworth 2006). Alternative approaches for fracture simulations use a statistical representation of roughness and asperity interactions combined with homogenization to predict the evolution of a single fracture (Ameli et al. 2014; Bond et al. 2017).
In all of these models, contacts along rough surfaces cannot be accurately captured because the geometric features (such as rough boundaries along fractures, pores and grains) are not explicitly represented, or they are approximated by spheres or rectangular grids. Chemical reactions are simplified with a set of empirical equations, such as an assumed mineral saturation state, instead of incorporating the full multicomponent reaction network. Numerical modeling of chemical–mechanical processes has never been attempted at the pore to microscopic scale where discontinuities are important. We term this approach reactive mechanics and transport (RMT) modeling at the discontinuum pore scale
Reactive mechanics and transport (RMT) modeling at the discontinuum pore scale
At discontinuum pore scales, numerical modeling of coupled chemical-mechanical processes faces two challenges. First, pore-scale models have more complicated geometric features that cannot be simplified easily, and therefore these geometric features lead to discontinuities in each physical field (i.e., flow, transport, mechanics).
Secondly, pore-scale processes are described by a different set of governing equations for flow, transport, chemical reactions and mechanics, as listed in Table 1. Darcy’s law is not sufficient to describe flow in the pores/channels, and the Navier-Stokes equation is required. For chemical reaction, transport to and from discontinuous surfaces should be considered in the mass balance instead of simply using an averaged, continuum reaction rate. For mechanics, it is not accurate to rely on the small-strain assumption and continuum constitutive laws, translational and rotational displacements need to be considered.
The most challenging additional equations involve the calculation of contact forces between multiple blocks. It is a grand challenge in the field to identify when and where contacts occur among many blocks, which is made complicated by the fact these blocks are moving, deforming, and in some cases breaking apart. In turn, motion, deformation and breakage of blocks are impacted by contact forces, which constitute a serial process. Thus, inaccurate calculation of contacts can lead to completely different overall system behavior. In order to carry out contact calculations, simplifications are often made, either assuming contact pairs are fixed, or by using simple shapes such as spheres or rectangles to approximate them (Fig. 1).
The numerical manifold method (NMM) for mechanics
The numerical manifold method (NMM) was proposed by Shi (1992) for continuum and discontinuum mechanical analysis. Subsequently, this Task Lead, developed multi-scale hydro-mechanical models and a series of computer code for fractures and porous media with rigorous geometrical and physical representations, including those at continuum (Hu et al. 2017b), discrete fracture (Hu et al. 2017a), and discontinuum pore (Hu and Rutqvist 2020) scales. At discontinuum pore scales, an algorithm was developed that rigorously incorporates contact detection, contact enforcement and open-close iteration. In each time step, contact detection is carried out to determine when and where contact occurs among elements divided by discontinuities.
As shown in Figure 2 considering the multiple blocks that are present, possible contact blocks are calculated. For every two possible contacting blocks, the contact pairs are calculated to determine on which edge or vertex contact could occur exactly. Figure 2 shows all the possibilities of contact lines (marked with thick lines) for two possible contacting angles, depending on the degrees of these two angles and their relative locations and relative orientations of their bisectors . For each contact pair, their contact state can be open (not in contact), closed (attached to each other), or sliding. Correspondingly, their contact forces are imposed in distinct ways due to different constitutive laws that apply. Within each time step, open-close iteration may be carried out several times until the enforced contacts reach convergence.
Coupling between the NMM and Crunch codes
The software packages CrunchFlow (Steefel et al, 2015) and CrunchClay (Steefel & Tournassat, 2020) offer a powerful approach to incorporate geochemistry into a microscale coupled mechanical-chemical model. Currently, a new model linking NMM to CrunchFlow is under development.
The two codes are sequentially coupled (Fig. 3). NMM calculates the displacement, strain and stress and updates the geometry each time. Then the calculated stress and updated geometry is transferred to CrunchFlow. CrunchFlow updates properties associated with the updated geometry from NMM, and uses the stress components to calculate free energy, and update the rate constants and solubilities. In turn, chemical dissolution or precipitation calculated by CrunchFlow may lead to changes in the geometry and/or mechanical properties, and those are calculated in NMM. As mechanical deformation and movement will occur much faster than chemical reactions and fluid transport, all the components will be flexibly treated as coupled or decoupled depending on the rates of processes and the time step chosen.
The challenges of linking of NMM to CrunchFlow are primarily associated with differences of meshes and interpolations in space and time. In particular, NMM uses Lagrangian mesh and allows for boundaries to across meshes, while CrunchFlow uses Eulerian mesh and the mesh needs to conforms to boundaries, which makes it challenging to capture large deformation and moving boundaries when contact occurs among interfaces and particles. New approaches are currently being explored for possible implicit representation of boundaries across meshes in CrunchFlow. As mechanical and chemical processes can occur at different time scales, differing number of time steps within each coupling iteration may be conducted in NMM and CrunchFlow until a converged solution is achieved.
Preliminary results of coupled processes induced by loading on the top surface of a single asperity are shown in Figure 4. Concentration of vertical normal stress occurs at the contact interface (Fig. 4b), and that causes strong pressure dissolution (Fig. 4d). Shear stress is significant at the four corners of asperity, which contributes to stronger pressure solution at those locations. At the free surfaces (shown clearly at Fig. 5b), mineral precipitation occurs (Fig. 5e). Such a simple example demonstrates the promise of the NMM-CrunchFlow, and this lays the basis for the proposed work.
Proposed Work
The long-term vision for this effort is to create a versatile pore- and grain-scale simulation tool capable of representing all mechanisms of brittle and plastic mineral deformation in carbonate and clay-bearing rocks, including temperature and chemical effects on the rates of those processes.
The code will be capable of numerical predictions of important questions at two spatial scales. Interface-scale simulations will investigate and predict how interfacial roughness, flow, normal and shear stresses, and chemical reaction affect interfacial properties (e.g., permeability and strength) and interface evolution through mechanisms such as pressure dissolution, asperity fracturing and healing. Mesoscale simulations will investigate and predict how interfacial properties and interface evolution contribute to continuum properties of granular rock and fractures. The models incorporated in the code will be extensively tested by experiment and will be used to design experimental studies.
Forces and Transport of Mineral Grain Boundaries
The structure and transport properties at grain and asperity contacts determine true contact area and contact stress and control the rates of solute redistribution. NMM-Crunch code will incorporate realistic descriptions of the mineral–mineral interaction forces at stressed interfaces and the coupled structure and transport properties.
The interface between rough mineral grains will be composed of asperity contacts, where the normal stress causes overlap of the EDL, and larger pores through which normal diffusion and flow can occur (Fig. 5). The thickness of the intervening aqueous film will be determined by the balance between the normal stress and the overlapping EDL disjoining pressure, which will be calculated by the continuum model developed in CrunchClay (extended version of CrunchFlow for clay) that can simulate both the Nernst-Planck equation (with explicit treatment of ions diffusing at different rates and the associated diffusion potential) and the electrical double layer (EDL).
The treatment of molecular diffusion, which regulates the transfer of mass from the stressed contact surfaces to free faces, using the Nernst-Planck equation is much more rigorous than is possible with a simple Fickian approach (Tournassat and Steefel, 2015). In addition, the water films separating stressed contacted grains may be as thin as 1-10 nanometers, in which case EDL overlap may occur. The EDL models considered in this implementation of CrunchFlow (Steefel et al, 2015) will be based on both the Poisson-Boltzmann equation and the Mean Electrostatic Potential approach (Tournassat and Steefel, 2015).
Lateral displacements caused by shear stress at mineral interfaces can be modeled using different approaches. Displacement of hydrated asperity contacts can be explicitly modeled using lubricated slip models and experimental or simulation values for the viscosity of confined aqueous films. Preliminary simulation suggest that calcite–calcite interfaces may lose water under relatively low stress, thereby the development of a friction law for anhydrous interfaces may be required. The NMM-Crunch model can incorporate any friction rate law from the literature or developed in this proposal.
Settling and Compaction of Clays and Carbonates
Because NMM has rigorous treatment for multi-body systems, it can be used to study reorganization of grain aggregates due to settling in the first regime of Carbonate Compaction. Thus, the code will be used to simulate the settling of inorganic and biogenic carbonates in order to establish the limits of compaction observed in shallow sediments and test the predictions of molecular simulation concerning the sliding and shearing along rough surfaces. The geometry will be explicitly represented and more detailed interfacial constitutive laws for slip will be incorporated as needed.
The code can also be applied to the settling of highly anisotropic minerals such as clays.
Pressure Dissolution at Ideal and Rough Interfaces
Fully coupled chemical–mechanical simulation capability will be developed through NMM-Crunch simulations of Pressure Dissolution at stressed mineral surfaces through subtasks that will be tightly integrated with experiment and simulation. For explicit simulation geometries, including a single rough interface and a granular rock, NMM will calculate the local stress tensor throughout the solid. Thermodynamic modeling will be used to determine the departure from local equilibrium and a chemical kinetic rate law will be used to determine the dissolution rate as a function of stress, surface speciation and local aqueous solution chemistry. A diffusive transport step will redistribute solutes between compressed thin films and connected pores. Where transient aqueous conditions exceed mineral saturation states, models for growth or nucleation will be implemented. There are numerous questions about the controls on each of these processes and this NMM–Crunch RMT model will enable explicit numerical investigate of the most appropriate models informed by experiment and simulation.
Extensions and Applications
Complementary numerical studies will be performed to further evaluate the chemical–mechanical processes likely to be important in rock creep. In particular, different loading scenarios will be explored to study instability behaviors such as Asaro-Tiller-Grinfeld instability.
References
Ameli P, Elkhoury JE, Morris JP, & Detwiler RL (2014) Fracture permeability alteration due to chemical and mechanical processes: A coupled high-resolution model. Rock Mechanics and Rock Engineering 47(5):1563-1573.
Evans O, Spiegelman M, Kelemen PB (2018) A poroelastic model of serpentinization: Exploring the interplay between rheology, surface energy, reaction, and fluid flow. Journal of Geophysical Research: Solid Earth 123: 8653–8675.
Evans O, Spiegelman M, Kelemen PB (2019) Phase-field modeling of reaction-driven cracking: Determining conditions for extensive olivine serpentinization. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1029/2019JB018614
Gan, Q., & Elsworth, D. (2016). A continuum model for coupled stress and fluid flow in discrete fracture networks. Geomechanics and Geophysics for Geo-Energy and Geo-Resources, 2(1), 43-61.
Hu, M., & Rutqvist, J. (2020). Finite volume modeling of coupled thermo-hydro-mechanical processes with application to brine migration in salt. Computational Geosciences.
Hu M, Rutqvist J, Wang Y (2017a) A numerical manifold method model for analyzing fully coupled hydro-mechanical processes in porous rock masses with discrete fractures. Advances in water resources.102: 111-126.
Hu M, Wang Y, Rutqvist J (2017b) Fully coupled hydro-mechanical numerical manifold modeling of porous rock with dominant fractures. Acta Geotechnica 12(2): 231-252.
McDermott C, Bond A, Harris AF, Chittenden N, & Thatcher K (2015) Application of hybrid numerical and analytical solutions for the simulation of coupled thermal, hydraulic, mechanical and chemical processes during fluid flow through a fractured rock. Environmental Earth Sciences 74(12):7837-7854.
Ogata S, Yasuhara H, Aoyagi K, Kishida K (2019) Coupled THMC analysis for predicting hydro-mechanical evolution in siliceous mudstone. 53rd U.S. Rock Mechanics/Geomechanics Symposium, 23-26 June, New York City, New York.
Rutqvist, J., Leung, C., Hoch, A., Wang, Y., & Wang, Z. (2013). Linked multicontinuum and crack tensor approach for modeling of coupled geomechanics, fluid flow and transport in fractured rock. Journal of Rock Mechanics and Geotechnical Engineering, 5(1), 18-31.
Shi G., 1992. Manifold method of material analysis. Transaction of the 9th Army Conference on Applied Mathematics and Computing. U.S. Army Research Office.
Steefel, C. I., Appelo, C. A. J., Arora, B., Jacques, D., Kalbacher, T., Kolditz, O., . . . Yeh, G. T. (2015). Reactive transport codes for subsurface environmental simulation. Computational Geosciences, 19(3), 445-478.
Tournassat, C., & Steefel, C. I. (2015). Ionic Transport in Nano-Porous Clays with Consideration of Electrostatic Effects. Reviews in Mineralogy and Geochemistry, 80(1), 287-329.
Ulven, O. I., Jamtveit, B., & Malthe-Sørenssen, A. (2014). Reaction-driven fracturing of porous rock. Journal of Geophysical Research: Solid Earth, 119(10), 7473-7486.
Yasuhara, H., & Elsworth, D. (2008). Compaction of a Rock Fracture Moderated by Competing Roles of Stress Corrosion and Pressure Solution. Pure and Applied Geophysics, 165, 1289-1306.
Yasuhara H, Kinoshita N, Ogata S, Cheon DS, Kishida K (2016) Coupled thermo-hydro-mechanical-chemical modeling by incorporating pressure solution for estimating the evolution of rock permeability. Int. J. Rock Mech. Min. Sci. 86: 104-114.