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