A mixed least-squares finite element model with spectral/hp approximations was developed to analyze a three-dimensional natural convection of non-Newtonian fluid, which obeys the Carreau-Yasuda constitutive model. The finite element model consists of velocity, pressure, stress, temperature, and heat flux as the field variables. The least-squares formulation provides a variational framework for the Navier-Stokes equations and no compatibility of the approximation spaces for field variables is imposed. Also, the use of spectral/hp elements in conjunction with the least-squares formulation alleviates various forms of locking which often appear in low-order least-squares finite element models for incompressible viscous flows and yields accurate results with exponential convergence. To verify and validate the code for Newtonian fluids, the current results are compared with the numerical and experimental results in the literature. Then, the parametric effects of the Carreau-Yasuda model on the flow characteristics are studied.