With the significant increase in WQI applications worldwide and lack of specific application guidelines, accuracy and reliability of WQI models is a major issue. It has been reported that WQI models produce significant uncertainties during the various stages of their application including: (i) water quality indicator selection, (ii) sub-index (SI) calculation, (iii) water quality indicator weighting and (iv) aggregation of sub-indices to calculate the overall index. This research provides a robust statistically sound methodology for assessment of WQI model uncertainties. Eight WQI models are considered. The Monte Carlo simulation (MCS) technique was applied to estimate model uncertainty, while the Gaussian Process Regression (GPR) algorithm was utilised to predict uncertainties in the WQI models at each sampling site. The sub-index functions were found to contribute to considerable uncertainty and hence affect the model reliability – they contributed 12.86% and 10.27% of uncertainty for summer and winter applications, respectively. Therefore, the selection of sub-index function needs to be made with care. A low uncertainty of less than 1% was produced by the water quality indicator selection and weighting processes. Significant statistical differences were found between various aggregation functions. The weighted quadratic mean (WQM) function was found to provide a plausible assessment of water quality of coastal waters at reduced uncertainty levels. The findings of this study also suggest that the unweighted root means squared (RMS) aggregation function could be potentially also used for assessment of coastal water quality. Findings from this research could inform a range of stakeholders including decision-makers, researchers, and agencies responsible for water quality monitoring, assessment and management.