This paper presents a time-domain boundary element method (BEM) for viscoelastic wave scattering by many cavities in 3D infinite space. The convolution quadrature method (CQM) is applied to the convolutions of the time-domain boundary integral equation in order to improve the numerical stability of BEM. Scattering of an incident wave by cavities in a viscoelastic solid is solved by the developed time-domain BEM. The incident wave propagating in a viscoelastic solid in time-domain is obtained using the inverse Fourier transform of that in frequency-domain. In addition, the hybrid parallelization using OpenMP and MPI is used to save the computational time and memory.