In this work, the effects of thermal buoyancy and rotation angle on the hydrodynamics and mixed convection heat transfer characteristics over four heated identical square cylinders placed in a parallel horizontal channel are numerically studied. For the flow field, the Lattice Boltzmann method (LBM) with multiple relaxation time collision operator (MRT) is implemented and the Bhatnagar-Gross-Ktook collision operator (BGK) is adopted for the thermal field simulation. The distance between two adjacent square cylinders is fixed for a parametric study and a uniform velocity is imposed at the inlet region. Ranges of the influence parameters are defined as: a ¼ 45 , 90 , 120 ; 50 Re 150; 0 Ri 0.5. Several flow and heat transfer parameters including the time-averaged lift/drag coefficients, the global time-averaged Nusselt number and the distributions of local Nusselt number are analyzed, respectively. The vorticity and isotherm patterns are also plotted to illustrate the visual disturbance of the thermal buoyancy and the rotation angle on the flow and thermal fields. Numerical results show that the flow and thermal patterns for all arrangements become disordered and nonsystematic with the introduction of thermal buoyancy, and the thermal pattern for the downstream cylinders has an upward drift in the wake region. The time-averaged lift coefficient for each square cylinder is relatively sensitive to the influence parameters and the arrangement, and it is always less than zero when the Richardson number is not equal to zero. The global timeaveraged Nusselt number is directly proportional to the Reynolds number and Richardson number, and the local Nusselt numbers on the front and rear faces for each square cylinder reach a maximum and minimum, respectively. It is also observed that the distribution of local Nusselt number for all rotation angles is intensively affected by the thermal buoyancy.