CSMP seems to have problems creating divergence free velocity fields for quadrilateral FE meshes. Has anyone realised this? We solve, after all

(1)or

(2)and the resulting divergence, as computed with the FV scheme (Stephan's generic and my old one) must be numerically zero (this is the case for triangular meshes, see below). Now, I only get this for a homogeneous permeability field and quadrilateral finite elements, as shown in this image:

If I include a low-permeability region in the centre such that flow is around it, then the divergence becomes immediately non-zero.

Has somebody noticed this and/or a good explanation? We only realised this problem when Christine's Newton iterator got stuck trying to resolve these local source/sink terms that arise from the divergence.