Algebraic closures for the turbulent scalar fluxes were evaluated for a discrete hole film cooling geometry using the results from a high-fidelity large eddy simulation (LES). Several models for the turbulent scalar fluxes exist, including the widely used gradient diffusion hypothesis (GDH), the generalized GDH (GGDH), and the higher-order GDH (HOGGDH). By analyzing the results from the LES, it was possible to isolate the error due to these turbulent mixing models. Distributions of the turbulent diffusivity, turbulent viscosity, and turbulent Prandtl number were extracted from the LES results. It was shown that the turbulent Prandtl number varies significantly spatially, undermining the applicability of the Reynolds analogy for this flow. The LES velocity field and Reynolds stresses were fed into a Reynolds-averaged Navier–Stokes (RANS) solver to calculate the fluid temperature distribution. This analysis revealed in which regions of the flow various modeling assumptions were invalid and what effect those assumptions had on the predicted temperature distribution.