We report the results for two buoyancy-driven benchmark cases, heat-driven square cavity at Ra = 10¹¹, and Rayleigh- Bénard convection at Ra = 10⁹, by using the elliptic-blending eddy viscosity k-ε-ζ-f model (or just ζ-f), and three different formulations of the turbulent heat flux, namely the Simple Gradient Diffusion Hypothesis (SGDH), the General Gradient Diffusion Hypothesis (GGDH) and the Algebraic heat Flux Model (AFM). The ζ-f model is well-posed for computing turbulent heat transfer since it contains an approximation of the normal Reynolds stress in the wall-normal direction that is needed in GGDH and AFM formulations. Furthermore, the modelling of the wall-blocking effect by using the elliptic-relaxation approach is physically more sound than the commonly used damping functions. This work is motivated by the recently held 17^th ERCOFTAC SIG15/MONACO2025 workshop on turbulent natural convection flows in differentially heated cavities, [1], which demonstrated superior performance of the ζ-f model in predicting the main flow features for the selected cases.