Seismic imaging of onshore targets is very challenging due to the 3D complex near-surface-related effects. In such areas, the seismic wavefield is dominated by elastic and visco-elastic effects such as highly energetic and dispersive surface waves. The interaction of elastic waves with the rough topography and shallow heterogeneities leads to significant converted and scattering energies, implying that both accurate 3D geometry representation and correct physics of the wave propagation are required for a reliable structured imaging. In this manuscript, we present an efficient and flexible full waveform inversion (FWI) strategy for velocity model building in land, specifically in foothill areas.Viscoelastic FWI is a challenging task for current acquisition deployment at the crustal scale. We propose an efficient formulation based on a time-domain spectral element method (SEM) on a flexible Cartesian-based mesh, in which the topography variation is represented by an accurate high-order geometry interpolation. The wave propagation is described by the anisotropic elasticity and isotropic attenuation physics. The numerical implementation of the forward problem includes efficient matrix-vector products for solving second-order elastodynamic equations, even for completely deformed 3D geometries. Complete misfit gradient expressions including attenuation contribution spread into density, elastic parameters and attenuation factors are given in a consistent way. Combined adjoint and forward fields recomputation from final state and previously saved boundary values allows the estimation of gradients with no I/O efforts. Two-levels parallelism is implemented over sources and domain decomposition, which is necessary for 3D realistic configuration. The gradient preconditioning is performed by a so-called Bessel filter using an efficient differential implementation based on the SEM discretization on the forward mesh instead of the costly convolution often-used approach. A non-linear model constraint on the ratio of compressional and shear velocities is introduced into the optimization process at no extra cost.The challenges of the elastic multi-parameter FWI in complex land areas are highlighted through synthetic and real data applications. A 3D synthetic inverse-crime illustration is considered on a subset of the SEAM phase II Foothills model with 4 lines of 20 sources, providing a complete 3D illumination. As the data is dominated by surface waves, it is mainly sensitive to the S-wave velocity. We propose a two-steps data-windowing strategy, focusing on early body waves before considering the entire wavefield, including surface waves. The use of this data hierarchy together with the structurally-based Bessel preconditioning make possible to reconstruct accurately both P- and S-wavespeeds. The designed inversion strategy is combined with a low-to-high frequency hierarchy, successfully applied to the pseudo-2D dip-line survey of the SEAM II Foothill dataset. Under the limited illumination of a 2D acquisition, the model constraint on the ratio of P- and S-wavespeeds plays an important role to mitigate the ill-posedness of the multi-parameter inversion process. By also considering surface waves, we manage to exploit the maximum amount of information in the observed data to get a reliable model parameters estimation, both in the near-surface and in deeper part.The developed FWI frame and workflow are finally applied on a real foothill dataset. The application is challenging due to sparse acquisition design, especially noisy recording and complex underneath structures. Additional prior information such as the logs data is considered to assist the FWI design. The preliminary results, only relying on body waves, are shown to improve the kinematic fit and follow the expected geological interpretation. Model quality control through data-fit analysis and uncertainty studies help to identify artifacts in the inverted models.