This should be a comment but I do not have the right to.

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}\|_{L^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}\|_{L^q} + \|f_n - f_m\|_{L^p(L^q)}\right),$$
and thus you have a Cauchy sequence.

The problem is that this is not the right space for the initial condition. To have this kind of estimate you need $u_{0,n}$ in some Besov (don't ask me, don't know them, look in the link).

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