Direct numerical simulations of fully developed turbulent flows in a horizontal square duct heated from below are performed at bulk Reynolds numbers Reb = 3000 and 4400 (based on duct width H) and bulk Richardson numbers 0≤Ri≤1.03. The primary objective of the numerical simulations concerns the characterization of the mean secondary flow that develops in this class of flows. On one hand, it is known that turbulent isothermal flow in a square duct presents secondary mean motions of Prandtl's second kind that finds its origin in the behavior of turbulence structures. On the other hand, thermal convection drives a mean secondary motion of Prandtl's first kind directly induced by buoyancy. As far as the mean structure of the cross-stream motion is concerned, it is found that different types of secondary flow regimes take place when increasing the value of the Richardson number. The mean secondary flow in the range 0.025≲Ri≲0.25 is characterized by a single large-scale thermal convection roll and four turbulence-driven corner vortices of the opposite sense of rotation to the roll, as contrasted with the classical scenario of the eight-vortex secondary flow pattern typical of isothermal turbulent square-duct flow. This remarkable structural difference in the corner regions can be interpreted in terms of combined effects, on instantaneous streamwise vortices, of the large-scale circulation and of the geometrical constraint by the duct corner. When further increasing the Richardson number, i.e., Ri ≳ 0.25, the structure of the mean secondary flow is solely determined by the large-scale circulation induced by the buoyancy force. In this regime, the additional mean cross-stream motion is characterized by the presence of two distinct buoyancy-driven vortices of opposite sense of rotation to the circulation only in two of the four corner regions. With increasing Ri, the large-scale circulation is found to enhance the wall skin friction and heat transfer. In the significant-buoyancy regime Ri ≳ 0.25, the mean cross-stream motion and its rms fluctuations are found to scale, respectively, with the buoyancy-induced velocity ug=√gβΔTH (g, β, and ΔT being the gravity acceleration, the volumetric coefficient of thermal expansion, and the temperature difference across the duct, respectively) and with the mixed velocity scale (√ν/H)ug (ν being the kinematic viscosity). It is suggested that the probable scalings for the rms of streamwise velocity component and of temperature fluctuation are related with the friction velocity uτ and friction temperature Tτ according to the magnitudes uτ2/√(ν/H)ug and Tτuτ/√(ν/H)ug, respectively.