We have developed a wavelet-multiscale adjoint scheme for the elastic full-waveform inversion of seismic data, including body waves (BWs) and surface waves (SWs). We start the inversion on the SW portion of the seismograms. To avoid cycle skipping and reduce the dependence on the initial model of these dispersive waves, we commence by minimizing an envelope-based misfit function. Subsequently, we proceed to the minimization of a waveform-difference (WD) metric applied to the SWs only. After that, we fit BWs and SWs indiscriminately using a WD misfit metric. In each of these three steps, we guide the iterative inversion through a sequence of nested subspace projections in a wavelet basis. SW analysis preserves a wealth of near-surface features that would be lost in conventional BW tomography. We used a toy model to illustrate the dispersive and cycle-skipping behavior of the SWs, and to introduce the two ways by which we combat the nonlinearity of waveform inversions involving SWs. The first is the wavelet-based multiscale character of the method, and the second the envelope-based misfit function. Next, we used an industry synthetic model to perform realistic numerical experiments to further develop a strategy for SW and joint SW as well as BW tomography. The effect of incorrect density information on wave-speed inversions was also evaluated. We ultimately formalize a flexible scheme for full-waveform inversion based on adjoint methods that includes BWs and SWs, and also considers P- and S-wave speeds, as well as density. Our method is applicable to waveform inversion in exploration geophysics, geotechnical engineering, regional, and global seismology.