  The construction of finite element approximations in $\mathbf{H}(\mbox{div},
{\Omega})$ usually requires the Piola transformation to map vector polynomials
from a master element to vector fields in the elements of a partition of the
region {\Omega}. It is known that degradation may occur in convergence order if
non affine geometric mappings are used. On this point, we revisit a general
procedure for the improvement of two-dimensional flux approximations discussed
in a recent paper of this journal (Comput. Math. Appl. 74 (2017) 3283-3295).
The starting point is an approximation scheme, which is known to provide
$L^2$-errors with accuracy of order $k+1$ for sufficiently smooth flux
functions, and of order $r+1$ for flux divergence. An example is $RT_{k}$
spaces on quadrilateral meshes, where $r = k$ or $k-1$ if linear or bilinear
geometric isomorphisms are applied. Furthermore, the original space is required
to be expressed by a factorization in terms of edge and internal shape flux
functions. The goal is to define a hierarchy of enriched flux approximations to
reach arbitrary higher orders of divergence accuracy $r+n+1$ as desired, for
any $n \geq 1$. The enriched versions are defined by adding higher degree
internal shape functions of the original family of spaces at level $k+n$, while
keeping the original border fluxes at level $k$. The case $n=1$ has been
discussed in the mentioned publication for two particular examples. General
stronger enrichment $n>1$ shall be analyzed and applied to Darcy's flow
simulations, the global condensed systems to be solved having same dimension
and structure of the original scheme.
