Coupled thermal, hydraulic, mechanical, and chemical (THMC) processes play a critical role in various energy and environmental engineering applications. The recently developed Darcy-Brinkman-Biot (DBB) framework has proven effective in modeling multiphase fluid flow in deformable solids across both pore and Darcy scales. In this study, we extend the DBB framework, originally designed for isothermal conditions, to address non-isothermal problems by incorporating an energy conservation equation. The resulting solver, hybridBiotThermalInterFoam, enables the simulation of multiphase fluid flow, heat transfer, and solid deformation in hybrid systems that encompass both solid-free regions and porous media. The new solver is validated through comparisons with analytical solutions for heat transfer in one-dimensional porous media with two fluids, as well as against established heat transfer solvers like chtMultiRegionFoam. Results demonstrate that the solver accurately reproduces consistent results in alignment with analytical and numerical benchmarks. Further, a series of 2D and 3D case studies, including two-phase flow and heat transfer in solid-free and deformable porous media, highlighted the solver’s capacity to simulate complex flow dynamics and heat transport in systems involving high mobility ratios, viscous fingering, and fracture propagation. The findings emphasize the importance of incorporating thermal effects when modeling multiphase flow, especially for applications like enhanced oil recovery, soil remediation, and enhanced geothermal systems. This work establishes hybridBiotThermalInterFoam as a powerful and reliable tool for investigating multiphase, multiscale, and multiphysics processes in porous materials. It offers significant potential for advancing our understanding of THMC coupled processes in energy geotechnics, geo-engineering, and environmental systems. 

Clay minerals – layered silicate nanoparticles that constitute roughly half of the mineral mass in soils, sediments, and sedimentary rocks – are a critical building block of the Earth’s subsurface. A key aspect of these materials is their control over fluid flow, solute migration, and solid mechanics in many systems. The aggregation of clay platelets strongly impacts associated transport and mechanical properties. Unfortunately, experimental and computational characterization of clay aggregation is inhibited by the delicate water-mediated nature of the process and by the wide range of spatial scales involved, from 1-nm-thick platelets to flocs with dimensions up to micrometers or more. In this work, we utilize a new coarse-grained molecular dynamics (CGMD) model to predict the microstructure, dynamics, and rheology of hydrated smectite clay assemblages as a function of particle diameter (6 to 25 nm) and the proportion of Na vs Ca exchangeable cations. Simulated systems have dimensions up to 0.1 m and contain up to 2,000 clay particles. We monitor the assembly and growth of clay tactoids (i.e., stacks of parallel clay platelets) and aggregates (i.e., assemblages of tactoids) as a function of simulation time, shear rate, exchangeable cation type, and clay platelet diameter. Equilibrated configurations are analyzed to predict structural properties including the number of platelets per tactoid, the size distribution of clay tactoids, clay basal spacing, the distribution of counterions within clay tactoids, and the modes of clay association. Simulation results are further analyzed to predict the rheological properties of smectite gels. Our results unlock new potential to characterize and understand clay aggregation in dilute suspension on a scale of thousands of particles with explicit representation of their counterion clouds.

The aggregation of clay minerals in liquid water exemplifies colloidal self-assembly in nature. These negatively charged aluminosilicate platelets interact through multiple mechanisms with different sensitivities to particle shape, surface charge, aqueous chemistry, and interparticle distance and exhibit complex aggregation structures. Experiments have difficulty resolving the associated colloidal assemblages at the scale of individual particles. Conversely, all-atom molecular dynamics (MD) simulations provide detailed insight on clay colloidal interaction mechanisms, but they are limited to systems containing a few particles. We develop a new coarse-grained (CG) model capable of representing assemblages of hundreds of clay particles with accuracy approaching that of MD simulations, at a fraction of the computational cost. Our CG model is parameterized based on MD simulations of a pair of smectite clay particles in liquid water. A distinctive feature of our model is that it explicitly represents the electrical double layer (EDL), i.e., the cloud of charge-compensating cations that surrounds the clay particles.Our new model captures the simultaneous importance of long-range colloidal interactions (i.e., interactions consistent with simplified analytical models, already included in extant clay CG models) and short-range interactions such as ion correlation and surface and ion hydration effects. The resulting simulations correctly predict, at low solid-water ratios, the existence of ordered arrangements of parallel particles separated by water films with a thickness up to ~10 nm and, at high solid-water ratios, the coexistence of crystalline and osmotic swelling states, in agreement with experimental observations.

