The present work presents a stable POD-Galerkin based reduced-order model (ROM) for two-dimensional Rayleigh- Bénard convection in a square geometry for three Rayleigh numbers: 104 (steady state), 3 × 105 (periodic), and 6 × 106 (chaotic). Stability is obtained through a particular (staggered-grid) FOM discretization that leads to a ROM that is pressure-free and has skew-symmetric (energy-conserving) convective terms. This yields long-time stable solutions without requiring stabilizing mechanisms, even outside the training data range. The ROM’s stability is validated for the different test cases by investigating the Nusselt and Reynolds number time series and the mean and variance of the vertical temperature profile. In general, these quantities converge to the FOM when increasing the number of modes convergence, and turn out to be a good measure of accuracy for the non- chaotic cases. However, for the chaotic case, convergence with increasing numbers of modes is not evident, and additional measures are required to represent the effect of the smallest (neglected) scales.