## 1 Introduction  
Thermally driven flows and buoyancy have important practical applications, such as ventilation problems, indoor comfort, thermally driven ocean currents, dynamics of plumes, and so on [1, 2].  
The dimensionless number that measures the importance of buoyancy forces due to thermal gradients compared to viscous forces is the Grashof number. Many times, the Grashof number is replaced by the Rayleigh number, which equals the Grashof number divided by the Prandtl number and represents the ratio of buoyancy forces to the viscous force and thermal diffusion [3].  
From a computational mechanics point of view, flows at large Rayleigh numbers are viewed as hard. As the Rayleigh number increases, the flow evolves from steady laminar to unsteady laminar conditions. A further increase of the Rayleigh number makes the flow turbulent. The efficient computation of high Rayleigh numbers requires being capable of taking large time-steps and that the CPU time per time-step be small.  
The literature on computational strategies to tackle thermally coupled low-speed flows is extremely vast. Almost every method and algorithm has been tested for this type of flow, both in segregated and the coupled formats. Regarding examples related to finite elements using the [PERSON] method, one can mention [4; 5]; using stabilized methods [6; 7; 8]; and using the Finite Increment Calculus (FIC) strategy [9]. [PERSON] and [PERSON] present a second-order-characteristics finite-element method for thermally coupled incompressible flow problems in [10]. Examples of finite-element methods used to study the turbulent regime through the inclusion of turbulence models can be found in [2], whereas [11; 12] employ quasi-static sub-scales to embrace turbulence, [13] represents the fine scales through dynamical bubbles and [14; 15] use transient sub-scales to model fine-scales. Weak imposition of boundary conditions is studied in [16].  
Other models used to simulate thermal convection are the anelastic model and the acoustically filtered Navier-Stokes equations (see [4; 17] and references therein). In the anelastic model, the temporal derivative of density in the continuity equation is ignored, and, thermodynamically, the density is considered to be only a function of temperature, so it does not depend on the dynamical pressure. In the filtered Navier-Stokes model, by taking the zero Mach number limit of the compressible equations, a set of elliptic equations is obtained where the acoustic waves are removed from the model; thermodynamically, the density does not depend on the dynamical part of the pressure.  
There are various unified approaches using stabilized finite-element methods to tackle compressible and incompressible flows, such as [18], but the work of this paper is based on [19; 20]. Ref. [21] comments on advantages of this formulation in terms of treating temperature dependent flows with a unified approach based on pressure or entropy variables.  
In this paper, two strategies to compute buoyant flows are compared: one based on the classical Boussinesq approximation and the other one based on considering the thermal compressibility of the liquid or gas. Furthermore, three stabilized methods based on a unified approach are applied to buoyancy flows and are compared and tested. Regarding the algorithms, a wide range of updated stabilization parameters are compared for this application, and conclusions are reached about the best choices.