[1] Measurements are often unable to uniquely characterize the subsurface at a desired modeling resolution. In particular, inverse problems involving the characterization of hydraulic properties are typically ill-posed since they generally present more unknowns than data. In a Bayesian context, solutions to such problems consist of a posterior ensemble of models that fit the data (up to a certain precision specified by a likelihood function) and that are a subset of a prior distribution. Two possible approaches for this problem are Markov chain Monte Carlo (McMC) techniques and optimization (calibration) methods. Both frameworks rely on a perturbation mechanism to steer the search for solutions. When the model parameters are spatially dependent variable fields obtained using geostatistical realizations, such as hydraulic conductivity or porosity, it is not trivial to incur perturbations that respect the prior spatial model. To overcome this problem, we propose a general transition kernel (iterative spatial resampling, ISR) that preserves any spatial model produced by conditional simulation. We also present a stochastic stopping criterion for the optimizations inspired from importance sampling. In the studied cases, this yields posterior distributions reasonably close to the ones obtained by a rejection sampler, but with a greatly reduced number of forward model runs. The technique is general in the sense that it can be used with any conditional geostatistical simulation method, whether it generates continuous or discrete variables. Therefore it allows sampling of different priors and conditioning to a variety of data types. Several examples are provided based on either multi-Gaussian or multiple-point statistics.