In this paper, we consider a weak Galerkin finite element method for the Kelvin-Voigt viscoelastic fluid flow model. Firstly, the weak Galerkin finite element method is used to approximate the spatial variable and we use piecewise polynomials of degrees k, k - 1and k - 1(k >= 1) to approximate the velocity, pressure, and the numerical trace of the velocity on the interfaces of elements, respectively. Secondly, the backward Euler difference method is adopted in the temporal discretization for the fully discrete scheme. Furthermore, the stability and optimal convergence of numerical solutions in L-infinity (L-2) and L-infinity (H-1)-norms of velocity as well as L-infinity (L-2)-norm of pressure were presented. Finally, numerical examples verify the effectiveness of the proposed method, which also obtain that the algorithm has convergence and robust for different retardation time. (c) 2022 IMACS. Published by Elsevier B.V. All rights reserved.