We investigate the behavior of the time derivatives of the solution to a linear timefractional, 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.