Terahertz pulse time-domain holography is the ultimate technique allowing the evaluating a propagation of pulse broadband terahertz wavefronts and analyze their spatial, temporal and spectral evolution. We have numerically analyzed pulsed broadband terahertz Gauss-Bessel beam’s both spatio-temporal and spatio-spectral evolution in the non-paraxial approach. We have characterized two-dimensional spatio-temporal beam behavior and demonstrated all stages of pulse reshaping during the propagation, including X-shape pulse forming. The reshaping is also illustrated by the energy transfer dynamics, where the pulse energy flows from leading edge to trailing edge. This behavior illustrates strong spatio-temporal coupling effect when spatio-temporal distribution of Bessel beam’s wavefront depends on propagation distance. The spatio-temporal and spatio-spectral profiles for different spectral components clearly illustrate the model where the Bessel beam’s wavefront at the exit from the axicon can be divided into radial segments for which the wave vectors intersect. Phase velocity via propagation distance is estimated and compared with existing experimantal results. Results of the phase velocity calculation depend strongly on distance increment value, thus demonstrating superluminal or subluminal behavior.