Computational four-equation turbulence model is developed and is applied to predict twodimensional unsteady thermal surface discharge into a reservoir. Turbulent stresses and heat fluxes in the momentum and energy equations are determined from transport equations for the turbulent kinetic energy (R), isotropic rate of kinetic energy dissipation (.epsilon.), mean square temperature variance (theta. over bar $^{2}$), and rate of destruction of the temperature variance (.epsilon. $_{ heta}$). Computational results by four-equation model are favorably compared with those obtained by an extended two-equation model. Added advantage of the four-equation model is that it yields quantitative information about the ratio between the velocity time scale and the thermal time scale and more detailed information about turbulent structure. Predicted time scale ratio is within experimental observations by others. Although the mean velocity and temperature fields are similarly predicted by both models, it is found that the four-equation model is preferably candidate for prediction of highly buoyant turbulent flows.