A key challenge in geoscience is that of combining different kinds of geo-information into one geo-model, typically describing the subsurface. This information can be direct information about geological processes, and spatial variability, or it can be indirect information from measurements of properties related to the subsurface, such as geophysical data. Ideally, when such a geo-model has been established, one should be able to quantify information about specific features related to the geo-model, consistent with all information.This integration of geo-information is typically solved using inverse problem theory (Menke, 2012;Tarantola & Valette, 1982a). Fast deterministic methods exist and have been widely used. For such methods, the goal is to obtain one optimal model, such as the simplest possible model, consistent with available information, typically in the form of observed data (Constable et al., 1987;Menke, 2012;Tikhonov, 1963). In practice, in part due to noise on data and model nonlinearities and imperfections, infinitely many models exist that will be consistent with data, and the deterministic approach can in general not account properly for such uncertainty.Probabilistic inversion methods can, in principle, take into account arbitrarily complex information, and integrate the information into one consistent model, as given by the posterior probability distribution. A full analytic expression of the posterior distribution is rarely possible. Instead, sampling methods can be used to generate a