  A numerical algorithm for solving mantle convection problems with strongly
variable viscosity is presented. Equations for conservation of mass and
momentum for highly viscous and incompressible fluids are solved iteratively by
a multigrid method in combination with pseudo-compressibility and local time
stepping techniques. This algorithm is suitable for large-scale
three-dimensional numerical simulations, because (i) memory storage for any
additional matrix is not required and (ii) vectorization and parallelization
are straightforward. The present algorithm has been incorporated into a mantle
convection simulation program based on the finite-volume discretization in a
three-dimensional rectangular domain. Benchmark comparisons with previous two-
and three-dimensional calculations including the temperature- and/or
depth-dependent viscosity revealed that accurate results are successfully
reproduced even for the cases with viscosity variations of several orders of
magnitude. The robustness of the numerical method against viscosity variation
can be significantly improved by increasing the pre- and post-smoothing
calculations during the multigrid operations, and the convergence can be
achieved for the global viscosity variations up to 1e10.
