  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 \v{e} 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.
