**Abstract** : H(div)-conforming finite element approximation spaces are usually formed by locally backtracking vector polynomial spaces defined on the master element by the Piola transformation. The main focus of the present paper is to study the effect of using non-affine elements on the accuracy of three dimensional flux approximations based on such spaces. For instance , instead of order k + 1 for flux and flux divergence obtained with Raviart-Thomas or Nédélec spaces with normal fluxes of degree k, based on affine hexahedra or triangular prisms, reduced orders k for flux and k − 1 for flux divergence may occur for distorted elements. To improve this scenario, a hierarchy of enriched flux approximations is considered, with increasing orders of divergence accuracy, holding for general space stable configurations. The original vector polynomial space is required to be expressed by a decomposition in terms of edge and internal shape functions. The enriched versions are defined by adding internal shape functions of the original family of spaces up to higher degree level k + n, n > 0, while keeping fixed the original border fluxes of degree k. This procedure gives approximations with the same original flux accuracy, but with enhanced divergence order k + n + 1, in the affine case, or k + n − 1 for elements mapped by general multi-linear mappings. The loss of convergence in the flux variable due to quadrilateral face distortions cannot be corrected by including higher order internal functions. Application of these enriched flux spaces to the mixed finite element formulation of a Darcy's model problem is discussed. The computational cost of matrix assembly increases with n, but the global condensed systems to be solved have same dimension and structure of the original scheme.