The sustainable management of groundwater demands a faithful characterization of the subsurface. This, in turn, requires information which is generally not readily available. To bridge the gap between data need and availability, numerical models are often used to synthesize plausible scenarios not only from direct information but also additional, indirect data. Unfortunately, the resulting system characterizations will rarely be unique. This poses a challenge for practical parameter inference: Computational limitations often force modelers to resort to methods based on questionable assumptions of Gaussianity, which do not reproduce important facets of ambiguity such as Pareto fronts or multi-modality. In search of a remedy, an alternative could be found in Stein Variational Gradient Descent, a recent development in the field of statistics. This ensemble-based method iteratively transforms a set of arbitrary particles into samples of a potentially non-Gaussian posterior, provided the latter is sufficiently smooth. A prerequisite for this method is knowledge of the Jacobian, which is usually exceptionally expensive to evaluate. To address this issue, we propose an ensemble-based, localized approximation of the Jacobian. We demonstrate the performance of the resulting algorithm in two cases: a simple, bimodal synthetic scenario, and a complex numerical model based on a real-world, pre-alpine catchment. Promising results in both cases - even when the ensemble size is smaller than the number of parameters - suggest that Stein Variational Gradient Descent can be a valuable addition to hydrogeological parameter inference.