ALBERT

Description: We developed a numerical thermodynamics laboratory called “Thermolab” to study the effects of the thermodynamic behavior of nonideal solution models on reactive transport processes in open systems. The equations of the state of internally consistent thermodynamic data sets are implemented in MATLAB functions and form the basis for calculating Gibbs energy. A linear algebraic approach is used in Thermolab to compute Gibbs energy of mixing for multicomponent phases to study the impact of the nonideality of solution models on transport processes. The Gibbs energies are benchmarked with experimental data, phase diagrams, and other thermodynamic software. Constrained Gibbs minimization is exemplified with MATLAB codes and iterative refinement of composition of mixtures may be used to increase precision and accuracy. All needed transport variables such as densities, phase compositions, and chemical potentials are obtained from Gibbs energy of the stable phases after the minimization in Thermolab. We demonstrate the use of precomputed local equilibrium data obtained with Thermolab in reactive transport models. In reactive fluid flow the shape and the velocity of the reaction front vary depending on the nonlinearity of the partitioning of a component in fluid and solid. We argue that nonideality of solution models has to be taken into account and further explored in reactive transport models. Thermolab Gibbs energies can be used in Cahn‐Hilliard models for nonlinear diffusion and phase growth. This presents a transient process toward equilibrium and avoids computational problems arising during precomputing of equilibrium data.

Description: Plain Language Summary: The behavior of Earth materials, rocks, minerals, melts, fluids, and gases is important to predict physical processes in the Earth with computer models. The purpose of this is to study how the changes of variables such as fluid and solid composition influence the diffusion, fluid flow, and reaction in rocks. Here, we present a set of computer codes, called Thermolab, to calculate important physical properties such as density and chemical composition of solids, fluids, and melts in chemical equilibrium. The calculations are based on the Gibbs energy that exists for every material. We use computer codes, written in MATLAB/OCTAVE language, to show how this Gibbs energy is calculated and used to compute chemical equilibrium and find the physical properties such as density and chemical composition. We discuss techniques for accurate calculation of chemical equilibrium and physical properties in real rocks. Finally, we use Thermolab to formulate a computer model of fluids reacting with rocks. We find that chemical composition of the fluid and rock strongly affects the speed and shape of the boundary between reacted and unreacted rock. Thermolab can be used in phase growth models to investigate the way in which rocks develop toward equilibrium.

Description: Key Points: Thermolab: a set of MATLAB codes is presented to perform equilibrium and nonequilibrium thermodynamic calculations. Local thermodynamic equilibrium is used to study effects of nonideality of solution models on nonlinear transport processes. Nonlinear diffusion processes are investigated with Thermolab providing a transient natural physical process toward equilibrium.

Description: Deutsche Forschungsgemeinschaft (DFG) http://dx.doi.org/10.13039/501100001659

Description: Russian Ministry of Science and Higher Education

Description: https://hansjcv.github.io/Thermolab/

Description: https://doi.org/10.5281/zenodo.6383253

Description: 〈span〉〈div〉SUMMARY〈/div〉Fingering, veining, channelling and focussing of porous fluids are widely observed phenomena in the Earth’s interior, driving a range of geo-processes across all scales. While observations suggest fairly localized flow patterns induced by fractures, the classical Darcian model predicts diffusive behaviour that leads to never-ending spreading and delocalization. We here investigate an alternative physical mechanism without the need to involve fractures. Decompaction weakening leads to the formation and propagation of localized flow-pathways in fluid-saturated porous media. We numerically solve the coupled equations using high-resolution 2-D and 3-D numerical modelling to predict non-linear porous flow in a non-linearly viscously deforming matrix. We show that high-porosity channels may be a dynamic and natural outcome of sufficiently resolved hydromechanical coupling and decompaction weakening. We propose an efficient solution strategy that involves an iterative pseudo-transient numerical method to solve the coupled system of equations in a matrix-free fashion. We discuss benefits and limitations of this approach that performs optimally on hardware accelerators such as graphical processing units and is well-suited for supercomputing. We benchmark the pseudo-transient routines against commonly used direct-iterative solving strategies and show convergence towards identical results. Furthermore, we use the fast solver to systematically study in 2-D the high-porosity channel propagation velocity as a function of bulk and shear viscosity ratios and report discrepancy between 2-D and 3-D configurations. We conclude that the fluid-flow rate in the channels is up to three orders of magnitude higher than expected by pure Darcian flow regimes and show that the high-porosity channels occurrence remains with strain rate dependant shear viscosity. We provide both the 2-D MATLAB-based direct-iterative and pseudo-transient routines for full reproducibility of the presented results and suggest our model configuration as a key benchmark case to validate the implementation of hydromechanical coupling in 2-D and 3-D numerical codes. The routines are available from Bitbucket and the Swiss Geocomputing Centre website, and are also supporting informationsupporting information to this paper.〈/span〉

