Searching for the precise solution of free surface has remained the main bottleneck in analyzing the unconfined seepage problem for earth-rock dams. To seek a tradeoff between searching efficiency and accuracy, and as a first attempt to explore the full impacts of seepage anomalies and its varying location on the multiple seepage parameters for the typical heterogeneity issues, this paper proposes a novel approach by combining the smoothed finite element method (S-FEM) and gradient smoothing technique (GST). Considering the practicality limitation due to huge calculation burden resulted from solving the internal integral of elements, the proposed method optimizes the shape functions calculation and explores the specific influence of division manner of smoothing cells on the computational efficiency and accuracy. To demonstrate the novelty and computational efficiency of the GST-based methodology, we performed a series classical model tests involving classic homogeneous and complicated inhomogeneous media models. We focused on the specific impacts of anomalies with varying degrees of permeability on the multiple parameters for seepage inversion, including water head, free surface solution, overflow point, seepage velocity and fluid pressure. The evidence implies that the significant influence of anomalies existence deserves high concern. More importantly, the results also demonstrate that our method improves the capacity of maintaining the high accuracy of seepage solutions with a preferable efficiency in handling the heterogeneous cases of earth dam, which provides great potential for developing multi-parameter seepage inversion, the accuracy of our GST-based seepage simulation is a significant step towards commercially-viable applications of multi-parameter seepage inversion.