For geophysical and environmental flows in the high-Reynolds-number regime, stable density stratification strongly affects the turbulent fluid motions. While turbulent flows, due to the presence of a cascade of spatial and temporal scales, present several challenges to accurate numerical approximation, density stratification in these flows exacerbates the situation further by suppressing vertical velocity fluctuations thereby enhancing the anisotropy of the turbulence. In this paper, based on the framework of residual based variational multiscale (RBVMS) methods, we design a new numerical formulation for incompressible stratified flows that introduces coupling between the velocity fine scales and densityequation residual, and gives improved numerical performance as evidenced by the results of two challenging turbulent flow simulations. (C) 2016 Elsevier Ltd. All rights reserved.