In this work, an optimization based approach presented in [5, 6, 7] for Discrete Fracture Network simulations is coupled with the Virtual Element Method (VEM) for the space discretization of the underlying Darcy law. The great flexibility of the VEM in handling rather general polygonal elements allows, in a natural way, for an effective description of irregular solutions starting from an arbitrary triangulation, which is built independently of the mesh on other fractures. Only partial conformity is in fact obtained with this approach. Numerical results performed on several DFN configurations confirm the viability and efficiency of the resulting method.

We investigate a new numerical approach for the computation of the threedimensional flow in a discrete fracture network that does not require a conforming discretization of partial differential equations on complex three-dimensional systems of planar fractures. The discretization within each fracture is performed independently of the discretization of the other fractures and of their intersections. An independent meshing process within each fracture is a very important issue for practical large-scale simulations, making mesh generation easier. Some numerical simulations are given to show the viability of the method. The resulting approach can be naturally parallelized for dealing with systems with a huge number of fractures. Introduction.Efficient numerical simulations of subsurface fluid flows in fractured rocks are of interest for many applications, including water resources management, contaminant transport and dissemination, oil prospecting, and enhanced oil/gas recovery. Among the major difficulties are intrinsic heterogeneity, directionality of the medium, and the multiscale nature of the phenomena, as well as uncertainty in the medium properties. A discrete fracture network (DFN) is a complex threedemensional (3D) structure obtained from intersecting planar fractures. DFN models are frequently preferred to more conventional continuum models as a basis for simulations. A classical approach to the problem is to model fractures as planar ellipses or polygons and stochastically generate DFNs with probabilistic distributions of density, aspect ratio, orientation, size, aperture of fractures, and hydrologic properties [9] and to simulate the flow through the obtained networks. Intensive numerical simulations with several configurations of DFNs and physical parameters are then performed in order to tackle the issue of uncertainty. The flow pattern strongly depends on density and size of fractures and for large-scale simulations different approaches are possible. For dense fracture networks and continuous distribution of size and aspect ratios, flow can be modeled as the flow in an equivalent continuous porous medium where the fracture network pattern leads to the definition of a suitable permeability tensor. For sparse fracture networks with some large fractures that discontinuously increase directionality of the flow, an explicit representation of the fracture network is more reliable. In both cases a stochastic approach to the uncertainty of the parameters is needed and this requires many simulations, so that efficiency and large applicability of numerical algorithms are fundamental issues.

A new approach for solving flow in Discrete Fracture Networks (DFN) is developed in this work by means of the Virtual Element Method. Taking advantage of the features of the VEM, we obtain global conformity of all fracture meshes while preserving a fracture-independent meshing process. This new approach is based on a generalization of globally conforming Finite Elements for polygonal meshes that avoids complications arising from the meshing process. The approach is robust enough to treat many DFNs with a large number of fractures with arbitrary positions and orientations, as shown by the simulations. Higher order Virtual Element spaces are also included in the implementation with the corresponding convergence results and accuracy aspects.

The most challenging issue in performing underground flow simulations in Discrete Fracture Networks (DFN), is to effectively tackle the geometrical difficulties of the problem. In this work we put forward a new application of the Virtual Element Method combined with the Mortar method for domain decomposition: we exploit the flexibility of the VEM in handling polygonal meshes in order to easily construct meshes conforming to the traces on each fracture, and we resort to the mortar approach in order to "weakly" impose continuity of the solution on intersecting fractures. The resulting method replaces the need for matching grids between fractures, so that the meshing process can be performed independently for each fracture. Numerical results show optimal convergence and robustness in handling very complex geometries.

In recent papers [1,2] the authors introduced a new method for simulating subsurface flow in a system of fractures based on a PDE-constrained optimization reformulation, removing all difficulties related to mesh generation and providing an easily parallel approach to the problem. In this paper we further improve the method removing the constraint of having on each fracture a non empty portion of the boundary with Dirichlet boundary conditions. This way, Dirichelet boundary conditions are prescribed only on a possibly small portion of DFN boundary. The proposed generalization of the method in [1, 2] relies on a modified definition of control variables ensuring the non-singularity of the operator on each fracture. A conjugate gradient method is also introduced in order to speed up the minimization process.