We extend the balancing domain decomposition by constraints (BDDC) method to flows in porous media discretised by mixed‐hybrid finite elements with combined mesh dimensions. Such discretisations appear when major geological fractures are modelled by one‐dimensional or two‐dimensional elements inside three‐dimensional domains. In this set‐up, the global problem and the substructure problems have a symmetric saddle‐point structure, containing a ‘penalty’ block due to the combination of meshes. We show that the problem can be reduced by means of iterative substructuring to an interface problem, which is symmetric and positive definite. The interface problem can thus be solved by conjugate gradients with the BDDC method as a preconditioner. A parallel implementation of this algorithm is incorporated into an existing software package for subsurface flow simulations. We study the performance of the iterative solver on several academic and real‐world problems. Numerical experiments illustrate its efficiency and scalability. Copyright © 2015 John Wiley & Sons, Ltd.