From the deal.II FAQ :...quadrilaterals and hexahedra typically provide a significantly better   approximation quality than triangular meshes with the same number of   degrees of freedom; you therefore get more accurate solutions for the   same amount of work...How can this point be made more precise?
How do hexahedral FEM meshes improve approximation quality per degree of freedom, compared to tetrahredal meshes?
pde;finite element;deal.ii
Having written this entry, let me also answer you question :-)The issue in essence boils down to the following: given a function $u(x)$ and its interpolation $u_h(x)$ onto either a triangular or quadrilateral mesh with the same number of unknowns, then which of the two is more accurate? In other words, is$$  \| u - I_h^\text{triangles} u \| < \| u - I_h^\text{quads} u \|,$$or is it the other way around. We know that they have the same asymptotic order ${\cal O}(h^{p+1})$ when using polynomials of order $p$ (in $L_2$). The question is simply which of the two has the smaller interpolation constant. The answer, in practice, is that the inequality above is wrong in general, and that the following is empirically true:$$  \| u - I_h^\text{triangles} u \| \approx 10 \| u - I_h^\text{quads} u \|.$$In other words, the extra polynomial $xy$ on quadrilaterals allows for a 10 times better approximation.(I'm sure you can find counterexamples; this is just an empirical observation for typical solutions of PDEs.)