The past ten years have been marked by dramatic advances in four seemingly isolated research fields: Thermodynamic modelling of minerals and rocks at high PT conditions, numerical simulation of the thermomechanical behaviour of the Earth’s interior, efficient decomposition techniques to solve complex simulation-based problems, and probabilistic data analysis and inversion methods. All these disciplines/techniques have individually created true “revolutions” in the way we understand and model natural systems, including the interior of the Earth. Nevertheless, a more profound understanding is still ahead of us from the formal combination of these disciplines/techniques into a single operational framework to study the physical state of the Earth’s interior. In this contribution, I will present and discuss the concept of multi-observable probabilistic tomography or “thermochemical tomography”. This new kind of tomography is particularly designed for studies of the fundamental thermodynamic variables of the Earth’s lithosphere, namely temperature, pressure and chemical composition. Once these variables are known, all physical parameters of interest (e.g. seismic velocities, density, viscosity, conductivity, etc), as well as traditional tomography images, are also retrieved in a thermodynamically-consistent way. The method is built on a simulation-based inversion technique where multiple satellite (e.g. gravity gradients, geoid height, etc) and land-based (e.g. seismic, plate motions, heat flow, etc) datasets can be jointly inverted to maximize the physical consistency of the resulting Earth model. Assembling this large problem required a collaborative effort between thermodynamicists, mineral physicists, geophysicists and geochemists, and marks the first step towards a full coupling between geophysics, geodynamics, thermodynamics, and geochemistry. I will present results for both synthetic and real case studies, which serve to highlight the advantages and limitations of this approach.