Abstract. Earth's ionosphere is an important medium of radio wave propagation in modern times. However, the effective use of the ionosphere depends on the understanding of its spatiotemporal variability. Towards this end, a number of ground- and space-based monitoring facilities have been set up over the years. The information from these stations has also been complemented by model-based studies. However, assessment of the performance of ionospheric models in capturing observations needs to be conducted. In this work, the performance of the IRI-2016 model in simulating the total electron content (TEC) observed by a network of Global Positioning System (GPS) receivers is evaluated based on the RMSE, the bias, the mean absolute error (MAE) and skill score, the normalized mean bias factor (NMBF), the normalized mean absolute error factor (NMAEF), the correlation, and categorical metrics such as the quantile probability of detection (QPOD), the quantile categorical miss (QCM), and the quantile critical success index (QCSI). The IRI-2016 model simulations are evaluated against gridded International Global Navigation Satellite System (GNSS) Service (IGS) GPS-TEC and TEC observations at a network of GPS receiver stations during the solar minima in 2008 and solar maxima in 2013. The phases of modeled and simulated TEC time series agree strongly over most of the globe, as indicated by a high correlations during all solar activities with the exception of the polar regions. In addition, lower RMSE, MAE, and bias values are observed between the modeled and measured TEC values during the solar minima than during the solar maxima from both sets of observations. The model performance is also found to vary with season, longitude, solar zenith angle, and magnetic local time. These variations in the model skill arise from differences between seasons with respect to solar irradiance, the direction of neutral meridional winds, neutral composition, and the longitudinal dependence of tidally induced wave number four structures. Moreover, the variation in model performance as a function of solar zenith angle and magnetic local time might be linked to the accuracy of the ionospheric parameters used to characterize both the bottom- and topside ionospheres. However, when the NMBF and NMAEF are applied to the data sets from the two distinct solar activity periods, the difference in the skill of the model during the two periods decreases, suggesting that the traditional model evaluation metrics exaggerate the difference in model skill. Moreover, the performance of the model in capturing the highest ends of extreme values over the geomagnetic equator, midlatitudes, and high latitudes is poor, as noted from the decrease in the QPOD and QCSI as well as an increase in the QCM over most of the globe with an increase in the threshold percentile TEC values from 10 % to 90 % during both the solar minimum and the solar maximum periods. The performance of IRI-2016 in simulating observed low (as low as the 10th percentile) and high (higher than the 90th percentile) TEC correctly over equatorial ionization anomaly (EIA) crest regions is reasonably good given that IRI-2016 is a climatological model. However, it is worth noting that the performance of the IRI-2016 model is relatively poor in 2013 compared with 2008 at the highest ends of the TEC distribution. Therefore, this study reveals the strengths and weaknesses of the IRI-2016 model in simulating the observed TEC distribution correctly during all seasons and solar activities for the first time.