The teleseismic P receiver functions are customarily inverted to attain the seismic velocities beneath a seismic station. Surface wave dispersion data are often added to reduce the effect of the non-uniqueness. The combination of P receiver function and surface wave works well in resolving the structures in the crust and uppermost mantle, but is less effective in characterizing greater (lithosphere and asthenosphere) depths due to the interference from crustal multiples. A solution to this problem is jointly to model teleseismic S receiver functions with surface wave and P receiver functions. This study adopts a fast, one-dimensional (1-D) inversion scheme. To avoid the effect of multidimensional structures away from the seismic station, we eliminate multiples that reverberate between the surface and interfaces below a restriction depth (RD), as well as S-to-P conversions below an inversion depth (ID). P-to-S conversions off the interfaces above the half-space and S-to-P conversions above the ID and multiples above the RD are properly modelled. This approach favours ray paths travelling close to stations and is, therefore, more suitable for 1-D inversions. We perform numerical experiments with and without noise and highlight the advantages of a joint receiver function and surface wave analysis. © 2010 Springer Science+Business Media B.V.