The numerical solution of a discrete Stokes flow problem is a common component of many computational geodynamic simulations. Such models typically utilize constitutive relationships which may contain large constrasts in material parameters (i.e O(10^8)). In this talk I will describe how to construct an effective preconditioner for the discrete stress tensor gradient and the Schur complement. Numerical experiments will be presented which demonstrate that the proposed preconditioners yield convergence rates which are largely insensitive to the grid size and viscosity contrast.