Summary: Multiscale Hybrid Mixed (MHM) method refers to a numerical technique targeted to approximate systems of differential equations with strongly varying solutions. For fluid flow, normal fluxes (multiplier) over macro element boundaries, and coarse piecewise constant potential approximations in each macro element are computed (upscaling). Then, small details are resolved by local problems, using fine representations inside the macro elements, setting the multiplier as Neumann boundary conditions (downscaling). In this work a variant of the method is developed, denoted by MHM- \( \mathbf{H}(\text{div})\), adopting mixed finite elements at the downscaling stage, instead of continuous finite elements used in all previous publications of the method. Thus, this alternative MHM method inherits improvements typical of mixed methods, as better flux accuracy, and local mass conservation at the micro scale level inside the macro elements which are important properties for multi-phase flows in rough heterogeneous media. Different two-scale stable space settings are considered. Vector face functions are supposed to have normal components restricted to a given finite dimensional trace space defined over the macro element boundaries. In each macro element, the internal flux components, with vanishing normal traces, and the potential approximations, may be enriched in different extents: with respect to internal mesh size, internal polynomial degree, or both, the choice being determined by the problem at hands. A unified general error analysis of the MHM-\( \mathbf{H}(\text{div})\) method is presented for all these two-scale space scenarios. Both MHM versions are compared for 2D test problems, with smooth solutions, for convergence rates verification, and for Darcy’s flow in heterogeneous media. MHM-\( \mathbf{H}(\text{div}) 3\) D simulations are presented for a known singular Darcy’s solution, using adaptive macro partitions, and for an oscillatory permeability scenario.