PurposeImportant statistical variations are likely to appear in the propagation of surface plasmon polariton waves atop the surface of graphene sheets, degrading the expected performance of real-life THz applications. This paper aims to introduce an efficient numerical algorithm that is able to accurately and rapidly predict the influence of material-based uncertainties for diverse graphene configurations.Design/methodology/approachInitially, the surface conductivity of graphene is described at the far infrared spectrum and the uncertainties of its main parameters, namely, the chemical potential and the relaxation time, on the propagation properties of the surface waves are investigated, unveiling a considerable impact. Furthermore, the demanding two-dimensional material is numerically modeled as a surface boundary through a frequency-dependent finite-difference time-domain scheme, while a robust stochastic realization is accordingly developed.FindingsThe mean value and standard deviation of the propagating surface waves are extracted through a single-pass simulation in contrast to the laborious Monte Carlo technique, proving the accomplished high efficiency. Moreover, numerical results, including graphene’s surface current density and electric field distribution, indicate the notable precision, stability and convergence of the new graphene-based stochastic time-domain method in terms of the mean value and the order of magnitude of the standard deviation.Originality/valueThe combined uncertainties of the main parameters in graphene layers are modeled through a high-performance stochastic numerical algorithm, based on the finite-difference time-domain method. The significant accuracy of the numerical results, compared to the cumbersome Monte Carlo analysis, renders the featured technique a flexible computational tool that is able to enhance the design of graphene THz devices due to the uncertainty prediction.