The numerical Solution of a discrete Stokes flow problem is a common component of many computational geodynamic models. We compare using the Schur complement reduction (SCR) and a fully coupled (FC) approach to solve geodynamic problems which are characterized by incompressible fluids containing large variations in viscosity. The scalability and robustness of these methods is assessed by examining their sensitivity to the discretization parameter h and the viscosity contrast Delta eta, respectively. We demonstrate that the scaled BFBt preconditioner [Elman, H.C., Howle, V.E., Shadid, J., Shuttleworth, R., Tuminaro, R., 2006. Block preconditioners; based on approximate commutators. SIAM J. Sci. Comput. 27 (5),1651-1668] can be an effective preconditioner for S = G(T)K(-1) G, even when the fluid viscosity varies in space, provided an appropriate scaling is used. The performance of this preconditioner is strongly linked to the scaling employed. Through a number of numerical experiments, we demonstrate that our new scaling strategy ensures that the scaled BFBt preconditioner remains effective in scenarios where the viscosity variations are either locally smooth or discontinuous. Of the solvers considered, applying FGMRES to the pressure Schur complement in combination with BFBt and Our choice of scaling matrices was the most scalable and robust approach. (C) 2008 Elsevier B.V. All rights reserved.