Surface-tension-driven Bénard convection in low-Prandtl-number fluids is studied by means of direct numerical simulation. The flow is computed in a three-dimensional rectangular domain with periodic boundary conditions in both horizontal directions and either a free-slip or no-slip bottom wall using a pseudospectral Fourier–Chebyshev discretization. Deformations of the free surface are neglected. The smallest possible domain compatible with the hexagonal flow structure at the linear stability threshold is selected. As the Marangoni number is increased from the critical value for instability of the quiescent state to approximately twice this value, the initially stationary hexagonal convection pattern becomes quickly time-dependent and eventually reaches a state of spatio-temporal chaos. No qualitative difference is observed between the zero-Prandtl-number limit and a finite Prandtl number corresponding to liquid sodium. This indicates that the zero-Prandtl-number limit provides a reasonable approximation for the prediction of low-Prandtl-number convection. For a free-slip bottom wall, the flow always remains three-dimensional. For the no-slip wall, two-dimensional solutions are observed in some interval of Marangoni numbers. Beyond the Marangoni number for onset of inertial convection in two-dimensional simulations, the convective flow becomes strongly intermittent because of the interplay of the flywheel effect and three-dimensional instabilities of the two-dimensional rolls. The velocity field in this intermittent regime is characterized by the occurrence of very small vortices at the free surface which form as a result of vortex stretching processes. Similar structures were found with the free-slip bottom at slightly smaller Marangoni number. These observations demonstrate that a high numerical resolution is necessary even at moderate Marangoni numbers in order to properly capture the small-scale dynamics of Marangoni convection at low Prandtl numbers.