The capability to predict the distribution of pollutants in water bodies is one of the most important issues in the design of jet outfalls. Three-dimensional computational fluid dynamics (CFD) model and multi-objective evolutionary polynomial regression (EPR-MOGA) are used and compared in modeling the temperature field in the side thermal buoyant discharge in the cross flow. The input variables used for training the EPR-MOGA models are spatial coordinates (x, y, z), jet to cross flow velocity ratio (R), depth of the channel (d), and the temperature excess (T0). A previous experimental study is used to verify and compare the performance of the EPR-MOGA and CFD models. The results show that the EPR-MOGA model predicts the thermal cross section of the flow and the spread of pollutants at the surface with a better accuracy than the CFD model. However, the CFD method performs significantly better than EPR-MOGA in predicting temperature profiles. The uncertainty analysis indicated that the EPR-MOGA model had lower mean prediction error and smaller uncertainty band than the CFD model. The relationships achieved by the EPR-MOGA model are very useful to predict temperature profiles, temperature half-thickness, and temperature spread on surface in practice.