We study the numerical solution of thermally convective viscoplastic fluids with yield stress. Following [12], a Bousinessq approximation of the convection effect is considered. The resulting coupled model is then regularized by means of a local regularization technique. After discretization in space, a second order BDF method is used for the time discretization of the regularized problem, leading, in each time iteration, to a nonsmooth system of equations, which is amenable to be solved by generalized Newton methods. A semismooth Newton algorithm with a modified Jacobian is constructed for the solution of the discrete systems. The paper ends with a detailed computational experiment that exhibits the main properties of the numerical approach. (C) 2012 Elsevier B.V. All rights reserved.