This paper introduces a newly developed wavelet technique for modeling of geophysical flow processes in multiple dimensions. The method utilizes the idea of collocation with the multilevel wavelet approximation. The multilevel structure of the algorithm facilitates the computational adaptation of the grid refinement in regions where sharp variations occur. We have tested this algorithm for viscoelastic flows with viscosity contrasts up to 1012. We have confirmed the findings of stress amplification in the thin highly viscous layer. The new method allows us to conduct the calculations for thin layers and high viscosity contrasts, close enough for realistic mantle-lithospheric interaction. Our results demonstrates the potential usefulness of the wavelet technique in large-scale numerical simulations in the geosciences. ¿ 1997 American Geophysical Union