Abstract
We develop a novel Hybrid High-Order method for the simulation of Darcy flows in fractured porous media. The discretization hinges on a mixed formulation in the bulk region and on a primal formulation inside the fracture. Salient features of the method include a seamless treatment of nonconforming discretizations of the fracture, as well as the support of arbitrary approximation orders on fairly general meshes. For the version of the method corresponding to a polynomial degree $k\ge 0$, we prove convergence in $h^{k+1}$ of the discretization error measured in an energy-like norm. In the error estimate, we explicitly track the dependence of the constants on the problem data, showing that the method is fully robust with respect to the heterogeneity of the permeability coefficients, and it exhibits only a mild dependence on the square root of the local anisotropy of the bulk permeability. The numerical validation on a comprehensive set of test cases confirms the theoretical results.