In contrast to seismic velocities, density is usually not well retrieved by seismic tomography. While gravity data are sensitive to density variations, inverting gravity data alone leads to poorly constrained solutions. Therefore, the joint inversion of gravity and seismic data is a promising way to reveal Earth’s 3D structure of density and seismic wave velocities. Seismic data considered in the joint inversion methods were so far often limited to simple observables (e.g. arrival times) and very few studies focused on full waveforms. In this work, we developped a workflow for the collaborative inversion of teleseismic waveforms and gravity anomalies. Through a series of synthetic tests using models derived from tomographic models of the western Pyrenees, we found that an efficient strategy for the joint inversion of seismic and gravity data is to first perform seismic inversion and then gravity inversion. The 3D full waveform inversion of teleseismic P waves is based upon the open source Python package TeleSeisFlows.