The characterization of hydraulic properties of aquifers and soils is essential to better predict water flow in the subsurface and the transport of heat or solutes. Typically, not enough direct data (e.g., hydraulic conductivity) are available to characterize the heterogeneous subsurface. Thus, additional indirect data (e.g., hydraulic heads) are important for improving characterization and, in turn, predictions by subsurface flow and transport models. In the geostatistical context, the resulting inverse problem for subsurface problems is typically underdetermined. Therefore, approaches were developed which limited the number of independent parameters to be estimated, either by defining a limited number of zones with constant parameters (Carrera & Neuman, 1986) or parameterizing the spatially variable parameter field by a geostatistical function with a few unknown parameters (Kitanidis & Vomvoris, 1983). Later, methods were formulated to estimate a series of equally likely solutions to the groundwater inverse problem, either by an ensemble-based variational data assimilation approach (Gómez-Hernández et al., 1997) or by a sequential data assimilation approach for an ensemble of random parameter fields (Chen & Zhang, 2006). In summary, the combination of regularization and casting the problem in a stochastic framework, helped to tackle groundwater inverse problems.Bayesian inversion has been widely used for model parameter inference (