This paper presents a finite volume discretization of two-phase Darcy flows in discrete fracture networks taking into account the mass exchange between the matrix and the fracture. We consider the asymptotic model for which the fractures are represented as interfaces of codimension one immersed in the matrix domain, leading to the so called hybrid dimensional Darcy flow model. The pressures at the interfaces between the matrix and the fracture network are continuous corresponding to a ratio between the normal permeability of the fracture and the width of the fracture assumed to be large compared with the ratio between the permeability of the matrix and the size of the domain. The discretization is an extension of the Vertex Approximate Gradient (VAG) scheme to the case of hybrid dimensional Darcy flow models. Compared with Control Volume Finite Element (CVFE) approaches, the VAG scheme has the advantage to avoid the mixing of the fracture and matrix rocktypes at the interfaces between the matrix and the fractures, while keeping the low cost of a nodal discretization on unstructured meshes. The convergence of the scheme is proved under the assumption that the relative permeabilities are bounded from below by a strictly positive constant. This assumption is needed in the convergence proof in order to take into account discontinuous capillary pressures in particular at the matrix fracture interfaces. The efficiency of our approach compared with CVFE discretizations is shown on two numerical examples of fracture networks in 2D and 3D.

Flow in fractured porous media occurs in the earth's subsurface, in biological tissues, and in man-made materials. Fractures have a dominating influence on flow processes, and the last decade has seen an extensive development of models and numerical methods that explicitly account for their presence. To support these developments, four benchmark cases for single-phase flow in three-dimensional fractured porous media are presented. The cases are specifically designed to test the methods' capabilities in handling various complexities common to the geometrical structures of fracture networks. Based on an open call for participation, results obtained with 17 numerical methods were collected. This paper presents the underlying mathematical model, an overview of the features of the participating numerical methods, and their performance in solving the benchmark cases.