In the present contribution we compare different mixed least-squares finite element formulations (LSFEMs) with respect to computational costs and accuracy. In detail, we consider an approach for Newtonian fluid flow, which is described by the incompressible Navier-Stokes equations. Starting from the residual forms of the equilibrium equation and the continuity condition, various first-order systems are derived. From these systems least-squares functionals are constructed by means of L 2-norms, which are the basis for the associated minimization problems. The first formulation under consideration is a div-grad first-order system resulting in a threefield formulation with stresses, velocities, and pressure as unknowns. This S-V-P formulation is approximated in H(div) × H 1 × L 2 on triangles and for comparison also in H 1 × H 1 × L 2 on quadrilaterals. The second formulation is the well-known div-curl-grad first-order velocityvorticity-pressure (V-V-P) formulation. Here all unknowns are approximated in H 1 on quadrilaterals. Besides some numerical advantages, as e.g. an inherent symmetric structure of the system of equations and a directly available error estimator, it is known that least-squares methods have also a drawback concerning mass conservation, especially when lower-order elements are used. Therefore, the main focus of the work is on performance and accuracy aspects on the one side for finite elements with different interpolation orders and on the other side on the usage of efficient solvers, for instance of Krylov-space or multigrid type. In order to demonstrate the capability of the formulations the results for some well-known benchmark problems are presented and discussed.