AbstractA theory of single-crystal gradient thermoplasticity that accounts for the stored energy of cold work and thermal annealing has recently been proposed by Anand et al. (Int J Plasticity 64:1–25, 2015). Aspects of the numerical implementation of the aforementioned theory using the finite element method are detailed in this presentation. To facilitate the implementation, a viscoplastic regularization of the plastic evolution equations is performed. The weak form of the governing equations and their time-discrete counterparts are derived. The theory is then elucidated via a series of three-dimensional numerical examples where particular emphasis is placed on the role of the defect-flow relations. These relations govern the evolution of a measure of the glide and geometrically necessary dislocation densities which is associated with the stored energy of cold work.