We are concerned with the Neumann type boundary value problem to parabolic systems ∂ t u − div(D ξ f (x, Du)) = −D u g(x, u), where u is vector-valued, f satisfies a linear growth condition and ξ → f (x, ξ) is convex. We prove that variational solutions of such systems can be approximated by variational solutions to ∂ t u − div(D ξ f p (x, Du)) = −D u g(x, u) with p > 1. This can be interpreted both as a stability and existence result for general flows with linear growth.