In many seismic applications, a constant quality factor is used to describe the constitutive laws of viscoelastic materials, characterized by frequency-independent attenuation characteristics. In such cases, the frequency dependence of the medium’s properties is not taken into account. To overcome this drawback, we proposed an elegant finite difference time domain implementation, with an auxiliary differential equation technique to explicitly solve any stress-strain relation. This scheme is inherited from the formalism of electromagnetism and is based on the separation of the propagation equations from the constitutive law defined in the stress-strain equation. The conventional assumption of a constant quality factor assumption can then be easily avoided in the modeling of seismic-wave propagation in viscoelastic media. We developed such a method and simulated synthetic traces over a simple, 2D viscoelastic homogeneous medium using a Zener model. Wave propagation phase velocities were estimated by means of dispersion analysis and appeared to match theoretical values over a reasonable frequency range. We also measured the material’s attenuation behavior by studying the quality factor, thereby validating our approach.