The propagation of anelastic surface waves presents substantial differences compared to the elastic case. Therefore, a forward modeling scheme to study surface waves in an anelastic earth is of some relevance to the geophysical problem.When propagating anelastic waves, accuracy is very important; in particular, numerical dispersion should not be confused with physical velocity dispersion. One modeling algorithm, based on the velocity-stress elastodynamic equations, uses a spectral method with a Chebychev expansion in the vertical direction. This approach allows the calculation of the spatial derivatives with high accuracy, and an effective incorpora-tion of the free surface boundary conditions since the method is not periodic. However, a direct application of the traction-free boundary conditions without regard to the other variables produces numerical instabilities. To solve this problem, a boundary treatment based on characteristics is implemented that results in a wave equation for the surface that automatically includes the boundary conditions.Two examples of wave propagation in an unconsolidated weathering zone are presented. The first model is homogeneous and the second structure contains a vertical interface separating an anelastic medium from an elastic region. The results indicate that the modeling correctly describes the anelastic properties of Rayleigh waves, even in the presence of a strong contrast in the material parameters.