The statistical behaviour and modelling of turbulent fluxes of the reaction progress variable and non-dimensional temperature in the context of Reynolds-Averaged Navier–Stokes (RANS) simulations have been analysed for flame–wall interactions within turbulent boundary layers. Three-dimensional Direct Numerical Simulation (DNS) databases of two different flame–wall interaction configurations—(i) statistically stationary oblique wall quenching (OWQ) of a V-flame in a turbulent channel flow and (ii) unsteady head-on quenching (HOQ) of a statistically planar flame propagating across a turbulent boundary layer—have been considered for this analysis. Scalar fluxes of both the temperature and reaction progress variable exhibit counter-gradient behaviour at all times during unsteady HOQ of statistically planar turbulent premixed flames considered here. In the case of statistically stationary V-flame OWQ, the scalar fluxes of both reaction progress variable and temperature exhibit counter-gradient behaviour before quenching, but gradient behaviour has been observed close to the wall once the flame begins to quench. The weakening of the effects of thermal expansion close to the wall as a result of flame quenching gives rise to a gradient type of transport for the streamwise component in the oblique quenching of the V-flame. It has been found that the relative orientation of the flame normal vector with respect to the wall normal vector needs to be accounted for in the algebraic scalar flux closure, which can be applied to different flame/flow configurations. An existing algebraic scalar flux model has been modified in this analysis for flame–wall interaction within turbulent boundary layers, and it has been demonstrated to capture the turbulent fluxes of the reaction progress variable and non-dimensional temperature reasonably accurately for both configurations considered here based on a priori DNS analysis.