Dual phase‐field model for immiscible fluid flow through fractured unsaturated porous media

The immiscible multiphase fluid flow in intact and fractured porous media is relevant to numerous engineering sectors, including industrial processes, geo‐engineering, and civil engineering. These encompass crucial applications that pose significant challenges in modeling, such as carbon dioxide (CO2${\rm CO}_2$) storage in underground reservoirs, contaminant transport, and desiccation‐induced hydraulic fracturing in unsaturated soils. In this work, a macroscopic framework is developed to describe the transport of immiscible and multiphase fluids in intact and fractured porous materials. The modeling approach utilizes the embedding of two phase‐field models within the continuum Theory of Porous Media (TPM). The case of negligible capillary pressure in the macroscopic scale is associated with many numerical challenges in solving this nonlinear DAEs system, which forms the major subject of the current work. Hence, the problem of the neglected capillary pressure is defined as a target problem. To deal with this, a thermodynamically consistent modified problem is defined, which incorporates an additional artificial macroscopic capillary pressure based on the Cahn–Hilliard phase‐field approach. To increase the robustness and stability of the numerical solution, the Newton–Raphson method is replaced by a homotopy method that allows a gradual transformation of a simpler, solvable problem into the original target problem that usually converges poorly. The sharp crack topology in the fracture regions is approximated by another phase‐field method, which forms a diffusive transition zone across the crack edges of the deformable porous material. For the solid matrix, we consider a small strain assumption with a heterogeneous distribution of the permeability parameter. The inclusion of heterogeneity allows for a more realistic modeling of crack propagation and fluid‐fluid interaction. For the numerical framework, the coupled DAE system is approximated by the finite element method (FEM), implemented in the open‐source project FEniCSx. Appropriate stabilization techniques are also discussed.