A novel data-driven approach is used to describe a spatially varying turbulent diffusivity coefficient for the Higher Order Generalised Gradient Diffusion Hypothesis (HOGGDH) closure of the turbulent heat flux to improve upon RANS cooling predictions in film cooling flows. Machine learning algorithms are trained on two film cooling flows and tested on a case of a different density and blowing ratio. The Random Forests and Neural Network algorithms successfully reproduced the LES described coefficient and the magnitude of the turbulent heat flux vector. The Random Forests model was implemented in a steady RANS solver with a k-ω SST turbulence model and applied to four cases. All cases saw improvements in the predicted Adiabatic Cooling Effectiveness (ACE) over the cooled surface compared to the standard Gradient Diffusion Hypothesis (GDH) approach, but only minor improvements in the centreline and lateral spread are seen compared to a HOGGDH model with a constant cθ of 0.6. Further improvements to cooling predictions are highlighted by extending these data-driven approaches into turbulence modelling to improve flow field predictions.