The time-dependent Navier-Stokes equations and the energy balance equation for an incompressible, constant property fluid in the Boussinesq approximation are solved by a least-squares finite element method based on a velocity-pressure-vorticity-temperature-heat-flux (u-P-omega-T-q) formulation discretized by backward finite differencing in time. The discretization scheme leads to the minimization of the residual in the l2-norm for each time step. Isoparametric bilinear quadrilateral elements and reduced integration are employed. Three examples, thermally driven cavity flow at Rayleigh numbers up to 10(6), lid-driven cavity flow at Reynolds numbers up to 10(4) and flow over a square obstacle at Reynolds number 200, are presented to validate the method.