This paper is devoted to the mathematical modelling and numerical simulation of basic mechanisms that drive multicontinuum systems in two- and three-dimensional porous media. We state a general mathematical model within the framework of mixture theory, capable to describe and predict flow patterns occurring in permeable media, a problem exhibiting high degrees of heterogeneity and large disparities in physical scales. The model equations result in a strongly coupled and nonlinear system of partial differential equations (PDE) that are written in terms of phase and barycentric mixture velocities, phase pressure, and saturation. We construct accurate, robust and reliable hybrid (or combined) mixed finite element -- primal finite volume-element methods for the discretization of the underlying equations on unstructured meshes. These schemes rely on mixed Brezzi-Douglas-Marini approximations of phase and total velocities, on piecewise constant elements for the approximation of phase or total pressures, and on a primal formulation using discontinuous finite volume elements defined on a dual diamond mesh to approximate scalar mixture constituents of interest (such as volume fraction, total density, saturation, etc.). A second order backward difference formula is employed in the approximation of time derivatives. Several numerical test cases are presented and discussed in detail, illustrating the validity of the approach proposed herein along with the accuracy of the mixed--primal method. We comment as well on the applicability of the model and numerical scheme to the study of other related systems of wide interest.