The wave equation describing shear wave propagation in three-dimensional (3-D) viscoelastic media is solved numerically with a finite differences time domain (FDTD) method. Solutions are simulated in terms of scatterer velocity waves and verified via comparison to 3-D experimentally acquired wave fields in a heterogenous hydrogel phantom. The numerical algorithm is used as a tool to study wave refraction occurring at the surface of heterogeneities and its effect on complex shear modulus estimation. We used an algebraic reconstruction technique for direct inversion of the wave equation to image the shear modulus and study artifacts produced when reconstructing moduli from 2-D and 3-D velocity data. Although 3-D velocity estimates are required in general, there are object geometries where 2-D reconstructions provide accurate estimations of the material properties.