Abstract: The existing semi-analytical numerical method can not be used to simulate unconfined flow of groundwater. Therefore, a semi-analytical numerical approach for the problems of unsteady groundwater flow in unconfined aquifers is proposed by using the Nueman's model for the response of the water table, and the decoupled weighted residual equations are derived by means of the Galerkin's method and the orthogonality of the trigonometric series. A computer program is developed for the semi-analytical numerical analysis of three-dimensional groundwater flow in unconfined aquifers, and the validities of the present method and program are verified by comparisons with the existing analytical solutions. The three-dimensional pattern of unconfined flow due to a fully penetrating well and the delayed response of the water table are studied based on the numerical results. Moreover, the proposed method can be reduced to the existing method for confined flow if the specific yield of unconfined aquifers is ignored. Finally, the applicability of the semi-analytical numerical method for simulating the three-dimensional flow in the aquifer systems consisting of an unconfined aquifer, an aquitard and a confined aquifer is demonstrated through an additional numerical application. The proposed method can provide a high efficient approach for the three-dimensional simulation of groundwater flow and land subsidence due to groundwater withdrawal in layered aquifer systems.