Land seismic velocity modeling is a difficult task largely related to the description of the near-surface complexities. Full-waveform inversion (FWI) is the method of choice for achieving high-resolution velocity mapping, but its application to land seismic data faces difficulties related to the complex physics, unknown and spatially varying source signatures, and low signal-to-noise ratio (S/N) in the data. Large parameter variations occur in the near surface at various scales causing severe kinematic and dynamic distortions of the recorded wavefield. Some of the parameters can be incorporated in the inversion model, whereas others, due to subresolution dimensions or unmodeled physics, need to be corrected through data preconditioning; a topic not well described for land data FWI applications. We have developed novel algorithms and workflows for surface-consistent data preconditioning using the transmitted portion of the wavefield, S/N enhancement by generation of common midpoint (CMP)-based virtual supershot gathers (VSGs), and robust 1.5D Laplace-Fourier full-waveform inversion (LF-FWI). Our surface-consistent scheme solves residual kinematic corrections and amplitude anomalies via scalar compensation or deconvolution of the near-surface response. S/N enhancement is obtained through the statistical evaluation of volumetric prestack responses at the CMP position, or VSGs. These are inverted via a novel 1.5D acoustic LF-FWI scheme using the Helmholtz wave equation and Hankel domain forward modeling. Inversion is performed with nonlinear conjugate gradients. The method is applied to a complex structure-controlled wadi area exhibiting faults, dissolution, collapse, and subsidence in which the high-resolution FWI velocity modeling helps clarify the geologic interpretation. The developed algorithms and automated workflows provide an effective solution for massive FWI of land seismic data that can be embedded in typical near-surface velocity analysis procedures.