Electrical resistivity tomography is an ill‐posed and nonlinear inverse problem commonly solved through deterministic gradient‐based methods. These methods guarantee a fast convergence towards the final solution, but the local linearization of the inverse operator impedes accurate uncertainty assessments. On the contrary, numerical Markov chain Monte Carlo algorithms allow for accurate uncertainty appraisals, but appropriate Markov chain Monte Carlo recipes are needed to reduce the computational effort and make these approaches suitable to be applied to field data. A key aspect of any probabilistic inversion is the definition of an appropriate prior distribution of the model parameters that can also incorporate spatial constraints to mitigate the ill conditioning of the inverse problem. Usually, Gaussian priors oversimplify the actual distribution of the model parameters that often exhibit multimodality due to the presence of multiple litho‐fluid facies. In this work, we develop a novel probabilistic Markov chain Monte Carlo approach for inversion of electrical resistivity tomography data. This approach jointly estimates resistivity values, litho‐fluid facies, along with the associated uncertainties from the measured apparent resistivity pseudosection. In our approach, the unknown parameters include the facies model as well as the continuous resistivity values. At each spatial location, the distribution of the resistivity value is assumed to be multimodal and non‐parametric with as many modes as the number of facies. An advanced Markov chain Monte Carlo algorithm (the differential evolution Markov chain) is used to efficiently sample the posterior density in a high‐dimensional parameter space. A Gaussian variogram model and a first‐order Markov chain respectively account for the lateral continuity of the continuous and discrete model properties, thereby reducing the null‐space of solutions. The direct sequential simulation geostatistical method allows the generation of sampled models that honour both the assumed marginal prior and spatial constraints. During the Markov chain Monte Carlo walk, we iteratively sample the facies, by moving from one mode to another, and the resistivity values, by sampling within the same mode. The proposed method is first validated by inverting the data calculated from synthetic models. Then, it is applied to field data and benchmarked against a standard local inversion algorithm. Our experiments prove that the proposed Markov chain Monte Carlo inversion retrieves reliable estimations and accurate uncertainty quantifications with a reasonable computational effort. [ABSTRACT FROM AUTHOR]