Archean tectonic evolution and metallogenic endowment in cratonic areas are intrinsically linked to lithological diversity and complex architectures. Unfortunately, reliable models of the deep structure of such areas are scarce due to poor surface exposure and limited availability of drill holes and high-resolution geophysical information. By combining multiple geophysical data sets, the strategy of cross-gradient joint geophysical inversion permits improved subsurface imaging of geological heterogeneity and structure. In this work we showcase the application of cross-gradient joint inversion of gravity and magnetic data to characterize the heterogeneity of the Leonora district in the eastern Yilgarn Craton, Western Australia. Our results harmonize both geophysical data types into integrated images of the top 15 km of the crust and successfully predict observations made by a seismic experiment and surface geological mapping.