Resolving hydromechanical coupling in two and three dimensions: spontaneous channelling of porous fluids owing to decompaction weakening

Abstract

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 information to this paper.

Modelling thermomechanical ice deformation using a GPU-based implicit pseudo-transient method (FastICE v1.0)

“…Velocity components are located at the cell mid-faces ( Figure 1), while shear stress components are located at the cell vertices in 2-D (e.g., Harlow and Welch, 1965). The resulting algorithms are well suited for taking advantage of modern many-core parallel accelerators, such as graphical processing units (GPUs) (Omlin, 2017;Räss et al, 2018;Duretz 165 et al, 2019;Räss et al, 2019a). Efficient parallel solvers utilising modern hardware provide a viable solution to resolve the computationally challenging coupled thermomechanical full Stokes calculations in 3-D.…”

“…(1 ≤ ν ≤ 10), and it is usually problem-dependent. This approach was successfully implemented in recent PT developments by Räss et al (2018Räss et al ( , 2019a and Duretz et al (2019).…”

“…An added challenge in full Stokes models is ice's strongly non-linear thermomechanics. Ice's 30 viscosity significantly depends on both temperature and strain-rate (Robin, 1955;Hutter, 1983;Morland, 1984), which can lead to spontaneous localisation of shear (e.g., Duretz et al, 2019;Räss et al, 2019a). Particularly challenging is the scale separation associated with localisation, which leads to micro-scale physical interaction generating meso-scale features such as thermally-activated shear zones or preferential flow paths in macro-scale ice domains.…”

“…Each operation on the GPU assigns one thread to compute the update of a given grid-point. Since on the GPU device, one core can simultaneously execute several threads, the operation set is executed on the entire computational domain almost concurrently.We tailor our numerical method to optimally exploit the massive parallelism of GPU hardware (Omlin, 2017; Räss et al, 55 2018;Duretz et al, 2019;Räss et al, 2019a). Our numerical implementation relies on an iterative and matrix-free method to solve the mechanical and thermal problems using a finite-difference discretisation on a Cartesian staggered grid.…”

“…We tailor our numerical method to optimally exploit the massive parallelism of GPU hardware (Omlin, 2017; Räss et al, 55 2018;Duretz et al, 2019;Räss et al, 2019a). Our numerical implementation relies on an iterative and matrix-free method to solve the mechanical and thermal problems using a finite-difference discretisation on a Cartesian staggered grid.…”