Numerical investigation of a side thermal buoyant discharge in the cross-flow is presented based on unsteady Reynolds-averaged Navier–Stokes equations closed with the realizable k- $$\upvarepsilon $$ turbulence model. Emphasis is placed on the detailed three-dimensional flow evolution and scalar mixing in an incompressible turbulent environment. The present study covers the cases with different jets to cross-flow velocity ratios (R) and initial temperatures. Moreover, various flow characteristics, including vortical structures, jet trajectories, jet streamlines, and intrinsic instabilities, are examined. Mixing ability is quantified by the decay rate of scalar temperatures and velocity magnitude, the probability density function, the spatial mixing deficiencies (SMDs), power spectral density analysis, and temporal mixing deficiencies (TMDs). Comparing the simulation results with the experimental data of Abdelwahed (Surface jets and surface plumes in cross-flows. Ph.D. thesis, Department of Civil Engineering and Applied Mechanics, McGill University, Montreal, Quebec, 1981), it shows good agreement in terms of the temperature half-width and half-thickness profiles, recirculation zone size, and jet trajectory. The numerical results of three-dimensional structures indicate a shear-layer, horseshoe, and surface roller vortices near the side-channel exit and secondary flows after the recirculation zone at the free surface of the main channel. The instantaneous temperature contours exhibit a vortex shedding phenomenon and gaps between the cross-flow and discharge jet at the shear layer. The maximum velocity magnitude location approaches the outer wall toward the main-channel downstream by increasing R. It is found that as the densimetric Froude number ( $$\hbox {Fr}_{{0}})$$ and R increase, the temperature dilution (S) generally decreases. The statistical analysis of TMD and SMD indicates a direct relationship between the mixing efficiency and buoyancy flux ( $${F}_{{0}})$$ .