Edited to make it correct for future reference.

You might want to look into this expository paper on maximal regularity for linear parabolic equations:
https://people.math.ethz.ch/~salamon/PREPRINTS/parabolic.pdf

Then, if you had an estimate of the form
$$\|u_n\|_{W^{1,p}_t(L^q_x)} \le C\left(\|u_{n,0}\|_{W^{2,q}} + \|f_n\|_{L^p(L^q)}\right),$$

for the same price, using that the equation is linear, you could use that
$$\|u_n - u_m\|_{W^{1,p}_t(L^q_x)} \le C\left(\|u_{n,0} - u_{m,0}\|_{W^{2,q}} + \|f_n - f_m\|_{L^p(L^q)}\right),$$
and thus you have a Cauchy sequence.

Note that $W^{2,q}$ is enough initial regularity to carry over, the optimal space is a Besov (check it in the link), but you can have the same kind of estimate with a more restrictive class.

I don't think you are going to have a stability/convergence result without having some existence + quantitative estimate.