[1] The observed geoid, dynamic topography, and surface plate velocities are controlled by various factors such as density and viscosity variations in the Earth's mantle and strength of the lithosphere. Previous studies have shown that the geoid signal cannot be resolved in details within the framework of a simplified model of the mantle flow considering only radial viscosity variations. Thus, a modeling technique handling both radial and lateral variations of viscosity and other parameters should be used. The spectral method provides a high-accuracy semianalytical solution of the Navier-Stokes and Poisson equations when viscosity is only depth (radially) dependent. In this study, we present the numerical approach, built up on the substantially revised method originally proposed by Zhang and Christensen (1993), for solving the Navier-Stokes equation in the spectral domain with lateral variations of viscosity (LVV). This approach incorporates a number of numerical algorithms providing efficient calculations of the instantaneous Stokes flow in the sphere and taking into account the effects of LVV, self gravitation, and compressibility. In contrast to the traditionally used propagator method, our approach suggests a continuous integration over depth without introducing internal interfaces. Various numerical tests have been employed to test accuracy and efficiency of the proposed technique. Benchmarking of the code shows its ability to solve the mantle convection problems implying strong LVV with high resolution.Components: 16,999 words, 5 figures.