A model for saturated-unsaturated flow with fractures acting as capillary barriers

High-resolution modeling of the flow dynamics in fractured soils is highly complex and computationally demanding as it requires precise geometrical description of the fractures in addition to resolving a multiphase free-flow problem inside the fractures. In this paper, we present an idealized model for saturated-unsaturated flow in fractured soils that preserves the core aspects of fractured flow dynamics using an explicit representation of the fractures. The model is based on Richards' equation in the matrix and hydrostatic equilibrium in the fractures. While the first modeling choice is standard, the latter is motivated by the difference in flow regimes between matrix and fractures, that is, the water velocity inside the fractures is considerably larger than in the soil even under saturated conditions. On matrix/fracture interfaces, the model permits water exchange between matrix and fractures only when the capillary barrier offered by the presence of air inside the fractures is overcome. Thus, depending on the wetting conditions, fractures can either act as impervious barriers or as paths for rapid water flow. Since in numerical simulations each fracture face in the computational grid is a potential seepage face, solving the resulting system of nonlinear equations is a nontrivial task. Here, we propose a general framework based on a discrete-fracture matrix approach, a finite volume discretization of the equations, and a practical iterative technique to solve the conditional flow at the interfaces. Numerical examples support the mathematical validity and the physical applicability of the model. Due to the presence of air in dry fractured soils, fractures of large aperture act as capillary barriers. If the capillary barriers are overcome, fractures become fast-flowing paths for water to travel downward. A flow model based on Richards' equation in the soil and instantaneous ponding in the fractures is proposed. The model is numerically consistent and its physical applicability is showcased in two-dimensional simulations.