Three-field finite-element modelling of coupled two-phase flow for geological problems: Towards the zero-porosity limit

Two-phase flow, a system where Stokes flow and Darcy flow are coupled, is of great importance in the Earth's interior, such as in subduction zones, mid-ocean ridges, and hotspots. However, it remains challenging to solve the two-phase equations accurately in the zero-porosity limit, for example when melt is fully frozen below solidus temperature. Here we propose a new three-field formulation of the two-phase system and present a robust finite-element implementation, which can successfully solve for the system where zero and non-zero porosity domains are both present. The reformulated equations, with solid velocity (v_s), total pressure (P_t), and fluid pressure (P_f) as unknowns, include penalty and regularization to avoid singularities, which exactly recover to the standard single-phase Stokes with penalty at zero porosity. The new formulation is implemented using a 2-D finite-element discretization with Q1P0Q1 elements. We demonstrate the correctness of our implementation based on benchmarks against analytical solutions, which gives expected convergence rates in both space and time. Example experiments, such as self-compaction, falling block, and mid-ocean ridge spreading, show that this formulation can robustly resolve zero- and non-zero-porosity domains simultaneously, and be used for a large range of applications in various geodynamic settings.