Investigating the near-surface structure of the Earth through geophysical methods is a crucial aspect in many areas of study in geotechnical engineering, environmental science and exploration. Among several geophysical methods, seismic-based ones are widely used for characterizing near-surface features that are distinguished by contrasts in elastic parameters. In this paper, we model 3D elastodynamic wave propagation and scattering using a method based on domain-type integral representation with Born approximation. We calculate the scattered wavefield by considering the first-order perturbations in density and Lame parameter contrasts of scatterers. Contrasts in Lame parameters can be useful for determining the material properties of subsurface structures in cases of weak contrasts in velocities accompanying considerable Lame parameter variations. We examine the effects of each parameter contrast on a series of models involving a subsurface scatterer. We also compare the seismograms obtained from our method with those from a 3D finite-difference wavefield modelling program, where we observe good agreement between the modelling results. Sensitivity of the wavefield to the perturbation in each model parameter is also examined by calculating and analysing the Frechet derivatives. In general, the method discussed here can provide a solid foundation for prospective imaging studies involving density and Lame parameters simultaneously.