change very little over typical integration spans (\(\sim\)1 mm), and models for non-tidal effects contain large uncertainty and are usually less accurate. Therefore, non-tidal mass variations are generally suggested by the International Earth Rotation Service (IERS) not to be applied in computing conventional instantaneous positions in GNSS routine data processing procedures [1, 11]. The corresponding non-tidal mass variations, which are typically with amplitudes of a few centimeters or less, still remain as signals embedded in the GNSS geodetic coordinate solutions [12, 13, 14]. Previous studies have established the significance of non-tidal loadings, especially the comprehensive effects induced by changes in atmospheric pressure loading (ATML), continental water storage (CWS), and non-tidal ocean loading (NTOL) [15, 16, 17]. The non-tidal loading effects can be modeled mainly by two different methods--mass variations from the Gravity Recovery and Climate Experiment (GRACE) satellite mission [18, 19, 20, 21] and solutions from geophysical models (GMs) [22, 23]. Many researchers have investigated the relationships among the different techniques, especially in regions with acute environmental mass variations [24, 25, 26, 27, 28]. The research of [25] used the homogeneously reprocessed 115 GNSS stations with high quality and GRACE level 2 gravity products from Release 04 from September 2002 to April 2009 for comparisons. Their GNSS data were completely reprocessed in the context of the GGOS-D project (Global Geodetic-Geophysical Observing System-Deutschland) [27], using the Bernese GNSS Software 5.1. The research of [26] compared eleven available solutions from the International GNSSS Service (IGS) Analysis Centers (ACs) with seven GRACE-modeled vertical deformation (VD) solutions based on the globally distributed 110 IGS08/IGS14 core stations, and suggested that GRACE data could provide VD results of good internal consistency, while GNSS solutions from multiple institutions implied large uncertainty. By investigating the predicted deformations with the position time series of 689 globally distributed continuous GNSS stations, research of [13] found that while GNSS vertical deformations are well explained by the model at a global scale, horizontal deformations were systematically underpredicted and out of phase with GNSS station position time series. The research of [13] also determined the estimates of degree-1 spherical harmonic coefficients from GRACE and global GNSS data, and demonstrated improved comparison of the computed deformations from their degree-1 solutions with GNSS measurements, particularly in the horizontal components. In terms of the comparisons between GNSS and GMs. The research of [28] comprehensively compared three different environmental loading methods to estimate surface deformations and correct the nonlinear variations in a set of GNSS weekly height time series. They demonstrated that there are differences remaining in different GMs, and pointed out that for different GMs, the consistency of the VD signals showed discrepancies. Although GNSS, GRACE and GMs can capture the seasonal vertical non-tidal deformations, there are still systematic discrepancies among the three techniques [29, 30, 31, 32]. It has been reported that the weighted root mean square (WRMS) of the GNSS residuals only decrease slightly after subtracting the GRACE-modeled or GM-derived VD series from the GNSS vertical coordinate time series, indicating non-tidal amplitudes differences existing in different techniques [25, 31]. Further analysis should be undertaken to analyze the relationships of these techniques and identify the possible reasons for the differences. GNSS data processing methods and theories have been under continuous refinement in the past 30 years [1, 11, 33, 34, 35]. The newest products come from the second data reprocessing campaign implemented by the IGS. Hereafter, we call it the IG2 products.