The flow and thermal field characteristics underlying turbulent convective heat transfer in two reference configurations, representing a jet exiting a planar slot and a circular tube and impinging on a heated wall, are studied computationally. The performance of a conventional Reynolds Stress Model (RSM) and its scale-resolving variant, employed in the context of a time-accurate Reynolds-Averaged Navier–Stokes (RANS) simulation, is investigated in conjunction with different modeling approaches for the turbulent heat flux uj′θ′¯, whose gradient represents the turbulence-induced heat source in the temperature transport equation. In this so-called Sensitized-RANS modeling strategy, employing an improved eddy-resolving instability-sensitive Reynolds Stress Model (IISRSM, Jakirlić and Maduta (2015)), the unresolved ‘sub-scale’ turbulent correlations represent the solutions of a set of corresponding transport equations for the entire stress tensor ui′uj′¯. The turbulent scalar-flux models follow from the classical gradient-diffusion approach with the diffusion-like coefficients formulated as a function of both the scalar eddy-viscosity and the anisotropic Reynolds stress tensor. The latter model expressions also take into account various non-linear relationships in terms of turbulent stress components. The computational results of the fully-developed channel flow, featured only by the mean shear and constant wall heat flux, representing a preliminary investigated flow configuration, agree very well with the corresponding DNS (Direct Numerical Simulation) database for both RSM versions in conjunction with all the turbulent heat flux models adopted. However, the study of the two flow impingement cases with the different structure of the impinging jet and nozzle-to-wall distance variations, characterized by a much more complex flow straining, using the considered conventional Reynolds stress model shows significantly different results for the Reynolds stress and thermal fields compared to the available DNS data, especially in the immediate impact region manifested by a significant overprediction of the turbulence intensity. In contrast, the IISRSM-relevant results for the mean flow and temperature fields, the associated integral characteristics, and the corresponding turbulence quantities are in very good agreement with the reference DNS data, due to the much higher predictive capabilities of this eddy-resolving turbulence model compared to its baseline counterpart.