This work is to extend the elliptic operator, which has been already adopted in turbulent stress model, to fully developed turbulent buoyant channel flows with changing the orientation of the buoyancy vector to be perpendicular to the channel walls. The turbulent heat flux models based on the elliptic concept are employed and closely linked to the elliptic blending second moment closure which is used for the prediction of Reynolds stresses. In order to reflect the stable or unstable stratification conditions, the present model introduces the gradient Richardson number into the thermal to mechanical time scale ratio and model coefficients. The present model has been applied for the computation of stably and unstably stratified turbulent channel flows and the prediction results are directly compared to the DNS data.