The convective flow in a thin liquid layer with a free surface heated from below is studied using a combination of accurate experiments with silicone oil (v=0.1 cm2 s−1) and high-resolution direct numerical simulations of the time-dependent governing equations. It is demonstrated that above a certain value εs of the threshold of primary instability, ε=0, square convection cells rather than the seemingly all-embracing hexagons are the persistent dominant features of Bénard convection. The transition from hexagonal to square cells sets in via a subcritical bifurcation and is accompanied by a sudden rapid increase of the Nusselt number. This implies that square cells are the more efficient mode of heat transport. Their wavenumber exceeds that of hexagonal cells by about 8%. The transition depends on the Prandtl number and it is shifted towards higher εs if the Prandtl number is increased. The replacement of hexagonal by square cells is mediated by pentagonal cells. In the transitional regime from hexagonal to square cells, characterized by the presence of all three planforms, the system exhibits complex irregular dynamics on large spatial and temporal scales. The time dependence becomes more vivid with decreasing Prandtl number until finally non-stationary square cells appear. The simulations agree with the experimental observations in the phenomenology of the transition, and in the prediction of both the higher Nusselt number of square Bénard cells and the subcritical nature of the transition. Quantitative differences occur with respect to the values of εs and the Prandtl number beyond which the time dependence vanishes. These differences are the result of a considerably weaker mean flow in the simulation and of residual inhomogeneities in the lateral boundary conditions of the experiment which are below the threshold of control.