In this paper we introduce a purely variational approach to the gradient flows, naturally arising in image denoising models, yielding the existence of global parabolic minimizers, in the sense that $ \int_0^{T}\big[\int_\Omega u \partial_t\varphi\, dx + F(u)\big] \,dt \le \int_0^{T}F(u+\varphi) \,dt, $ whenever $T>0$ and $\varphi\in C_0^\infty(\Omega\times(0,T))$. Our method applies to a wide class of nonparametric regression models in image restoration analysis, such as quantile, robust, and logistic regression. A prototype functional $F$ is the by now classical ${\rm TV}(L^2)$-functional (i.e., the pure ${\rm TV}$-denoising case in image reconstruction) introduced by Rudin, Osher and Fatemi [Phys. D, 268 (1992), pp. 259--268]: $ \boldsymbol F(u):= \mathbf{TV}(u)+\frac{\kappa}{2}\int_\Omega |u-u_o|^2\, dx, $ where $u_o\colon\Omega\to [0,1]$ is a noisy, monochromatic image and $\kappa\gg 1$ a large penalization parameter. The evolutionary variational solutions are obtained as limits of maps, minimizing a convex variational functional in n+1 dimensions with domain $\Omega_T:=\Omega \times (0,T)$. Our approach yields a new way of proving the existence of global weak solutions to the associated Cauchy--Dirichlet problem, $\partial _{t}u- {\rm div} \big(\frac{D u}{|D u|}\big)=\kappa (u-u_o)$ in $\Omega \times (0,\infty)$ and $u=u_o$ on the parabolic boundary. Our approach also applies in situations where the considered functionals do not allow the derivation of the associated parabolic equation. We are able to deal with Dirichlet and Neumann type boundary conditions on the lateral boundary, and furthermore with the gradient flow associated to functionals modeling image deblurring, such as $ F(u)=\mathbf{TV}(u)+\frac\kappa2\int_\Omega\big| \mathbf K[u]-u_o\big|^2\, dx, $ where $\mathbf K\colon L^1(\Omega) \to L^2(\Omega)$ is a bounded, linear, injective operator satisfying the DC-condition $\mathbf K[1]=1$.