The thermochemical structure of the lithosphere exerts control on melting mechanisms in the mantle as well as the location of volcanism and ore deposits. Imaging the complex interactions between the lithosphere and asthenospheric mantle requires the joint inversion of multiple data sets and their uncertainties. In particular, the combination of seismic velocity and electrical conductivity with data proxies for bulk composition and elusive minor phases is a crucial step towards fully understanding large-scale lithospheric structure and melting. We apply a novel probabilistic approach for joint inversions of 3D magnetotelluric and seismic data to image the lithosphere beneath southeast Australia. Results show a highly heterogeneous lithospheric structure with deep conductivity anomalies that correlate with the location of Cenozoic volcanism. In regions where the conductivities have been at odds with sub-lithospheric temperatures and seismic velocities, we observe that the joint inversion provides conductivity values consistent with other observations. The results reveal a strong relationship between metasomatized regions in the mantle and i) the limits of geological provinces in the crust, which elucidates the subduction-accretion process in the region; ii) distribution of leucitite and basaltic magmatism; iii) independent geochemical data, and iv) a series of lithospheric steps which constitute areas prone to generating small-scale instabilities in the asthenosphere. This scenario suggests that shear-driven upwelling and edge-driven convection are the dominant melting mechanisms in eastern Australia rather than mantle plume activity, as conventionally conceived. Our study offers an integrated lithospheric model for southeastern Australia and provides insights into the feedback mechanism driving surface processes.