Three different iteration methods for a three-dimensional coordinate-transformed saturated–unsaturated flow model are compared in this study. The Picard and Newton iteration methods are the common approaches for solving Richards’ equation. The Picard method is simple to implement and cost-efficient (on an individual iteration basis). However it converges slower than the Newton method. On the other hand, although the Newton method converges faster, it is more complex to implement and consumes more CPU resources per iteration than the Picard method. The comparison of the two methods in finite-element model (FEM) for saturated–unsaturated flow has been well evaluated in previous studies. However, two iteration methods might exhibit different behavior in the coordinate-transformed finite-difference model (FDM). In addition, the Newton–Krylov method could be a suitable alternative for the coordinate-transformed FDM because it requires the evaluation of a 19-point stencil matrix. The formation of a 19-point stencil is quite a complex and laborious procedure. Instead, the Newton–Krylov method calculates the matrix–vector product, which can be easily approximated by calculating the differences of the original nonlinear function. In this respect, the Newton–Krylov method might be the most appropriate iteration method for coordinate-transformed FDM. However, this method involves the additional cost of taking an approximation at each Krylov iteration in the Newton–Krylov method. In this paper, we evaluated the efficiency and robustness of three iteration methods—the Picard, Newton, and Newton–Krylov methods—for simulating saturated–unsaturated flow through porous media using a three-dimensional coordinate-transformed FDM.