<p>We study a stationary model of doubly diffusive flows with temperature-dependent viscosity on bounded Lipschitz domains in two and three dimensions. A new well-posedness and regularity analysis of weak solutions under minimal assumptions on domain geometry and data regularity are established. A fully non-conforming finite element method based on Crouzeix-Raviart elements, which ensures locally exactly divergence-free velocity fields is explored. Unlike previously proposed schemes, this discretization enables to establish uniqueness of the discrete solutions. We prove the well-posedness of the discrete problem and derive a priori error estimates. An accuracy test is conducted to verify the theoretical error decay rates in flow, Stokes and Darcy regimes.</p>