I have the equation $$u_t - \Delta u = f\text{ on $\Omega$}$$ $$\partial_\nu u = g\text{ on $\partial\Omega$}$$ $$u(0) = u_0$$ for $f \in L^2(0,T;H^1)$, $g \in L^2(0,T;H^1(\partial\Omega))$ and $u_0 \in L^\infty$.
Is it possible to represent the solution $u$ in terms of semigroups and the Neumann heat kernel? Eg. the weak solution $u \in H^1(0,T;H^{-1}) \cap L^2(0,T;H^1)$ can be written as $$u(t) = A(t)u_0 + B(t)g(t) + C(t)f(t)$$,
I couldn't find anything after searching online... all seem to discuss only the Dirichlet case.
The usual trick when you have boundary data is to get rid of it by shifting it to the inside of your domain. For instance in your case, you would want to find a function $$G : \Omega \times (0,T) \to \mathbb{R}$$ such that $$G(x,t) = g(x,t)$$ for $$x \in \partial \Omega$$. Then instead of your problem you consider $$ \begin{cases} v_t - \Delta v = f - G_t + \Delta G &\text{in }\Omega \\ \partial_\nu v =0 &\text{on } \partial \Omega, \end{cases} $$ and you attack this with semigroup methods. Then $$u := v+ G$$ should solve your problem.
To actually pull this off you might need to tweak the regularity you have for $g$. It seems a bit off from what I would expect for a parabolic problem of this type.
One more remark: I think your formula for the semigroup solution is off in the homogeneous boundary data case. It should read $$ u = A(t) u_0 + \int_0^T A(t-s) f(s)ds, $$ which is the usual "Duhamel formula."