Bentonite clay is considered for use as a protective buffer in the isolation of high-level radioactive waste (HLRW) because of its low permeability and high swelling pressure. However, the decay of radioactive waste can increase the temperature of bentonite barriers above 100 °C, resulting in changes in coupled thermal, hydraulic, mechanical, and chemical (THMC) properties of bentonite clay. This study utilizes large-scale all-atom molecular dynamics (MD) simulations of hydrated Na-montmorillonite assemblages to predict the THMC properties of smectite at temperatures above 100 °C using an approach previously validated below 100 °C. The dehydration process of bentonite is modeled by progressively removing water molecules from the system, followed by equilibration to high dry densities. Data analysis yields various material properties of bentonite clay as functions of temperature and dry density, including hydraulic conductivity, self-diffusivity of water, heat capacity, thermal expansion coefficient, and suction. Predictions are compared with experimental results on compacted bentonite properties measured on scales of centimeters and days. Cross-validations between MD simulation results and experimental observations help bridge knowledge gaps across different scales, facilitating the use of nanoscale simulations to parameterize bentonite properties in large-scale simulators. This work advances the understanding of the THMC properties of bentonite at high temperatures and supports the evaluation of bentonite performance in isolating high-level radioactive waste.

Coupled thermal, hydraulic, mechanical, and chemical (THMC) processes, such as desiccation-driven cracking or chemically driven fluid flow, significantly impact the performance of composite materials formed by fluid-mediated nanoparticle assembly, including energy storage materials, ordinary Portland cement, bio-inorganic nanocomposites, liquid crystals, and engineered clay barriers used in the isolation of hazardous wastes. These couplings are particularly important in the isolation of high-level radioactive waste (HLRW), where heat generated by radioactive decay can drive temperature up to at least 373 K in the engineered barrier. Here, we use large-scale all-atom molecular dynamics (MD) simulations of hydrated smectite clay nanoparticle assemblages to predict the fundamental THMC properties of hydrated compacted clay over a wide range of temperatures (up to 373 K) and dry densities relevant to HLRW management. Equilibrium simulations of clay-water mixtures at different hydration levels are analyzed to quantify material properties including thermal conductivity, heat capacity, thermal expansion, suction, water and ion self-diffusivity, and hydraulic conductivity. Predictions are validated against experimental results on the properties of compacted bentonite clay. Our results demonstrate the feasibility of using atomistic-level simulations of assemblages of clay nanoparticles on scales of tens of nanometers and nanoseconds to infer the properties of compacted bentonite on scales of centimeters and days, a direct upscaling over 6 orders of magnitude in space and 15 orders of magnitude in time. 

Bentonite, a fine-grained geologic material rich in smectite clay, is considered for use in the isolation of high-level radioactive waste (HLRW) because of its low hydraulic permeability, high swelling pressure, and geochemical stability. A complicating factor in this application is that heat released by nuclear waste can trigger complex coupled thermal-hydraulic-mechanical-chemical (THMC) phenomena within the barrier. Prediction of these phenomena using large-scale simulators, which typically examine problems on scales of 10-2 to 104 m, is inhibited by insufficient knowledge of the material properties of bentonite and their dependence on temperature. Here, these properties were evaluated using replica-exchange molecular dynamics (REMD) simulations of a clay assemblage containing 27 Na-smectite nanoparticles with full atomistic-level resolution solvated using 187,131 water molecules. The simulations yielded predictions of heat capacity, thermal conductivity, thermal expansivity, hydraulic conductivity, and water and ion diffusivity at temperatures of 298 to 373 K. Results showed that temperature modulates the capacity of clay barriers to transfer heat, fluids, and chemical species to different degrees. Material properties of hydrated smectite predicted on scales of tens of nanometers and nanoseconds were consistent with the properties of bentonite measured on scales of centimeters and days.

Structural fault trapping of buoyant fluids, such as hydrocarbons and CO2, relies on fault sealing capacity. The traditional approach to quantifying the sealing capacity of fault zone materials is Shale Gouge Ratio (SGR), which implicitly homogenizes fault properties. However, a large range of fault throw and the presence of fault heterogeneity can lead to different trapping capacities for the same structure. This paper presents a stochastic study to statistically determine the possible range of CO2 column height at a normal fault in sand-shale sequences. We present the procedures of two stochastic models including the continuous shale gouge model (CSGM) and the discrete smear model (DSM). The models are applied in an example case for the High Island field in the Gulf of Mexico (GOM) basin combining borehole geophysical data and measurements from laboratory experiments. Both one-dimensional and two-dimensional cases are discussed. The results show that more than 50% of the realizations predict negligible fault sealing capacity. The CSGM method shows that CO2 column height is overall proportional to SGR and breakthrough pressure. The DSM model suggests that clay smear continuity may significantly affect fault sealing capacity. Large fault throws (> 100 m) result in smear breaches and therefore in highly variable maximum CO2 column heights, varying from 0 m to a maximum of ~95 m. Both the DSM method and the CSGM method proposed in this paper permit explaining geologic uncertainties and the resulting variability of column heights observed in the field. Our results together with the field data from the literature demonstrate the utility of combining both methods to help improve column height prediction. Uncertainty quantification of clay ductility, smear location, and fault heterogeneity is important for determining fault sealing capacity and improving reservoir risk management associated with carbon dioxide geological storage.

