ABSTRACT The present article numerically investigates the turbulent buoyancy-driven (natural convection) flow along a vertical plate with a low Reynolds turbulence two-equation k-ϵ model. The deployed turbulence model is appropriate for low Reynolds number (LRN) turbulent flow adjacent to a solid boundary, and the turbulence kinetic energy (TKE) and dissipation rate of TKE are estimated using the momentum equations and are solved simultaneously with the mean flow conservation equations. Two-dimensional time-dependent viscous incompressible turbulent flow is simulated. This flow domain is governed by a highly non-linear group of partial differential equations, namely the time-averaged continuity, momentum, and energy, and also the flow property is determined through TKE, and dissipation rate of TKE equations. Since these equations are not solvable using analytical methods, an implicit second-order finite difference method is employed to solve the governing turbulent flow equations numerically. The simulated time-averaged velocity, temperature, TKE, and dissipation rate of TKE profiles along with friction factor and heat transfer rate are computed for different values of turbulent Reynolds and Prandtl numbers. Average velocity, temperature, turbulence energy, and dissipation rate under both transient and steady-state conditions are decreased with increment in . There is a decrement in average transient velocity and turbulence energy rates with increasing values, whereas the average temperature profile increases. Average steady temperature and turbulence energy profiles are also enhanced with rising values of . The majority of researchers focus on laminar flow and thus far turbulent flow has mainly been addressed through experiments only. To better elucidate the flow in this condition, the novelty of the present study is to utilize a numerical FDM procedure to probe more deeply into the turbulence characteristics in buoyancy-driven (natural convection) flow along a vertical plate with a low Reynolds turbulence two-equation k-ϵ model, as this topic is fundamental to many applications in science and engineering. Furthermore, the results of the turbulent flow simulations obtained from the current model are corroborated with existing results for the special case of laminar flow and a good correlation is achieved.