A viscoacoustic numerical modeling algorithm is developed, tested, and implemented for two-dimensional media. The loss mechanism is simulated according to a standard linear solid model. With anelastic effects defined through superposition of relaxation mechanisms, and incorporating them into the wave equation in terms of memory variables, the time-domain computation of responses of viscoacoustic media with arbitrary spatial variations is feasible. A numerical solution of viscoacoustic responses is obtained by solving either two-dimensional homogeneous or heterogeneous viscoacoustic wave equations. The model is parametrized through stress, strain relaxation times, density, relaxed and unrelaxed modulus, and velocity. On constructing P and S velocity, Q(p), and Q(x) models separately, attenuated and decoupled propagating P and S wave fields are synthesized. According to this scheme one can model widely varying Q in highly heterogeneous media. Simulation of a pure acoustic wave field is regarded as a special case when temporally invariant stress and strain variables are considered. Through systematic analysis, a scheme with a staggered grid is highly desirable for stable computation and reduction of numerical artifacts such as grid dispersion and grid anisotropy. Considerable Nyquist error may arise when a regular grid for a first-derivative operator is used; a staggered grid enables significant improvement. Implementation of a pseudospectral technique with staggered grid provides a more general and more accurate description of wave propagation phenomena. Results from viscoacoustic simulation show that considerable variation of dynamic properties such as a decreased amplitude due to intrinsic Q and temporal shift due to physical wave dispersion phenomena can be analyzed quantitatively. Combined effects of source directivity, complex velocity and attenuation structures may alter significantly in the computed synthetic seismograms.