A hybrid upwind scheme for two-phase flow in fractured porous media

TL;DR

This work tackles the computational challenge of simulating two-phase flow in fractured porous media where fractures form a mixed-dimensional network. It extends the hybrid upwind (HU) discretization to a mixed-dimensional framework, enabling a more robust Newton solver by representing the total flux and gravity-driven deviations instead of treating phases separately as in phase-potential upwinding (PPU). Across 2D and 3D cases with single fractures and complex networks, HU consistently reduces Newton iterations and improves convergence, particularly in strong gravity or high-contrast permeability regimes. The approach is implemented in PorePy and demonstrates potential for enabling larger, more accurate simulations of fractured systems, with future work including additional phases, capillary effects, and gravity-consistent transmissibilities.

Abstract

Simulating the flow of two fluid phases in porous media is a challenging task, especially when fractures are included in the simulation. Fractures may have highly heterogeneous properties compared to the surrounding rock matrix, significantly affecting fluid flow, and at the same time hydraulic aperture that are much smaller than any other characteristic sizes in the domain. Generally, flow simulators face difficulties with counter-current flow, generated by gravity and pressure gradients, which hinders the convergence of non-linear solvers (Newton). In this work, we model the fracture geometry with a mixed-dimensional discrete fracture network, thus lightening the computational burden associated to an equi-dimensional representation. We address the issue of counter-current flows with appropriate spatial discretization of the advective fluid fluxes, with the aim of improving the convergence speed of the non-linear solver. In particular, the extension of the hybrid upwinding to the mixed-dimensional framework, with the use of a phase potential upstreaming at the interfaces of subdomains. We test the method across several cases with different flow regimes and fracture network geometry. Results show robustness of the chosen discretization and a consistent improvements, in terms of Newton iterations, compared to use the phase potential upstreaming everywhere.

Figures (14)

• Figure 1: Mixed-dimensional domain. (a) Domains $\Omega_1$ and $\Omega_2$ are connected through the mortar interface $\Gamma_j$. (b) Example of three branches, in 3D fracture intersections and in 2D fractures, that intersect to a point $\Omega_4$.
• Figure 2: Grid elements of a 1D domain, $\Omega_1$ facing a 2D domain, $\Omega_2$, and mortar, $\Gamma$. Any cell is denoted by $\mathfrak{c}$, while the faces, i.e., the boundary of $\mathfrak{c}$ are $\mathfrak{f}$.
• Figure 3: Case 1. Domains and grids for the example in \ref{['sec:case_1']}. The fracture is located in three different positions, (a) horizontally with tips touching the borders, (b) vertically with one tip touching the border and the other immersed in the matrix, (c) with an angle and the tips touching the borders. Different type of grids are used: (a) structured with quadrilateral elements, (b) unstructured with triangular elements, (c) unstructured non-conforming at the fracture interface with triangular elements.
• Figure 4: Case 1. Characteristics regarding the iterative method.
• Figure 5: Case 1. Horizontal. Saturation profile along a vertical line. The left panel illustrates $S_0$ at time $t = 6.8$, during the countercurrent flow of the two phases. A more diffusive trend is observable for the HU scheme than the PPU. In the right panel, showing the stationary solution at $t = 20$, the lines nearly coincide. No visible differences in the numerical diffusion are appreciable at the discontinuity since the interface fluxes are discretized with the same scheme.