  We consider a generic dispersive massive gravity theory and numerically study
its resulting modified energy and strain spectra of tensor gravitational waves
(GWs) sourced by (i) fully developed turbulence during the electroweak phase
transition (EWPT) and (ii) forced hydromagnetic turbulence during the QCD phase
transition (QCDPT). The GW spectra are then computed in both spatial and
temporal Fourier domains. We find, from the spatial spectra, that the slope
modifications are weakly dependent on the eddy size at QCDPT, and, from the
temporal spectra, that the modifications are pronounced in the $1$--$10{\rm
nHz}$ range -- the sensitivity range of the North American Nanohertz
Observatory for Gravitational Waves (NANOGrav) -- for a graviton mass $m_{\rm
g}$ in the range $2\times10^{-23}{\rm eV}\lesssim m_{\rm
g}c^2\lesssim7\times10^{-22}{\rm eV}$.
