Numerical models based on seepage-deformation coupling governing equations are often used to simulate soil hydrodynamics and deformation in unsaturated porous media. Among them, Picard iteration method with pressure head as the main variable is widely used because of its simplicity and ability to deal with partial saturation conditions. It is well known that the method is prone to convergence failure under some unfavorable flow conditions and is also computationally time-consuming. In this study, the soil–water characteristic curve (SWCC) of unsaturated soil described by the exponential function is used to linearize the coupling equations to overcome the repeated assembly of nonlinear ordinary differential equations. The finite element method with six-node triangular element is used to discretely linearize the coupling governing equations. Further, the classical Gauss–Seidel iterative method (GS) can be used to solve the linear equations generated from the linearized coupling equations. However, the convergence rate of GS seriously restricts the ill-condition of the linear equations, especially when the condition number of linear equations is much larger than 1.0. Thus, we propose an improved Gauss–Seidel iterative methods MP(m)-GSCMGI by combining multistep preconditioning and cascadic multigrid. The applicability of the proposed methods in simulating variably saturated flow and deformation in unsaturated porous media is verified by numerical examples. The results show that the proposed improved methods have faster convergence rate and computational efficiency than the conventional Picard and GS. The hybrid improved method MP(m)-GSCMGI can achieve more robust convergence and economical simulation.