  We present the Wavelet-based Edge Multiscale Parareal (WEMP) Algorithm,
recently proposed in [Li and Hu, {\it J. Comput. Phys.}, 2021], for efficiently
solving subdiffusion equations with heterogeneous coefficients in long time.
This algorithm combines the benefits of multiscale methods, which can handle
heterogeneity in the spatial domain, and the strength of parareal algorithms
for speeding up time evolution problems when sufficient processors are
available. Our algorithm overcomes the challenge posed by the nonlocality of
the fractional derivative in previous parabolic problem work by constructing an
auxiliary problem on each coarse temporal subdomain to completely uncouple the
temporal variable. We prove the approximation properties of the correction
operator and derive a new summation of exponential to generate a single-step
time stepping scheme, with the number of terms of
$\mathcal{O}(|\log{\tau_f}|^2)$ independent of the final time, where $\tau_f$
is the fine-scale time step size. We establish the convergence rate of our
algorithm in terms of the mesh size in the spatial domain, the level parameter
used in the multiscale method, the coarse-scale time step size, and the
fine-scale time step size. Finally, we present several numerical tests that
demonstrate the effectiveness of our algorithm and validate our theoretical
results.
