SUMMARY We present a method to jointly invert surface wave dispersion data and gravity measurements for 3-D shear wave velocity and density models. We implemented a petrophysical approach to combine the kernels of both methodologies in a single process. The synthetic experiments show that jointly inverted models recover shear wave velocity and density better than separate inversions. In particular, density models benefit from the good vertical resolution of surface wave dispersion data, while shear velocity models benefit from the good lateral resolution of gravity data. We also proposed two methods to stabilize the solution when using high-grade polynomials. We applied the methodology to the Los Humeros Geothermal area to demonstrate its applicability in a complex geological scenario. Compared with separate inversion, the joint inversion contributes to enhancing key aspects of the geothermal system by (i) delimitating better the geometry of the caldera deposits in the first 0–2.8 km deep by increasing the vertical resolution in density, (ii) delimitating better the lateral borders of low-Vs bodies at different depths interpreted as a part of a complex magmatic chamber system and (iii) estimating the local shear wave velocity–density relationship that conforms to other known relationships for sedimentary and igneous rocks but with some differences that bring us additional information.