  This paper presents an algorithm for calculating an ensemble of solutions to
natural convection problems. The ensemble average is the most likely
temperature distribution and its variance gives an estimate of prediction
reliability. Solutions are calculated by solving two coupled linear systems,
each involving a shared coefficient matrix, for multiple right-hand sides at
each timestep. Storage requirements and computational costs to solve the system
are thereby reduced. Moreover, this paper addresses a need for higher order
methods to solve natural convection problems. Stability and convergence of the
method are proven under a timestep condition involving fluctuations of the
velocity. Numerical tests are provided which confirm the theoretical analyses.
