In the framework of underground flow simulations in fractured media, fractures may act as preferential paths, and may have a strong impact on the flow. The discrete fracture network (DFN) model allows for an explicit representation of the interconnected fractures. Flux on each fracture is assumed to be driven by Darcy Law, and suitable matching conditions are imposed along fracture intersections, ensuring flux balance and head continuity. The exact displacement of fractures in the network is usually not known, and networks are typically generated sampling geometrical and hydro-geological properties from probabilistic distributions; this stochastic generation is likely to generate geometrical configurations very challenging for the meshing process, which is a major issue in the framework of DFN simulations. Stochasticity in geometrical parameters may also result in a nonsmooth behavior of the quantity of interest with respect to the stochastic parameters. In order to quantify the influence of these stochastic parameters on the output of DFN models, we propose an approach based on the geometric Multi Level Monte Carlo (MLMC) method, applied in conjunction with a well assessed underlying solver for performing DFN flow simulations. Key properties of the solver are its capability of circumventing the need of conforming meshes, and its consequent extreme robustness with respect to geometrical complexities in the network. These features make the solver quite suitable to be used in conjunction with MLMC for the effective application of uncertainty quantification strategies, as they allow to tackle complex geometrical configurations, also with very coarse meshes.