In this paper, we introduce a novel numerical approach combining the Crank–Nicolson time-stepping scheme with the weak Galerkin finite element method (WG-FEM) for solving a class of singularly perturbed, time-dependent convection–diffusion equations that include a nonlinear reaction component. To accurately resolve steep gradients and boundary layers induced by small perturbation parameters, we utilize a specially constructed layer-adapted mesh. The underlying solution structure and its asymptotic behavior with respect to the perturbation parameter \(\varepsilon\) are systematically analyzed, which informs the design of an effective interpolation technique tailored to the mesh. The convergence theory is established through a two-tier discretization framework, wherein spatial variables are first discretized using WG-FEM followed by temporal discretization via the Crank–Nicolson scheme. We derive uniform a priori error estimates in an energy-like norm, proving that the proposed method attains a convergence rate of order \(\mathcal{O}\big((\Delta t)^2 + N^{-r}\ln^{r}N\big)\) that is independent of \(\varepsilon\) . To substantiate the theoretical findings, we conduct a series of benchmark numerical experiments. The numerical investigations confirm that the method achieves the expected second-order accuracy in time and optimal convergence of order \(\mathcal{O}\big((\Delta t)^2 + N^{-(r+1)}\big)\) in the \(L^2\) -norm. Furthermore, the scheme demonstrates a superconvergent rate of \(\mathcal{O}\big((\Delta t)^2 + N^{-2r}\ln^{2r}N\big)\) in the discrete \(L^{\infty}\) -norm, underscoring the method’s robustness and precision in handling singular perturbations and nonlinearities simultaneously.
