In this paper, we investigate a variational discretization approximation of parabolic bilinear optimal control problems with control constraints. For the state and co-state variables, triangular linear finite element and difference methods are used for space and time discretization, respectively, superconvergence in H 1 -norm between the numerical solutions and elliptic projections are derived. Although the control variable is not discrete directly, convergence of second order in L 2 -norm is obtained. These theoretical results are confirmed by two numerical examples.