We design a numerical algorithm for wave simulation in anelastic media in the presence of free surface, which can be used to model seismic waves at the Earth's surface and ultrasonic waves in heterogeneous materials. The stress-strain relation is based on the Kelvin-Voigt mechanical model, which has the advantage of not requiring additional field variables. The model requires two anelastic parameters and twice the spatial derivatives of the lossless case. The high-frequency components of the wave field are more attenuated than the low-frequency components, with the attenuation factors being approximately proportional to the square of the frequency. The modeling simulates 3-D waves by using the Fourier and Chebyshev methods to compute the spatial derivatives along the horizontal and vertical directions, respectively. We stretch the mesh in the vertical direction to increase the minimum grid spacing and reduce the computational cost. Instabilities of the Chebyshev differential operator due to the implementation of the free-surface boundary conditions are solved with a characteristic approach, where the characteristic variables are evaluated at the source central frequency. The results of the modeling are verified by comparisons to the analytical solutions for Lamb's problem and propagation in unbounded homogeneous media. Examples illustrating the propagation of Rayleigh and Love waves are presented.