This paper presents the results of large-eddy simulation (LES) of a downward wall jet that encounters an upward flow in a vertical rectangular channel, named negative buoyant fountain (NBF). The downward jet fluid has higher temperature than the opposing upward stream so that the influence of the buoyancy force on the physics of the mixing process can be investigated. Two subgrid-scale models available in the OpenFOAM platform, i.e. the dynamic k-equation model (DKM) and the wall-adapting local eddy-viscosity (WALE) model, are used for cross-comparison. The numerical results of the time-averaged and fluctuating velocity fields agree satisfactorily with the experimental data. The LES results form an accurate, high-resolution, three-dimensional database of buoyant turbulent mixing, which has expanded the sparse flow field measurements provided by the experiment, hence the numerical database could be used to validate the Reynolds-averaged Navier-Stokes (RANS) simulation method and to contribute to the turbulence modeling. In particular, the statistics of turbulent heat flux have been analyzed and various turbulent heat flux models have been a priori tested based on the LES numerical database. Such information could be useful for RANS simulation guidelines or for deriving anisotropic turbulent heat flux models applicable for turbulent buoyant flows sharing the similar principal flow structures in complex geometries.