The information provided by different geophysical data sets (gravimetric, magnetic, seismic, etc.) can be used, together with petrophysical and geostatistical information, to estimate the major lithological properties of the rocks within the studied volume. The formalization of this inverse problem requires a joint representation and parameterization of the different media properties in the model. The information relating rock properties together couples the inversion of the plural geophysical data sets and allows one to relate the observations with the lithological parameters of the model. The representation by probability density functions (pdfs) of the different types of information entering the problem is also required and provides the mathematical framework to formulate their combination. The resulting joint posterior pdf is composed of two factors: the joint likelihood function, which is the product of independent likelihood functions associated with each geophysical data set, and the joint prior pdf. The latter decomposes, following a partition of the model parameter space in a primary (lithological) subspace and secondary (physical) subspace, as a marginal pdf over the primary model parameter space and a conditional pdf over the secondary model parameter space. A Markov‐chain Monte Carlo method was adapted to sample joint models from the posterior pdf: (1) the method starts with a Markov‐chain sampling primary models from the marginal prior pdf, (2) the chain is extended to the joint model space, by sampling from the conditional pdf of the secondary parameters with respect to the primary parameters, and (3) it is modified to sample from the posterior pdf, by applying the Metropolis rule, which uses the evaluation of the joint likelihood function to accept or reject model transitions in the sampling chain. Finally, posterior marginal or posterior conditional pdfs for the model parameters or the model properties can be straightforwardly calculated from the set of joint models sampled by the chain.