<p>This article develops the primal hybrid finite element method with Lagrange multipliers to approximate nonlinear parabolic initial-boundary value problems with gradient type non-linearity. A modified elliptic projection is used to produce optimal order error estimates for the semi-discrete and backward Euler-based complete discrete schemes. In addition, error estimates in <i>L</i>^∞-norm are established which are optimal in nature. Superconvergence result of the gradient in <i>L</i>^∞-norm is discussed for the error between the primal hybrid solution and elliptic projection. As a bi-product, the proposed analysis provides optimal error analysis for non-conforming CR-elements. Finally, numerical tests are performed to validate the theoretical findings.</p>