This paper considers a Leray regularization model of incompressible, nonisothermal fluid flows which uses nonlinear filtering based on indicator functions and introduces an efficient numerical method for solving it. The proposed method uses a multistep, second-order temporal discretization with a finite element (FE) spatial discretization in such a way that the resulting algorithm is linear at each time level and decouples the evolution equations from the velocity filter step. Since the indicator function chosen in this model is mathematically based on approximation theory, the proposed numerical algorithm can be analyzed robustly, i.e., the stability and convergence of the method is provable. A series of numerical tests are carried out to verify the theoretical convergence rates and to compare the algorithm with direct numerical simulation and the usual Leray-model of the flow problem.