The chemical composition of erupted basalts provides a record of the thermo-chemical state of their source region and the melting conditions that lead to their formation. Here we present the first probabilistic inversion framework capable of inverting both trace and major element data of mafic volcanic rocks to constrain mantle potential temperature, depth of melting, and major and trace element source composition. The inversion strategy is based on the combination of (i) a two-phase multi-component reactive transport model, (ii) a thermodynamic solver for the evolution of major elements and mineral/liquid phases, (iii) a disequilibrium model of trace element partitioning and (iv) an adaptive Markov chain Monte Carlo algorithm. The mechanical and chemical evolution of melt and solid residue are therefore modelled in an internally- and thermodynamically-consistent manner. We illustrate the inversion approach and its sensitivity to relevant model parameters with a series of numerical experiments with increasing level of complexity. We show the benefits and limitations of using major and trace element compositions separately before demonstrating the advantages of a joint inversion. We show that such joint inversion has great sensitivity to mantle temperature, pressure range of melting and composition of the source, even when realistic uncertainties are assigned to both data and predictions. We further test the reliability of the approach on a real dataset from a well-characterised region: the Rio Grande Rift in western North America. We obtain estimates of mantle potential temperature ( ∼ 1340 °C), lithospheric thickness ( ∼ 60 km) and source composition that are in excellent agreement with numerous independent geochemical and geophysical estimates. In particular, this study suggests that the basalts in this region originated from a moderately hot upwelling and include the contribution from a slightly depleted source that experienced a small degree of melt or fluid metasomatism. This component is likely associated with partial melting of the lower portions of the lithosphere. The flexibility of both the melting model and inversion scheme developed here makes the approach widely applicable to assessing the thermo-chemical structure and evolution of the lithosphere-asthenosphere system and paves the way for truly joint geochemical-geophysical inversions.