Numerical simulation of saturated-unsaturated subsurface flows is widely used in many branches of science and engineering, and rapid developments in computer technology have enabled not only one-dimensional but also multi-dimensional simulations on a personal computer. However, a multidimensional subsurface flow simulation still incurs heavy load on computational resources, particularly for simulations in wide regions with long periods. An Iterative Alternating Direction Implicit (IADI) scheme has certain advantages in terms of computational cost and algorithmic simplicity. However, it is barely used at present because it occasionally incurs numerical instabilities and convergence difficulties. Another reason is that three-dimensional simulations cannot be performed by the original IADI scheme. This study has proposed an advanced IADI algorithm for solving the saturated-unsaturated flow equation: this advanced scheme is more numerically stable than the original IADI scheme and can be used for three-dimensional simulations. The performance of the proposed scheme was assessed through test simulations. In all the simulations, the new method was shown to be faster than the fully implicit scheme linearized by the modified Picard iteration method while still yielding very similar results. (C) 2011 Elsevier By. All rights reserved.