Described here are systems and methods for efficiently simulating wave propagation in one or more viscoelastic media. The systems and methods described here implement an elastodynamic model that includes memory variables for one or more viscoelastic properties. Boundary conditions can also be generated based on a perfectly matched layer (PML) condition. A hybrid implicit/explicit method is used to generate wave simulation data. In this method, equations related to the memory variables and PML are solved using an implicit integrator, and other equations in the elastodynamic model are solved using an explicit method.