The objective of this study was to model the operation of a vertical-flow constructed wetland (VF-CW) for domestic wastewater, using Monte-Carlo simulations and selected probability distributions of various random variables. The analysis was based on collected wastewater quality data, including the values of the pollutant indicators BOD5 (biochemical oxygen demand), CODCr (chemical oxygen demand), and TSS (total suspended solids), in the 2017–2020 period. Anderson–Darling (A–D) statistics were applied to assess the fit of the theoretical distributions to the empirical distributions of the random variables under study. The selection of the best-fitting statistical distributions was determined using the percentage deviation (PBIAS) criterion. Based on the analyses that were performed, the best-fitting statistical distributions for the pollution indicators of the raw wastewater were the generalised extreme value distribution for BOD5, the Gaussian distribution for CODCr, and the log-normal distribution for TSS. For treated effluent, the log-normal distribution was the best fit for BOD5 and CODCr; the semi-normal distribution, for TSS. The new data generated using the Monte-Carlo method allowed the reliability of the VF-CW operation to be assessed by determining the reliability indices, i.e., the average efficiency of the removal of pollutants (η), the technological efficiency index (R), the reliability index (CR), and the risk index of the negative control of the sewage treatment plant operation (Re). The obtained results indicate that only in the case of CODCr, the analysed treatment facility may fail to meet the requirements related to the reduction of organic pollutants to the required level, which is evidenced by the values of the indicators CR = 1.10, R = 0.49, and η = 0.82. In addition, the risk index of the negative operation of the facility (Re) assumes a value of 1, which indicates that during the period of its operation, the VF-CW system will not operate with the required efficiency in relation to this indicator. The novelty of this work is the implementation of the indicated mathematical simulation methods for analysing the reliability of the operation of the domestic wastewater treatment facility.