Abstract In this work we present and analyze a finite element scheme yielding discontinuous Galerkin approximations to the solutions of the stationary Boussinesq system for the simulation of non-isothermal flow phenomena. The model consists of a Navier–Stokes-type system, describing the velocity and the pressure of the fluid, coupled to an advection-diffusion equation for the temperature. The proposed numerical scheme is based on the standard interior penalty technique and an upwind approach for the nonlinear convective terms and employs the divergence-conforming Brezzi–Douglas–Marini (BDM) elements of order k for the velocity, discontinuous elements of order k - 1 {k-1} for the pressure and discontinuous elements of order k for the temperature. Existence and uniqueness results are shown and stated rigorously for both the continuous problem and the discrete scheme, and optimal a priori error estimates are also derived. Numerical examples back up the theoretical expected convergence rates as well as the performance of the proposed technique.