A Three‐Field Formulation for Two‐Phase Flow in Geodynamic Modeling: Toward 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, with solid velocity (vs), total pressure (Pt), and fluid pressure (Pf) as unknowns, and present a robust finite‐element implementation, which can be used to solve problems in which domains of both zero porosity and non‐zero porosity are present. The reformulated equations include regularization to avoid singularities and exactly recover to the standard single‐phase incompressible Stokes problem at zero porosity. We verify the correctness of our implementation using the method of manufactured solutions and analytic solutions and demonstrate that we can obtain the 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 can be used for a large range of applications in various geodynamic settings. Plain Language Summary: The Earth's interior is hot and consists of highly viscous rocks that are deformable over geological time scales. In some regions, such as mid‐ocean ridges or hotspots, mantle temperature is high enough such that mantle rocks become partially molten, forming low‐viscosity magma. Migration of magma toward Earth's surface is a process of a low‐viscosity fluid flowing in a deformable high‐viscosity matrix, which is a two‐phase system. Investigating magma migration requires to solve the fully coupled two‐phase system. However, at lower temperatures magma disappears because of solidification, such that the system reduces to a single (solid) phase, which is mathematically problematic because the linear system for the coupled two‐phase problem is singular and under‐determined in this case. In this work, we propose a new formulation for the coupled two‐phase system, which works for both single‐phase and two‐phase cases. The new formulation is implemented and benchmarked against a set of analytical solutions. Using selected example experiments, which have well‐known solutions and contain zero‐porosity domains, we demonstrate that the model works consistently and robustly for a large range of geodynamic settings. Key Points: We propose a new three‐field two‐phase formulation that works consistently in the zero‐porosity single‐phase limitBenchmarks demonstrate that our implementation on a Q1‐P0 element has the expected spatial and temporal convergence rateThe approach is robust and suitable to resolve geological problems in various settings including vanishing porosity [ABSTRACT FROM AUTHOR]