A numerically exact approach is proposed to study heat transport through a molecular junction that is coupled to two phonon reservoirs at different temperatures. Within a spin-boson type model, the junction is represented by a two level subsystem and the phonon reservoirs are described by two independent harmonic baths. Quantum simulations of the thermal flux (heat current) are performed employing the multilayer multiconfiguration time-dependent Hartree theory. Results for the transient and stationary heat current are discussed for different physical regimes. Furthermore, the validity of Red field theory, as an example of an approximate approach applicable for weak system-bath coupling, is examined using the proposed numerically exact approach. (C) 2008 Elsevier B.V. All rights reserved.