Epidemic threshold has always been a very hot topic for studying epidemic dynamics on complex networks. The previous studies have provided different theoretical predictions of the epidemic threshold for the susceptible-infected-recovered (SIR) model, but the numerical verification of these theoretical predictions is still lacking. Considering that the large fluctuation of the outbreak size occurs near the epidemic threshold, we propose a novel numerical identification method of SIR epidemic threshold by analyzing the peak of the epidemic variability. Extensive experiments on synthetic and real-world networks demonstrate that the variability measure can successfully give the numerical threshold for the SIR model. The heterogeneous mean-field prediction agrees very well with the numerical threshold, except the case that the networks are disassortative, in which the quenched mean-field prediction is relatively close to the numerical threshold. Moreover, the numerical method presented is also suitable for the susceptible-infected-susceptible model. This work helps to verify the theoretical analysis of epidemic threshold and would promote further studies on the phase transition of epidemic dynamics.