Let $\Omega$ be a bounded domain, and consider the equation $$u_t - \Delta u = f$$ $$u(0) = u_0 \in L^1(\Omega)$$ with Neumann BCs (or Dirichlet if convenient) where $f$ is smooth.
Using energy/variational methods, is it true that there is a solution $u$ in some Bochner space? What spaces does the solution lie in? What can I expect? I would like something like $u \in L^2(0,T;H^1)$ with $u' \in L^2(0,T;H^{-1})$ which is standard for $u_0 \in L^2$ but not sure about the $L^1$ case. In particular, I would like some regularity on the gradient and the time derivative (so I do not want "very weak solutions").
I am well-versed in the theory for $u_0 \in L^2$ and am not interested in that case. Please note that I want to do this in a variational setting so semigroup/classical methods are out of the questions. Thanks for any help.