The prevalent multifractal characteristics of turbulent velocity fluctuations in the atmosphere are important for estimating various wind effects in wind engineering. Here, the multifractal characteristics of turbulent velocity fluctuations, including the small-scale multiscaling, the long-tail distributions and the intermittency, are thoroughly investigated by using a high-frequency dataset of three-dimensional velocities (100 Hz) collected at three levels during one month. To reduce uncertainties in the estimate of multiscaling exponents, a new method, the sequential extended self-similarity, is proposed. Based on this method, we obtain the multiscaling exponents of qth-order moments of velocity increments as a function of q, that is the so-called multifractal spectrum. The multifractal random walk (MRW) model is then shown to describe the various multifractal spectra of turbulent velocity fluctuations. With the help of this model, we find a connection between the small-scale multiscaling and the long-tail distributions, which is generally observed in our dataset, again validating the MRW model. A non-linear multifractal spectrum is commonly considered to be related to the intermittency of turbulent velocity fluctuations at small scales and its curvature is usually used as a quantification of intermittency intensity. However, we suggest that models capturing the non-linear multifractal spectrum may fail to represent the long-tail distribution, which is a more direct quantification of intermittency. Finally, qualitative variations of validated indicators with specific boundary-layer parameters are investigated. Results show that the intermittency of turbulent velocity fluctuations is more relevant to the friction velocity, compared with the average wind speed, the average temperature, and the surface-layer stability.