  A systematic numerical study on weak Galerkin (WG) finite element method for
second order linear parabolic problems is presented by allowing polynomial
approximations with various degrees for each local element. Convergence of both
semidiscrete and fully discrete WG solutions are established in
$L^{\infty}(L^2)$ and $L^{\infty}(H^1)$ norms for a general WG element $({\cal
P}_{k}(K),\;{\cal P}_{j}(\partial K),\;\big[{\cal P}_{l}(K)\big]^2)$, where
$k\ge 1$, $j\ge 0$ and $l\ge 0$ are arbitrary integers. The fully discrete
space-time discretization is based on a first order in time Euler scheme. Our
results are intended to extend the numerical analysis of WG methods for
elliptic problems [J. Sci. Comput., 74 (2018), 1369-1396] to parabolic
problems. Numerical experiments are reported to justify the robustness,
reliability and accuracy of the WG finite element method.
