Abstract The interplay between compaction‐driven fluid flow and plastic yielding within porous media is investigated through numerical modeling. We establish a framework for understanding the dynamics of fluid flow in deforming porous materials that corresponds to the equations describing solitary porosity wave propagation. A concise derivation of the coupled fluid flow and poro‐viscoelastoplastic matrix behavior is presented, revealing a connection to Biot's equations of poroelasticity and Gassmann's theory in the elastic limit. Our findings demonstrate that fluid overpressure resulting from channelized fluid flow initiates the formation of new shear zones. Through three‐dimensional simulations, we observe that the newly formed shear zones exhibit a parabolic shape. Furthermore, plasticity exerts a significant influence on both the velocity of fluid flow and the shape of fluid channels. Importantly, our study highlights the potential of spontaneous channeling of porous fluids to trigger seismic events by activating both new and pre‐existing faults.
<p>Developing new numerical reactive transport models is essential for predicting and describing natural and technogenic petroleum and geological processes at different scales. Examples of such processes are pore fluid migration in subduction zones, causing seismic and volcanic activity, chemical and thermal enhanced oil recovery activities, etc. New numerical reactive transport models must be validated against analytical or semi-analytical solutions to ensure its correct numerical implementation. In this study, we construct thermo-hydro-chemo-mechanical model which takes into account multi-phase fluid flow in porous matrix associated with inter- and intra-phase chemical reactions with significant temperature and volume effect and treats porosity and permeability evolution. All equations are derived from basic principles of conservation of mass, energy, and momentum and the thermodynamic admissibility of all equations is verified. We solve the proposed system of equations both with a finite difference approach on a staggered grid and characteristic-based Lax-Friedrichs different order schemes to treat the disintegration of discontinuities.&#160; Resolving the problem of large discrepancies during the time evolution of coupled physical processes is challenging. For that, we use pseudo-iterations which force slow modes to attenuate quickly. Furthermore, we perform dimensionless analysis of the proposed model which allows us to detect proper dimensionally independent, dimensionally dependent and non-dimensional parameters. A new semi-analytical is derived which is based on a relaxation method of defining the stationary solution of system of partial differential equations, so detection specific regimes for reaction front propagation are possible. As a result, reaction front velocity dependence on Peclet, Damkohler and Lewis nondimensional parameters is obtained.</p>
Abstract Quantifying natural processes that shape our planet is a key to understanding the geological observations. Many phenomena in the Earth are not in thermodynamic equilibrium. Cooling of the Earth, mantle convection, mountain building are examples of dynamic processes that evolve in time and space and are driven by gradients. During those irreversible processes, entropy is produced. In petrology, several thermodynamic approaches have been suggested to quantify systems under chemical and mechanical gradients. Yet, their thermodynamic admissibility has not been investigated in detail. Here, we focus on a fundamental, though not yet unequivocally answered, question: which thermodynamic formulation for petrological systems under gradients is appropriate—mass or molar? We provide a comparison of both thermodynamic formulations for chemical diffusion flux, applying the positive entropy production principle as a necessary admissibility condition. Furthermore, we show that the inappropriate solution has dramatic consequences for understanding the key processes in petrology, such as chemical diffusion in the presence of pressure gradients.
<p>Many geodynamic processes are coupled. For example, in the partially molten mantle, the solid and molten mantle phases interact chemically during porous melt flow. For such two-phase reactive melt migration, solid and melt densities are functions of temperature, pressure, and chemical composition. Numerical models of such coupled physical-chemical systems require special treatment of the various couplings and concise numerical implementation. We elaborate a 2-D thermo-hydro-mechanical-chemical (THMC) numerical model for melt migration by porosity waves coupled to chemical reactions (Bessat et. al., 2021). We consider a simple ternary chemical system of forsterite-fayalite-silica to model melt migration within partially molten peridotite around the lithosphere-asthenosphere boundary. Our THMC model can simulate porosity waves of different shapes depending on the ratio of shear to bulk viscosity and the ratio of decompaction to compaction bulk viscosity. For an initial circular (blob-like) porosity perturbation, having a 2-D Gaussian shape, the geometry of the propagating reactive porosity wave remains blob-like if all viscosities are similar. If the decompaction bulk viscosity is smaller than the compaction bulk viscosity, so-called decompaction weakening, then the propagating porosity wave evolves into a channelized form. Our simulations quantify the variation from a blob-like to a channel-like porosity wave as a function of the viscosity ratios. We describe the 2-D THMC numerical algorithm which is based on the pseudo-transient finite difference method. Furthermore, we quantify the impact of channelization on the chemical differentiation during melt flow. Particularly, we quantify the evolution of the total silica concentration during melt migration as a function of the degree of channelization.</p><p>References</p><p>Bessat, A., Pilet, S., Podladchikov, Y. Y., & Schmalholz, S. M. (2022). Melt migration and chemical differentiation by reactive porosity waves. Geochemistry, Geophysics, Geosystems. In press. &#160;</p>
SUMMARY The efficient and accurate numerical modelling of Biot’s equations of poroelasticity requires the knowledge of the exact stability conditions for a given set of input parameters. Up to now, a numerical stability analysis of the discretized elastodynamic Biot’s equations has been performed only for a few numerical schemes. We perform the von Neumann stability analysis of the discretized Biot’s equations. We use an explicit scheme for the wave propagation and different implicit and explicit schemes for Darcy’s flux. We derive the exact stability conditions for all the considered schemes. The obtained stability conditions for the discretized Biot’s equations were verified numerically in one-, two- and three-dimensions. Additionally, we present von Neumann stability analysis of the discretized linear damped wave equation considering different implicit and explicit schemes. We provide both the Matlab and symbolic Maple routines for the full reproducibility of the presented results. The routines can be used to obtain exact stability conditions for any given set of input material and numerical parameters.
Summary High performance computing plays an important role in Earth sciences, especially in forward modeling of physical processes and applications to seismic imaging. We here solve anisotropic elastodynamic Biot's equations in the time domain using the finite volume method on a rectangular three-dimensional grid. The resolution of the model domain is more than 1 billion grid cells. We rely on the CUDA-C language and MPI to execute the code on multiple graphical processing units. We implemented overlap of computation and MPI communication in order to achieve a weak scaling parallel efficiency of 98 percent among 8 GPUs. Our results show that high-resolution 3D forward simulations are tractable nowadays and that the speed-up, compared to traditional CPU-based approaches, is significant.
Abstract. This study explores stress drops and earthquake nucleation within the simplest elasto-plastic media using two-dimensional simulations, emphasizing the critical role of temporal and spatial resolutions in accurately capturing stress evolution and strain fields during seismic cycles. Our analysis reveals that stress drops, triggered by plastic deformation once local stresses reach the yield criteria, reflect fault rupture mechanics, where accumulated strain energy is released suddenly, simulating earthquake behavior. Finer temporal discretization leads to sharper stress drops and lower minimum stress values, while finer spatial grids provide more detailed representations of strain localization and stress redistribution. Our analysis reveals that displacement accumulates gradually during interseismic periods and intensifies during major stress drops, reflecting natural earthquake cycles. Furthermore, the initial wave field patterns during earthquake nucleation are complex, with high-amplitude shear components. The histogram of stress drop amplitudes shows a non-Gaussian distribution, characterized by a sharp peak followed by a gradual decay, where small stress drops are more frequent, but large stress drops still occur with significant probability. This "solid turbulence" behavior suggests that stress is redistributed across scales, with implications for understanding the variability of seismic event magnitudes. Our results demonstrate that high-resolution elasto-plastic models can reproduce key features of earthquake nucleation and stress drop behavior without relying on complex frictional laws or velocity-dependent weakening mechanisms. These findings emphasize the necessity of incorporating plasticity into models of fault slip to better understand the mechanisms governing fault weakening and rupture. Furthermore, our work suggests that extending these models to three-dimensional fault systems and accounting for material heterogeneity and fluid interactions could provide deeper insights into seismic hazard assessment and earthquake mechanics.
<p>Preservation of mechanically-controlled microstructures can help us to unravel the long-term stress state in geological materials. To better understand the stress state in such a microstructure, we need to quantify the processes in a coupled, chemo-mechanical, point of view. One of such a mechanically-controlled microstructure is oscillatory zoning in high-temperature metamorphic rocks. The presented example is a sharp zoned plagioclase of 150 x 200 &#181;m size with thin compositional lamellae of 1-10 &#181;m alternating from the core towards the rim. This microstructure is interpreted to be mechanically-controlled, since conventional diffusion failed to preserve the observed microstructure within timescales that would be reasonable from a regional geology point of view. In contrast, considering that chemical diffusion is coupled to mechanical deformation the observed zoning can be maintained over the geologically-relevant timescales.</p><p>Despite of the recent valuable progress in our understanding of these microstructures, the mechanisms controlling its evolution from slowly cooled rocks are still not complete. Here, we numerically investigate the coupled, chemo-mechanical, effect that generates oscillatory zones mechanically maintained over geologically relevant timescales. We test the possibility of modelling oscillatory zoning in minerals that is similar to the exsolution process. We apply a classical Cahn-Hilliard-type equation where we add more complexity considering the impact of deformation during the process.</p>
Widely accepted model of Fickian linear diffusion of inert or trace-like elements is restricted to ideal solution models of components with equal molar mass. Simultaneous diffusion of multiple concentrations without mechanical stresses is well-described by the classical Maxwell-Stefan model, which is limited to the use of concentration gradients. Quantitative predictions of concentrations evolution in real mixtures should be treated instead by modified Maxwell-Stefan closure relations, which result in a correct equilibrium limit due to the use of the chemical potential gradients instead of concentration gradients. There is no linearity and tensorial homogeneity assumptions on flux-force relationships of classical irreversible thermodynamics. Coupling the multicomponent diffusion to mechanics results in pressure gradients that contribute to the Gibbs-Duhem relationship. Note, it was demonstrated that current models used for describing chemical diffusion in presence of stress gradient don’t remain invariance with respect to the choice of units, such as mole and mass, and the thermodynamic admissibility is doubted [1].We develop a new thermodynamically admissible model for multicomponent diffusion in viscously deformable rocks. Thermodynamical admissibility of this model ensures non-negative entropy production, while maintaining invariance with respect to the choice of units and reference frame. We demonstrate the correct Fickian limit and equilibrium limit with zero gradients of chemical potentials of individual components instead of concentration gradients in classical Maxwell-Stefan model. The model satisfies conventional Gibbs-Duhem and Maxwell relationships under pressure gradients and represents the natural coupling to the viscous multi-phase models featuring spontaneous flow localization.For numerical purposes, we develop the optimal pseudo-transient scheme for diffusion fluxes coupled to viscoelastic bulk deformation. This new effective damping techniques are compared to analytical solutions. The developed model is applied for radial garnet growth with multicomponent diffusion under pressure gradient, hydration porosity waves and melt transport in the Earth’s crust.1. Tajčmanová, L., Podladchikov, Y., Moulas, E., & Khakimova, L. (2021). The choice of a thermodynamic formulation dramatically affects modelled chemical zoning in minerals. Scientific reports, 11(1), 1-9.