One of the continuing trends in geodynamics is to develop codes that are suitable to model thermo-hydromechanical processes with an increasing level of self-consistency. Developing such models is particularly challenging as some coupling terms can introduce strong non-linearities. We demonstrate that a finite difference discretization combined with pseudo transient solvers are well suited for such problems. Moreover, due to the inherent parallelism of the algorithm and thanks to the novel Julia language and the ParallelStencil.jl package, GPU implementation is not only feasible but also straightforward.

Here, we consider fluid flow in a deformable porous medium coupled to thermo-mechanical processes. We present a thermodynamically self-consistent set of governing equations, describing such processes. The governing equations consist of the conservation of mass, momentum, and energy in two phases. One phase represents the solid skeleton, which deforms in a poro-(visco-)elasto-plastic manner. The second phase represents a low viscosity fluid (water, CH, melt), percolating through the solid skeleton, that is described by Darcy’s law. A special process we will investigate is brittle failure of the matrix due to high fluid pressure (hydro-fracturing, fault reactivation, diking).

The system of equations is solved numerically, using the pseudo transient method, that is well suited to solve highly non-linear problems, as solving the global equations and iterating the non-linearities can be done at the same time. Moreover, the algorithm requires large number of local and cheap operations, which is ideal for GPU implementation. We will describe the governing equations, their numerical implementation, and show examples of numerical simulations that include mode-1 and mode-2 fractures. We apply the numerical model to study the effects of pore pressure in a siliciclastic reservoir on fault reactivation and associated integrity of cap rocks, and to provide a continuum analogue to magmatic dike emplacement.