This paper concerns with the finite element approximation of a nonlinear second order parabolic system which describes the $L^2$-gradient flow for a class of linear growth energy functionals. Besides their appeals in differential geometry and calculus of variations, linear growth energy functionals and their gradient flows also arise naturally from emerging applications of image processing such as color image denoising. In this paper, we introduce a family of variational models for color image denoising which minimize linear growth energy functionals of maps into the unit sphere in $\rm{R}^3$. These models generalize the popular 1-harmonic map model which has been studied intensively in recent years. To compute the solutions of the variational models, we first derive their $L^2$-gradient flow equations and then introduce some fully discrete implicit finite element method for the gradient flow equations. It is proved that the proposed finite element method is uniquely solvable and absolutely stable, and the finite element solution converges to the PDE solution as the mesh sizes tend to zero. Numerical experiments are presented to demonstrate the effectiveness of the proposed variational models for color image denoising and to show the efficiency of the proposed finite element method. A numerical comparison of the proposed models with the channel-by-channel model is also presented.