Evaluating the distribution patterns of surface temperature data at Very Long Baseline Interferometry (VLBI)/Global Positioning System (GPS) co-located sites w.r.t. normality is one of the most important issues in modeling surface temperature data over long periods. Such evaluation can generate algorithms for filling in missing data at measurement sites. Some algorithms in the literature, such as those in the study of Cho et al. J Coast Res 65. doi: 10. 2112/SI65-321. 1, (2013), require trend, harmonic, and residual components to fill in the missing data. Trend and harmonic components estimate an optimal model that can be used to assist such algorithms when filling in missing data. The present study is based on the investigation of the normal distribution of the residuals of a surface temperature time series at VLBI/GPS co-located sites, after removing the trend and seasonal effects through harmonic components (inter-daily variations). This study uses surface temperature data collected from the VLBI/GPS co-located sites of two different regions in Europe: Matera (Italy) and Wettzell (Germany). The data collected from these sites form a time series, and time series analyses and conventional k-sigma outlier detection are implemented on these data sets before subjecting them to goodness of fit tests for normality. The residual components of the time series are acquired through a decomposing trend and signal effect from the original time series, assuming that the residuals of the time series are normally distributed. In testing the hypothesis that an observed frequency distribution fits the normal distribution, the following tests are used: Pearson chi (2), Kolmogorov-Smirnov, Anderson-Darling, Shapiro-Wilk or Shapiro-Francia, D'Agostino, Jarque-Bera, skewness, and kurtosis tests. Some graphical methods are also applied to support the results of the goodness of fit tests for normality. Some proposals on the application of the goodness of fit tests are put forward, such as the evaluation of the estimation model for trend and harmonic components by considering the properties of the implemented goodness of fit tests. The results of this study can be used to determine the optimal model for estimating trend and harmonic components. The output of the present study is expected to have an important role in modeling surface temperature distributions at co-located VLBI/GPS sites for filling in missing data. Above all, meteorological data, such as temperature, pressure, and humidity, are of specific interest for modeling tropospheric delay, the main error factor in positioning in space geodesy, which in turn makes investigations on the distribution of meteorological data more attractive in geoscience.