Geophysical methods provide remotely sensed data that are sensitive to subsurface properties and interfaces. Knowledge about discontinuities is important throughout the Earth sciences: for example, the saltwater/freshwater interface in coastal areas drive mixing processes; the temporal development of the discontinuity between frozen and unfrozen ground is indicative of permafrost development; and the regolith-bedrock interface often plays a predominant role in both landslide and critical-zone investigations. Accurate detection of subsurface boundaries and their geometry is challenging when using common inversion routines that rely on smoothness constraints that smear out any naturally occurring interfaces. Moreover, uncertainty quantification of interface geometry based on such inversions is very difficult. In this paper, we present a probabilistic formulation and solution to the geophysical inverse problem of inferring interfaces in the presence of significant subsurface heterogeneity. We implement an empirical-Bayes-within-Gibbs formulation that separates the interface and physical property updates within a Markov chain Monte Carlo scheme. Both the interface and the physical properties of the two sub-domains are constrained to favour smooth spatial transitions and pre-defined property bounds. Our methodology is demonstrated on synthetic and actual surface-based electrical resistivity tomography data sets, with the aim of inferring regolith-bedrock interfaces. Even if we are unable to achieve formal convergence of the Markov chains for all model parameters, we demonstrate that the proposed algorithm offers distinct advantages compared to manual- or algorithm-based interface detection using deterministic geophysical tomograms. Moreover, we obtain more reliable estimates of bedrock resistivity and its spatial variations.