Flow in fractured porous media is encountered in a broad spectrum of industrial, environmental, and engineering applications. This covers, for example, geothermal systems and hydraulic fracturing, carbon dioxide sequestration, management of karstic aquifers, and oil production. For those and other applications, modeling has become an indispensable tool for several purposes, such as understanding physical processes, developing monitoring systems, and providing predictions. Most of the existing studies in the literature deal with saturated porous media while flow in variably saturated fractured porous media is marginally investigated and the related processes are not fully understood.
The flow in variably saturated porous media is often modeled using Richards’ equation (RE) that incorporates nonlinear constitutive relations between pressure head, hydraulic conductivity, and water content. From a numerical point of view, accurate numerical solution of the RE remains a challenge, especially in the presence of formation fractures and sharp wetting fronts. Fractures typically introduce high complex geometry that requires dense meshes with fine elements to capture the high contrast in properties between the fractures and the adjacent matrix zones.
The Discrete Fracture Network (DFN) approach, in which the fractures are embedded as lower- dimensional elements, (n-1)-D, in higher-dimensional physical domain, n-D, is widely used to deal with the flow in saturated fractured porous media. To the best of our knowledge, the DFN approach has not been implemented for flow in unsaturated porous media, which is currently modeled using the traditional dual-porosity and equivalent porous media approaches. Thus the first objective of this work is to develop a new efficient model for variably saturated flow in fractured porous media based on the DFN approach.
In this model, we use Richard’s equation for both the matrix and fractures. This can be encountered in the case of sediment-ﬁlled fractures or thin high permeable layers in porous media. To address challenges related to space discretization, the mixed hybrid finite element method with a new mass lumping technique is implemented to avoid spurious oscillations. Time integration is performed using an advanced solver based on a higher-order backward differentiation with adaptivity in time stepping. The proposed solver has the privilege of selecting time step sizes within the fractures different from those in the matrix. The newly developed DFN model is verified against an industrial FE-based simulator and showed significant superiority in terms of efficiency and robustness.
Keywords: Discrete Fracture Network; Variably Saturated Flow; Mixed Hybrid Finite Element method; Mass Lumping