Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP)We develop two families of mixed finite elements on quadrilateral meshes for approximating (u,p) solving a second order elliptic equation in mixed form. Standard Raviart-Thomas (RT) and Brezzi-Douglas-Marini (BDM) elements are defined on rectangles and extended to quadrilaterals using the Piola transform, which are well-known to lose optimal approximation of del.u. Arnold-Boffi-Falk spaces rectify the problem by increasing the dimension of RT, so that approximation is maintained after Piola mapping. Our two families of finite elements are uniformly inf-sup stable, achieve optimal rates of convergence, and have minimal dimension. The elements for u are constructed from vector polynomials defined directly on the quadrilaterals, rather than being transformed from a reference rectangle by the Piola mapping, and then supplemented by two (one for the lowest order) basis functions that are Piola mapped. One family has full H(div)-approximation (u, p, and del.u are approximated to the same order like RT) and the other has reduced H(div)-approximation (p and del.u are approximated to one less power like BDM). The two families are identical except for inclusion of a minimal set of vector and scalar polynomials needed for higher order approximation of del.u and p, and thereby we clarify and unify the treatment of finite element approximation between these two classes. The key result is a Helmholtz-like decomposition of vector polynomials, which explains precisely how a divergence is approximated locally. We develop an implementable local basis and present numerical results confirming the theory.54633323356U.S. National Science Foundation [DMS-1418752]Sao Paulo Research Foundation (FAPESP) [2013/19021-8]Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP)