Large-eddy simulation, with recently developed dynamic subgrid-scale models, is used to study the effect of heterogeneous surface temperature distributions on regional-scale turbulent fluxes in the stable boundary layer (SBL). Simulations are performed of a continuously turbulent SBL with surface heterogeneity added in the form of streamwise transitions in surface temperature. Temperature differences between patches of 6 and 3 K are explored with patch length scales ranging from one-half to twice the equivalent homogeneous boundary layer height. The surface temperature heterogeneity has important effects on the mean wind speed and potential temperature profiles as well as on the surface heat flux distribution. Increasing the difference between the patch temperatures results in decreased magnitude of the average surface heat flux, with a corresponding increase in the mean potential temperature in the boundary layer. The simulation results are also used to test existing models for average surface fluxes over heterogeneous terrain. The tested models fail to fully represent the average turbulent heat flux, with models that break the domain into homogeneous subareas grossly underestimating the heat flux magnitude over patches with relatively colder surface temperatures. Motivated by these results, a new parameterization based on local similarity theory is proposed. The new formulation is found to correct the bias over the cold patches, resulting in improved average surface heat flux calculations.