Description: 〈span〉〈div〉Summary〈/div〉Fingering, veining, channelling and focussing of porous fluids are widely observed phenomena in the Earth’s interior, driving a range of geo-processes across all scales. While observations suggest fairly localised flow patterns induced by fractures, the classical Darcian model predicts diffusive behaviour that leads to never-ending spreading and delocalisation. We here investigate an alternative physical mechanism without the need to involve fractures. Decompaction weakening leads to the formation and propagation of localised flow-pathways in fluid-saturated porous media. We numerically solve the coupled equations using high-resolution 2-D and 3-D numerical modelling to predict nonlinear porous flow in a nonlinearly viscously deforming matrix. We show that high-porosity channels may be a dynamic and natural outcome of sufficiently resolved hydro-mechanical coupling and decompaction weakening. We propose an efficient solution strategy that involves an iterative pseudo-transient numerical method to solve the coupled system of equations in a matrix-free fashion. We discuss benefits and limitations of this approach that performs optimally on hardware accelerators such as graphical processing units and is well-suited for supercomputing. We benchmark the pseudo-transient routines against commonly used direct-iterative solving strategies and show convergence towards identical results. Furthermore, we employ the fast solver to systematically study in 2-D the high-porosity channel propagation velocity as a function of bulk and shear viscosity ratios and report discrepancy between 2-D and 3-D configurations. We conclude that the fluid flow rate in the channels is up to three orders of magnitude higher than expected by pure Darcian flow regimes and show that the high-porosity channels occurrence remains with strain rate dependant shear viscosity. We provide both the two-dimensional MATLAB-based direct-iterative and pseudo-transient routines for full reproducibility of the presented results and suggest our model configuration as a key benchmark case to validate the implementation of hydro-mechanical coupling in 2-D and 3-D numerical codes. The routines are available from Bitbucket and the Swiss Geocomputing Centre website, and are also supplementary material to this article.〈/span〉

Description: We present analytical derivations and 2-D numerical simulations that quantify magnitudes of deviatoric stress and tectonic overpressure (i.e. difference between the pressure, or mean stress, and the lithostatic pressure) by relating them to lateral variations in the gravitational potential energy ( GPE ). These predictions of tectonic overpressure and deviatoric stress associated with GPE differences are independent of rock rheology (e.g. viscous or elastic) and rock strength. We consider a simple situation with lowlands and mountains (plateau). We use a numerical two-layer model consisting of a crust with higher Newtonian viscosity than that in the mantle, and also a three-layer model in which the two-layer lithosphere overlies a much less viscous asthenosphere. Our results (1) explain why estimates for the magnitude of stresses in Tibet, previously published by different authors, vary by a factor of two, (2) are applied to test the validity of the thin sheet approximation, (3) show that the magnitude of the depth-integrated tectonic overpressure is equal to the magnitude of the depth-integrated deviatoric stress if depth-integrated shear stresses on vertical and horizontal planes within the lithosphere are negligible (the thin sheet approximation) and (4) show that under thin sheet approximation tectonic overpressure is required to build and support continental plateaus, such as in Tibet or in the Andes, even if the topography and the crustal root are in isostatic equilibrium. Under thin sheet approximation, the magnitude of the depth-integrated tectonic overpressure is equal to the depth-integrated horizontal deviatoric stress, and both are approximately 3.5 x 10 12 N m –1 for Tibet. The horizontal driving force per unit length related to lateral GPE variations around Tibet is composed of the sum of both tectonic overpressure and deviatoric stress, and is approximately 7 x 10 12 N m –1 . This magnitude exceeds previously published estimates for the force per unit length required to fold the Indo-Australian Plate south of India, and hence the uplift of the Tibetan plateau could have folded the Indian Plate. We also discuss the mechanical conditions that are necessary to achieve isostasy, for which the lithostatic pressure is constant at a certain depth. The results show that tectonic overpressure can exist at a certain depth even if all deviatoric stresses are zero at this depth, because this tectonic overpressure is related to horizontal gradients of vertical shear stresses integrated across the entire depth of the lithosphere. The magnitude of the depth-integrated tectonic overpressure of 3.5 x 10 12 N m –1 implies that the pressure estimated from observed mineral assemblages in crustal rocks is likely significantly different from the lithostatic pressure, and pressure recorded by crustal rocks is not directly related to depth. In case of significant weakening of the entire lithosphere by any mechanism our analytical and numerical studies provide a simple estimation of tectonic overpressure via variations in GPE .

Description: 〈span〉〈div〉SUMMARY〈/div〉Numerous geological processes are governed by thermal and mechanical interactions. In particular, tectonic processes such as ductile strain localisation can be induced by the intrinsic coupling that exists between deformation, energy and rheology. To investigate this thermo-mechanical feedback, we have designed 2-D codes that are based on an implicit finite-difference discretisation. The direct-iterative method relies on a classical Newton iteration cycle and requires assembly of sparse matrices, while the pseudo-transient method uses pseudo-time integration and is matrix-free. We show that both methods are able to capture thermo-mechanical instabilities when applied to model thermally activated shear localisation; they exhibit similar temporal evolution and deliver coherent results both in terms of non-linear accuracy and conservativeness. The pseudo-transient method is an attractive alternative, since it can deliver similar accuracy to a standard direct-iterative method but is based on a much simpler algorithm and enables high-resolution simulations in 3-D. We systematically investigate the dimensionless parameters controlling 2-D shear localisation and model shear zone propagation in 3-D using the pseudo-transient method. Code examples based on the pseudo-transient and direct-iterative methods are provided in the supplementary material and are available from Bitbucket and the Swiss Geocomputing Centre website.〈/span〉