The ground-based GNSS VTEC model can adequately capture the spatiotemporal characteristics of ionosphere activities. However, it is difficult to precisely model VTECs with a unified mathematical model. We dedicate to comprehensively assess the performance of four mathematical models for VTEC modeling, i.e., polynomial, trigonometric series, spherical harmonic and multi-surface function. To capture the varying ionosphere situations, three typical regions in Western Europe, Southeast China, and North America are chosen. To reflect the precision and accuracy from different aspects, four evaluation measures are defined based on the comparison with CODE GIM products, the poster-fitting residuals of VTEC modeling, the comparison with high-precision quasi-VTEC computed with L1–L2 phase observations as well as the analysis of precision and stability of receiver and satellite DCB by-products. The results show that, in small regions, the performances of all models are insensitive to the model orders but sensitive to the ionosphere activity. On the whole, the polynomial and spherical harmonic function are comparable in terms of their performance and computation efficiency; Trigonometric series is instable with systematic biases in large region due to its incapability of describing the spatial ionosphere variations; Multi-surface function outperforms the others thanks to its epoch-wise solution with low computation efficiency. In addition, the annual solutions of satellite and receiver DCBs indicate their sufficient stability, which can then be applied as known for single-epoch ionosphere modeling.