Developing a numerical and algorithmic tool which correctly identies unyielded regions in yield stress uid ow is a challenging task. Two approaches are commonly used to handle the singular behaviour at the yield surface, i.e. the Augmented Lagrangian approach and the regularization approach, respectively. Generally in the regularization approach, solvers do not perform eciently when the regularization parameter gets very small. In this work, we use a formulation introducing a new auxiliary stress. The three eld formulation of the yield stress uid corresponds to a regularization-free Bingham formulation. The resulting set of equations arising from the three eld formulation is solved eciently and accurately by a monolithic nite element method. The velocity and pressure are discretized by the higher order stable FEM pair Q_2/P_1^disc and the auxiliary stress is discretized by the Q_2 element. Furthermore, this problem is highly nonlinear and presents a big challenge to any nonlinear solver. Therefore, we developed a new adaptive discrete Newton method, which evaluates the Jacobian with the divided dierence approach. We relate the step length to the rate of the actual nonlinear reduction for achieving a robust adaptive Newton method. We analyse the solvability of the problem along with the adaptive Newton method for Bingham uids by doing numerical studies for a prototypical conguration "viscoplastic uid ow in a channel".