Stokes equations with viscoplastic rheologies play an important role in geodynamics. The nonlinearities of these rheologies make the numerical solution of the resulting systems challenging, and iterative methods often converge slowly or not at all. Yet accurate solutions are critical for representing the physics. We study a basic but representative viscoplastic rheology law which involves a yield stress that is independent of the dynamic pressure, referred to as von Mises yield criterion. We propose a new stress–velocity Newton solution algorithm that treats the stress as an independent variable during the Newton linearization but requires solution only of Stokes systems that are of the usual velocity–pressure form. Comparing the performance of the proposed Newton method with the standard Newton method and the Picard fixed-point method, we observe a significant reduction in the number of iterations and improved stability with respect to problem nonlinearity, mesh refinement, and the polynomial order of the discretization.