Building on a method of analysis for the Navier-Stokes equations introduced by Hopf [Math. Ann. 117, 764 (1941)], a variational principle for upper bounds on the largest possible time averaged convective heat flux is derived from the Boussinesq equations of motion. When supplied with appropriate test background fields satisfying a spectral constraint, reminiscent of an energy stability condition, the variational formulation produces rigorous upper bounds on the Nusselt number (Nu) as a function of the Rayleigh number (Ra). For the case of vertical heat convection between parallel plates in the absence of sidewalls, a simplified (but rigorous) formulation of the optimization problem yields the large Rayleigh number bound Nu\ensuremath{\le}0.167 ${\mathrm{Ra}}^{1/2}$-1. Nonlinear Euler-Lagrange equations for the optimal background fields are also derived, which allow us to make contact with the upper bound theory of Howard [J. Fluid Mech. 17, 405 (1963)] for statistically stationary flows. The structure of solutions of the Euler-Lagrange equations are elucidated from the geometry of the variational constraints, which sheds light on Busse's [J. Fluid Mech. 37, 457 (1969)] asymptotic analysis of general solutions to Howard's Euler-Lagrange equations. The results of our analysis are discussed in the context of theory, recent experiments, and direct numerical simulations. \textcopyright{} 1996 The American Physical Society.