  This paper is devoted to the study of Bingham flow with variable density. We
propose a local bi-viscosity regularization of the stress tensor based on a
Huber smoothing step. Next, our computational approach is based on a
second-order, divergence-conforming discretization of the Huber regularized
Bingham constitutive equations, coupled with a discontinuous Galerkin scheme
for the mass density. We take advantage of the properties of the divergence
conforming and discontinuous Galerkin formulations to incorporate upwind
discretizations to stabilize the formulation. The stability of the continuous
problem and the full-discrete scheme are analyzed. Further, a semismooth Newton
method is proposed for solving the obtained fully-discretized system of
equations at each time step. Finally, several numerical examples that
illustrate the main features of the problem and the properties of the numerical
scheme are presented.
