This article presents the local stability analysis of the boundary layer flows over a heated surface with variable viscosity using a non-modal approach. The three-different variable viscosity models have been considered (namely, linear, exponential, and inverse models) and they are a function of temperature only. The base flow profiles deflect towards the wall with increased sensitivity parameter \((\epsilon > 0)\) resulting in narrow boundary layers for all viscosity models. However, the reverse effect has been observed in all models for \(\epsilon < 0\) . The standard procedure has been adopted to obtain the linearized governing stability equations in the normal velocity-vorticity form by eliminating the pressure term. The heat balance equation is coupled with the normal velocity-vorticity equations by temperature-dependent viscosity. The governing stability equations are spatially discretized using the Chebyshev spectral collocation approach. The discretized governing equations along with appropriate boundary conditions form an initial value problem (IVP). The transient growth in energy is mathematically computed by linear superposition of the non-orthogonal eigenvectors and large energy amplification has been achieved for a short time even though all eigenmodes lie down in the lower stable half-plane. The maximum transient energy growth \((G_{\max })\) reduces with the increased sensitivity parameter \((\epsilon > 0)\) for all viscosity models resulting in an increased critical Reynolds number and flow becoming modally more stable. However, the reverse effect has been found for \(\epsilon < 0\) . Among all viscosity models, the linear viscosity model produced a lower energy amplification and temporal growth rate \((\omega _i)\) for a given value of \(\epsilon \) .