A non-linear subgrid-scale (SGS) heat flux model is introduced in large eddy simulation for turbulent thermal flows. Unlike the linear isotropic eddy diffusivity model, the proposed model accounts for the SGS heat flux in terms of the large-scale strain-rate tensor and the temperature gradients. This is equivalent to using a tensor diffusivity. The model is to some extent similar to a scale-similarity model subjected to a Taylor expansion for the filtering operation. The formulation leading to the present proposal is discussed. The model is examined in LES for a buoyant flow in an infinite vertical channel with two differentially heated side walls. It is shown that the proposed model reproduces reasonable results as compared with the isotropic SGS diffusivity model and DNS data.