Global Navigation Satellite Systems (GNSS) can sense deformations of the Earth’s crust. All components, but in particular the vertical component are affected by mass loading, i.e. external forces resulting from the redistribution and changes in fluid mass. These effects include non-tidal atmospheric loading (NTAL), non-tidal ocean loading (NTOL), and hydrological loading (HYDL). If these loading effects are not compensated in the processing of space geodetic data, the obtained results will be distorted. Thus, physics-based loading models exist that can be applied to correct station positions. This study investigates if machine learning (ML) in combination with environmental variables can replace or augment the existing physics-based models via a data-driven modelling of GNSS displacements. Therefore, vertical displacements of 3553 GNSS stations in Europe are utilized to train and validate XGBoost models. Three different strategies were tested, differing in the preprocessing of the GNSS data, i.e. whether or which physics-based loading models were applied beforehand. A significant improvement was achieved for all strategies ranging from 4.4% to 22.9%. The improvement is calculated based on the root mean squared error (RMSE) reduction of the GNSS residual coordinates w.r.t. a trajectory model, accounting for a linear trend, seasonal signals, and discontinuities in the GNSS coordinate time series. In addition to evaluating the ML models, a thorough feature importance analysis based on SHapley Additive exPlanations (SHAP) is carried out to better understand the driving factors of the model output and to gain insights into what signals could still be found to enhance existing physical models.