Direct numerical simulations have been performed of natural convection in a closed cavity with constant heat fluxes on the vertical sides. Results are presented for Rayleigh numbers in the range 3×1011–4×1014 and Prandtl number equal to 1. The flow is characterized by two boundary layers near the heated and cooled vertical sides and two intrusion layers near the horizontal walls, which complete a recirculation pattern. The net energy surplus transferred upwards by both boundary layers is balanced by heat conduction through the main body of the cavity. The vertical boundary layer thickness is of O(Ra−2/9) and the intrusion layers of O(Ra−1/6). Scalings for thicknesses, velocities and temperature gradients, obtained from order of magnitude estimates, are confirmed by the numerical simulations. The bulk temperature gradient is much larger than those across the thin boundary layers and, therefore, effective in delaying the onset of instability. Unstable flows are thus found only above a Rayleigh number of the order of 1012, in contrast to a cavity with fixed lateral wall temperatures where this is close to 108. The threshold for instability and the dominant frequencies are in satisfactory accord with relevant linear stability results from the literature. The bulk of the cavity responds passively and filters the excitation from the unstable boundary layers by retaining frequencies close to the Brunt-Vaisala frequencies in a density stratified fluid.