This paper presents the results of research and analyses related to the development of a new quasigeoid model fitted to GNSS/levelling data for the area of Poland (PL-geoid2021). The model was determined employing two procedures based on the Geophysical Gravity data Inversion technique (GGI method): procedure A consisted of the determination of the gravimetric quasigeoid model in the first step and its subsequent fitting to GNSS/levelling data in the second step, and procedure B consisted of a one-step determination of the model fitted to GNSS/levelling data. Both models were developed using the global geopotential model SGG-UGM-2 and gravity data covering the area of Poland, and slightly extend beyond Poland's southern and northern borders. The average model was adopted as the final model. It was demonstrated that the accuracy of the gravimetric quasigeoid model had a very low dependence on the reference topographic mass density model used. On the basis of this model, the GNSS/levelling datasets were also assessed and outliers were identified. The estimated accuracy of the gravimetric model, determined based on four GNSS/levelling datasets, was in the range of ± 1.2 to ± 1.7 cm, in terms of the standard deviation of the differences between the measured and model-determined height anomalies. Due to partial lack of gravity data just beyond the Polish border, the edge effect was also analysed. The accuracy of the final quasigeoid model (estimated in the same way as the gravimetric model) ranges from ± 1.0 to ± 1.2 cm. It should be noted, however, that this assessment is not fully independent because three of the four sets of GNSS/levelling points used for it, were also used to build the final model.