  We carry out direct numerical simulation (DNS) of flow in a turbulent square
duct by focusing on heat transfer effects, considering the case of unit Prandtl
number. Reynolds numbers up to $Re_\tau \approx 2000$ are considered which are
much higher than in previous studies, and which yield clear scale separation
between inner- and outer-layer dynamics. Close similarity between the behavior
of the temperature and the streamwise velocity fields is confirmed as in
previous studies related to plane channels and pipes. Just like the mean
velocity, the mean temperature is found to exhibit logarithmic layers as a
function of the nearest wall, however with a different slope. The most
important practical implication is the validity of the traditional hydraulic
diameter as the correct reference length for reporting heat transfer data, as
we rigorously show here. Temperature and velocity fluctuations also have
similar behavior, but apparently logarithmic growth of their inner-scaled peak
variances is not observed here unlike in canonical wall-bounded flows. Analysis
of the split contributions to the heat transfer coefficient shows that mean
cross-stream convection associated with secondary motions is responsible for
about $5\%$ of the total. Finally, we use the DNS database to highlight
shortcomings of traditional linear closures for the turbulent heat flux, and
show that substantial modeling improvement may be in principle obtained by
retaining at least the three terms in the vector polynomial integrity basis
expansion.
