  We introduce a method which provides accurate numerical solutions to
fractional-in-time partial differential equations posed on $[0,T] \times
\Omega$ with $\Omega \subset \mathbb{R}^d$ without the excessive memory
requirements associated with the nonlocal fractional derivative operator. Our
approach combines recent advances in the development and utilization of
multivariate sparse spectral methods as well as fast methods for the
computation of Gauss quadrature nodes with recursive non-classical methods for
the Caputo fractional derivative of general fractional order $\alpha > 0$. An
attractive feature of the method is that it has minimal theoretical overhead
when using it on any domain $\Omega$ on which an orthogonal polynomial basis is
already available. We discuss the memory requirements of the method, present
several numerical experiments demonstrating the method's performance in solving
time-fractional PDEs on intervals, triangles and disks and derive error bounds
which suggest sensible convergence strategies. As an important model problem
for this approach we consider a type of wave equation with time-fractional
dampening related to acoustic waves in viscoelastic media with applications in
the physics of medical ultrasound and outline future research steps required to
use such methods for the reverse problem of image reconstruction from sensor
data.
