The aim of the present paper is to evaluate the capabilities of the unsteady Reynolds-averaged Navier-Stokes (URANS) method to simulate a turbulent plume. Turbulent plume is mainly driven by thermal buoyancy and is characterized by unsteadiness, energy nonequilibrium, counter-gradient diffusion and lack of universal scaling. All these phenomena are believed to be associated with distinct large coherent eddy structures, which are hardly captured by steady RANS equations. In URANS approach, the large scale structures are fully resolved, whereas the stochastic small scales are modeled. Therefore large coherent eddy structure can be captured well by the URANS computations. The present simulation is carried out using an open source code called buoyantPimpleFOAM which is based on OpenFOAM platform. For modeling the effect of turbulence, the buoyancy-corrected standard k-epsilon model and RNG k-epsilon model are employed. Since buoyancy plays a significant role in production and destruction of turbulence in a plume, additional turbulence production term due to buoyancy effects is introduced on the basis of simple gradient diffusion hypothesis (SGDH) and generalized gradient diffusion hypothesis (GGDH). A comparison of computational results with the experimental data reported in the literature has been presented in detail. We are able to capture the well-known Gaussian nature of the plume and empirical plume spread coefficients. (c) 2014 Elsevier Ltd. All rights reserved.