  We propose two different discrete formulations for the weak imposition of the
Neumann boundary conditions of the Darcy flow. The Raviart-Thomas mixed finite
element on both triangular and quadrilateral meshes is considered for both
methods. One is a consistent discretization depending on a weighting parameter
scaling as $\mathcal O(h^{-1})$, while the other is a penalty-type formulation
obtained as the discretization of a perturbation of the original problem and
relies on a parameter scaling as $\mathcal O(h^{-k-1})$, $k$ being the order of
the Raviart-Thomas space. We rigorously prove that both methods are stable and
result in optimal convergent numerical schemes with respect to appropriate
mesh-dependent norms, although the chosen norms do not scale as the usual
$L^2$-norm. However, we are still able to recover the optimal a priori
$L^2$-error estimates for the velocity field, respectively, for high-order and
the lowest-order Raviart-Thomas discretizations, for the first and second
numerical schemes. Finally, some numerical examples validating the theory are
exhibited.
