The simulation of viscoplasitc flows is still attracting considerable attention in many industrial applications. However, the underlying numerical discretization and regularization may suffer from numerical oscillations, in particular for high Bingham and Reynolds numbers flows. In this work, we investigate the Variational Multiscale stabilized finite element method in solving such flows. We combined it with a posteriori error estimator for anisotropic mesh adaptation, enhancing the use of the Papanastasiou regularization. Computational results are compared to existing data from the literature and new results have demonstrated that the approach can be applied for Bingham numbers higher than 1000 yielding accurate predictions.