Purpose <p>The current study examines propagation of love-type waves with an initially stressed heterogeneous viscoelastic orthotropic medium is bonded between a heterogeneous dry sandy medium and an initially stressed heterogeneous transversely isotropic half-space. The half-space and layer interfaces are taken as sliding contacts with two distinct sliding parameters. The orthotropic medium is modelled using a Kelvin–Voigt type viscoelastic formulation. The vertical non-homogeneity is modelled as a quadratic function in the upper medium, and as linear functions in both the intermediate layer and the half-space, varying with depth. </p> Methods <p>Using the approach of variable separation and appropriate replacements, the governing partial differential equations of motion are solved analytically. At both interfaces, sliding-type boundary conditions are applied, permitting partial slip rather than perfect bonding. With MATHEMATICA software, dispersion and displacement were analysed numerically and graphically. A comparison with the classical Love wave equation and a number of specific situations have been used to test the validity of our model. </p> Results <p>The results show clear layer-dependent behavior: heterogeneity and initial stress either accelerate or delay Love waves depending on the layer in which they occur, whereas interfacial slip factors have opposite effects at different interfaces. In particular, heterogeneity in the viscoelastic layer increases phase velocity while decreasing damping, but heterogeneity in the sandy layer and half-space improves both. </p> Conclusion <p>These findings show that Love-wave features are influenced not only by material qualities, but also by the location of heterogeneity, stress, and slip, offering new physical insights for seismic site response research, fault zone studies, and near-surface wave-based characterization.</p>