I made a one line change to entropy_viscosity.cc to disable stabilization near the boundary for SUPG in addition to EV and the heat fluxes are closer to what I would expect for the first couple of timesteps. This suggests that yes, the issue was related to the inclusion of the SUPG parameter when calculating heat fluxes. I am testing this branch now to see whether this also solves the failure to converge issue but it will take a day or so to run.
I think that if the non-physical contribution to thermal conductivity is going to be included in the heat transport, it should only be applied in the direction of the streamlines, so effectively (kI + supg_parameterv/norm(v)) * grad(T) where k is the thermal conductivity, I is the identity matrix, supg_parameter is the stabilization parameter and v/norm(v) is a unit vector along a streamline.