  This paper considers a Leray regularization model of incompressible,
non-isothermal fluid flows which uses nonlinear filtering based on indicator
functions, and introduces an efficient numerical method for solving it. The
proposed method uses a multi-step, 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-$\alpha$ model of the flow problem.
