Temporal variations of geoid heights are vitally important in geodesy and Earth science. They are essentially needed for dynamic and kinematic updates of the static geoid model. These temporal variations, which substantially differ for different geographic locations, can successfully be determined using the Gravity Recovery And Climate Experiment (GRACE) mission data. So far, statistical decomposition methods, e.g. the Principal Component Analysis/Empirical Orthogonal Function (PCA/EOF) method, have not been implemented for the analysis and modelling of temporal mass variations within the Earth's system over the area of Poland. The aim of this contribution is to analyse and model temporal variations of geoid heights obtained from GRACE mission data over the area of Poland using the PCA/EOF method. Temporal variations of geoid heights were obtained from the latest release, i.e. release five, of monthly GRACE-based Global Geopotential Models. They can reach the level of 10 mm. The PCA modes and their corresponding EOF loading patterns were estimated using two different algorithms. The results obtained revealed that significant part of the signal of temporal variations of geoid heights over Poland can be obtained from the first three PCA modes and EOF loading patterns. They demonstrate the suitability of the PCA/EOF method for analysing and modelling temporal variations of geoid heights over the area investigated.