In this study, we explore an efficient stochastic approach for uncertainty quantification of unconfined groundwater flow in heterogeneous media, where a sparse polynomial chaos expansion (PCE) surrogate model is constructed with the aid of the feature selection method. The feature selection method is introduced to construct a sparse PCE surrogate model with a reduced number of basis functions, which is accomplished by the least absolute shrinkage and selection operator-modified least angle regression and cross-validation. The training samples are enriched sequentially with the quasi-optimal samples until the results are satisfactory. In this study, we test the performance of the sparse PCE method for unconfined flow with the presence of random hydraulic conductivity and recharge, as well as pumping well. Numerical experiments reveal that, even with large spatial variability and high random dimensionality, the sparse PCE approach is able to accurately estimate the flow statistics with greatly reduced computational efforts compared to Monte Carlo simulations.