This paper aims to study an H1-Galerkin mixed finite element method (FEM) for Klein–Gordon–Zakharov (KGZ) equations with power law nonlinearity. The bilinear element and zero-order Raviart–Thomas element (Q11/Q10×Q01) are taken as approximation spaces for the original variables p, ϕ and the auxiliary variables u=∇p and ψ=∇ϕ, respectively. The supercloseness and superconvergence estimates of order O(h2+τ2) are presented for p and ϕ in the H1 norm and u and ψ in the L2 norm, which improve the known results that the original variable ϕ can only reach optimal error estimate shown in numerical experiments without theoretical analysis. Here h is the subdivision parameter and τ is the time step. Finally, some numerical results are provided to support the theoretical analysis.