  We investigate the behavior of the time derivatives of the solution to a
linear time-fractional, advection-diffusion-reaction equation, allowing space-
and time-dependent coefficients as well as initial data that may have low
regularity. Our focus is on proving estimates that are needed for the error
analysis of numerical methods. The nonlocal nature of the fractional derivative
creates substantial difficulties compared with the case of a classical
parabolic PDE. In our analysis, we rely on novel energy methods in combination
with a fractional Gronwall inequality and certain properties of fractional
integrals.
