In the numerical simulation of groundwater flow, uncertainties often affect the precision of the simulation results. Stochastic and statistical approaches such as the Monte Carlo method, the Neumann expansion method and the Taylor series expansion, are commonly employed to estimate uncertainty in the final output. Based on the first-order interval perturbation method, a combination of the interval and perturbation methods is proposed as a viable alternative and compared to the well-known equal interval continuous sampling method (EICSM). The approach was realized using the GFModel (an unsaturated-saturated groundwater flow simulation model) program. This study exemplifies scenarios of three distinct interval parameters, namely, the hydraulic conductivities of six equal parts of the aquifer, their boundary head conditions, and several hydrogeological parameters (e.g. specific storativity and extraction rate of wells). The results show that the relative errors of deviation of the groundwater head extremums (RDGE) in the late stage of simulation are controlled within approximately ±5% when the changing rate of the hydrogeological parameter is no more than 0.2. From the viewpoint of the groundwater head extremums, the relative errors can be controlled within ±1.5%. The relative errors of the groundwater head variation are within approximately ±5% when the changing rate is no more than 0.2. The proposed method of this study is applicable to unsteady-state confined water flow systems.