Unsaturated soil behaviors characterize the failure mechanisms of geotechnical infrastructures with transient seepage conditions. Therefore, an accurate estimate of the unsaturated groundwater flow is vital in improving hazard management and assessment. This study attempts to develop a numerical scheme for 2-D transient analysis under unsaturated conditions. First, the unsaturated groundwater flow was described using the mass conservation law. Then, the Finite Difference Method and Backward Euler approximation were applied for space and time discretization, respectively. Furthermore, the simple Picard iteration was applied to linearize the governing equation. The reliability of the presented method was verified with the analytical solution. The evaluation results demonstrated the sufficiency of the proposed method, quantitatively expressed by the maximum error of 0.04% for opened boundary conditions and 0.15% for closed boundary conditions. The significant advantage of the proposed method is the flexibility with various soil-water characteristic curve models and associated hydraulic conductivity functions, which helps to improve the applicability in practice.