SUMMARYTo obtain the deformation of the solid Earth from a global navigation satellite system (GNSS)-observed velocity field, the loading effect of the surface mass variations should be effectively deducted. However, the GNSS-observed velocity field in mainland China is currently limited only to the loading correction calculated using Gravity Recovery and Climate Experiment (GRACE) spherical harmonic coefficients, which is equivalent to the approximately 300-km smoothed result in the spatial domain; thus, the derived tectonic deformation is inaccurate. Therefore, it is important to study and identify a reasonable method for calculating the loading effect of the surface mass change model and to carry out an effective loading correction of the GNSS velocity field. In this study, the performances of two calculation methods, namely the GRACE spherical harmonic coefficient and Green's function, were analyzed and compared. In addition, we constructed a comprehensive model of the global surface mass variations, calculated the vertical load velocity in mainland China using Green's function method and compared the results with those for the GRACE spherical harmonic products. We found that the difference between the results of the GRACE spherical harmonic coefficient and Green's function methods was more than 1 mm/yr in the North China Plain, implying that the GRACE spherical harmonic coefficient method cannot be used for loading correction of the observed GNSS vertical velocity field. In contrast, the loading effect calculated using Green's function method can be more effectively applied for loading correction of the GNSS vertical velocity field in mainland China. The GNSS-observed velocity exhibited a clear uplift in the North China Plain and the west glacier areas; however, the GNSS velocity fields were significantly reduced after the loading correction, indicating that the observed GNSS vertical velocity fields were mainly caused by the surface mass loading due to the negative correlation between the vertical load velocity and the surface mass changes. Moreover, we found that the loading correction accounted for more than 50 per cent of the GNSS vertical velocity field in most of the glaciated regions in eastern and western China, and the maximum value exceeded 300 per cent, indicating that the loading effect was large. Finally, we obtained the GNSS vertical velocity field for mainland China with a loading correction. Additionally, the spectral characteristics of the time-varying gravity field in mainland China were investigated. The results showed that clear annual, semi-annual and 10-year medium- and long-period signals exist.