The aim of the present work was the modeling of the wastewater treatment plant operation work using Monte Carlo method and different random variables probability distributions modeling. The analysis includes the following pollutants indicators; BOD5 (Biochemical Oxygen Demand), CODCr (Chemical Oxygen Demand), Total Suspended Solids (SSt), Total Nitrogen (TN), and Total Phosphorus (TP). The Anderson–Darling (A–D) test was used for the assessment of theoretical and empirical distributions compatibility. The selection of the best-fitting statistical distributions was performed using peak-weighted root mean square (PWRMSE) parameter. Based on the performed calculations, it was stated that pollutants indicators in treated sewage were characterized by a significant variability. Obtained results indicate that the best-fitting pollutants indicators statistical distribution is Gauss Mixed Model (GMM) function. The results of the Monte Carlo simulation method confirmed that some problems related to the organic and biogenic pollutants reduction may be observed in the Wastewater Treatment Plant, in Jaworzno.