The characteristics of wall heat flux (WHF) beneath a supersonic turbulent boundary layer interacting with an impinging shock wave with a 33.2° angle at Mach 2.25 are analyzed using direct numerical simulation. It is found that the QP85 scaling, defined as the ratio of the mean WHF and wall pressure, changes across the interaction. The probability density function of the WHF fluctuations normalized by the local root-mean-squared value is similar to that of wall shear stress. Comparing the WHF and wall pressure spectra shows that the low-frequency shock unsteadiness exhibits little influence on the spectrum. The space–time correlation of the fluctuating WHF reveals that both the streamwise correlation length scale and the convection velocity experience a sharp decrease in the separation region and subsequent recovery in the downstream region. Moreover, the mean WHF in an incident shock interaction is decomposed for the first time. An analysis of the velocity and temperature fluctuations based on bidimensional empirical mode decomposition is performed to evaluate the contribution of turbulent structures with specific spanwise length scales to the mean WHF generation. The decomposed results indicate that the contribution associated with the large-scale structures in the outer region is greatly amplified by the shock interaction and has the leading role in the generation downstream of the interaction.