  We propose a semismooth Newton method for non-Newtonian models of
incompressible flow where the constitutive relation between the shear stress
and the symmetric velocity gradient is given implicitly; this class of
constitutive relations captures for instance the models of Bingham and
Herschel-Bulkley. The proposed method avoids the use of variational
inequalities and is based on a particularly simple regularisation for which the
(weak) convergence of the approximate stresses is known to hold. The system is
analysed at the function space level and results in mesh-independent behaviour
of the nonlinear iterations.
