The surface displacement caused by hydrological loading makes an important contribution to the non-linear crustal movement observed at the International Global Navigation Satellite System Service (IGS) stations. In this paper, the amplitude, correlation, and root mean square (RMS) of the vertical displacement time series signals of 47 IGS stations are used to analyze which data of Gravity Recovery and Climate Experiment (GRACE) or Global Land Data Assimilation System (GLDAS) can better reflect the hydrological load effect in Europe. The results show that in Europe, the hydrological load effect calculated based on GRACE data is more accurate than that of GLDAS, which has not been reported before. Then, the relationship between the GPS height and GRACE load deformation in terms of annually-oscillating signals, correlation, and phase is analyzed by using singular spectrum analysis, the Pearson correlation coefficient, and wavelet coherence (WTC). It was found that GPS and GRACE agree at some stations (e.g., BOR1 and ZIMM), while they differ significantly in amplitude and phase at other stations (e.g., KIRU and NOT1), indicating that not all GRACE-derived displacements of IGS stations can clearly explain their nonlinear motion. and The correlation coefficients between GPS and GRACE are higher than 0.7 at 85 % of stations. Amongst them, the values are obviously greater than 0.8 (e.g., ZIMM and LAMA) around inland areas and high mountains, and even less than 0.6 (e.g., ANKR and KIRU) along the coast of the Mediterranean ocean, which more precisely shows that the hydrological load effect has obvious spatial and regional characteristics compared with previous studies. In addition, the relative phase of the WTC solution is basically consistent under non-detrend and detrend, which shows that the relative phase difference of each station is only related to the nonlinear movement and not to the linear trend caused by the tectonic deformation. Finally, we study the influence of GRACE hydrological load on the RMS of GPS height, which is reduced by 24.60 % on average, and the reduction rate distribution of the RMS is in good agreement with the spatial distribution of the correlation coefficient.