This work is devoted to an inverse problem of identifying a source term depending on both spatial and time variables in a parabolic equation from single Cauchy data on a part of the boundary. A Crank-Nicolson Galerkin method is applied to the least squares functional with quadratic stabilizing penalty term. The convergence of finite dimensional regularized approximations to the sought source as measurement noise levels and mesh sizes approach to zero with appropriate regularization parameter is proved. Moreover, under a suitable source condition, an error bound and corresponding convergence rates are proved. Finally, two numerical experiments are presented to illustrate the theoretical findings.