Commercial-scale development of CO2 geological storage necessitates robust and real-time monitoring methods to track the injected CO2 plume and provide assurance of CO2 storage. Pressure monitoring above the injection zone is a method to detect potential CO2 leaks into overlying formations. We present a generic CO2 storage model with a single injector to predict pressure changes above the caprock due to both fast hydraulic communication and partially undrained loading, the latter often neglected in reservoir simulation. The simulation used a compositional simulator coupled with geomechanics to solve the poroelastic equations in the entire storage complex. The results show that changes of pore pressure above the caprock caused by partially undrained loading reach up to ~15 kPa within ~10 days followed by a gradual decay with time. This is about 1% of the pressure increase in the injection zone. Furthermore, the pressure changes above the caprock are closely related to the advance of the CO2 plume. The results also include forward simulations considering the presence of: a fault either with high or low permeability, a poorly isolated abandoned well, a leaky injector, and a second injector. Fluid flow through high permeability paths across the caprock favors a ~one order of magnitude higher, yet more gradual pressure increase than the base case with a fully covering caprock. Pressure monitoring above the caprock is a feasible technology to track the CO2 plume, requires high precision pressure measurements, and must account for partially undrained poroelastic loading to interpret correctly measured pressure signals in the field.

The injection of fluids into a compartmentalized formation induces pore pressure buildup and may result in re-activation of sealing faults. Among other variables, the pore volume compressibility (PVC) can affect the amount of pore pressure change during injection. PVC has been traditionally measured with isotropic loading compressibility tests. However, long and thin reservoirs subjected to depletion or injection typically follow a uniaxial-strain stress path, rather than an isotropic stress path. Furthermore, injection unloads the reservoir rock by reducing effective stress, whereas depletion causes loading. This paper reports experimental measurements of the uniaxial strain unloading compressibility of Frio sand, a member of Tertiary strata in the Gulf of Mexico Basin. The uniaxial strain unloading compressibility increases non-linearly from 0.29 to 1.45 GPa-1 (2 to 10 µsip) as the mean effective stress is reduced from 26 to 5 MPa. The uniaxial strain unloading compressibility of Frio sand is about one-third of the uniaxial strain loading compressibility at comparable levels of effective stress. The uniaxial strain compressibility of Frio sand is roughly one half of the isotropic compressibility. Reservoir simulation highlights that using incorrect pore compressibility values considerably underestimates the expected increase of pore pressure in a compartmentalized formation during injection. Uniaxial strain unloading compressibility of target reservoir rocks should be accurately estimated or measured in order to prevent excessive pressure build-up in target storage  formations during injection of CO2 or any other fluid.

Faults are key components in defining fluid migration pathways and seals in sedimentary basins. The sealing capacity of faults is closely related to the petrophysical and geomechanical properties of fault gouge. Clay smear, cataclasis, and diagenesis favor a high capillary breakthrough pressure and low permeability in clastic sediments, and therefore fault gouge seal. However, significant uncertainty remains in accurately predicting the sealing capacity of faults for CO2 storage. We conducted a series of experiments, including absolute permeability, breakthrough pressure, and post-breakthrough CO2 permeability measurements on synthetic fault gouge samples, made from homogeneous mixtures of Frio sand and Anahuac shale, lithofacies of tertiary sediments in the Gulf of Mexico basin. The results show that the permeability of synthetic fault gouge decreases by about one order of magnitude with increments of 10 wt% of clay (mostly smectite) from the Anahuac shale. Independent of clay content and bulk porosity, all permeability measurements scale with the void ratio of the clay fraction. The breakthrough pressure of synthetic fault gouge increases by approximately half order of magnitude with increments of 10 wt% clay. The samples with clay content above 40 wt% reach a breakthrough pressure equivalent to a (supercritical) CO2 column height of more than 100 m. The measurements on fault gouge properties are meaningful to quantitatively evaluate fault sealing capability and migration of buoyant fluids through faults in sand-shale sequences.