Turbulent secondary flows are defined as flows where the mean vorticity vector is aligned with the mean wind direction. They are typically classified as Prandtl’s secondary flow of the first or second kind and can be produced by stretching and/or tilting of vorticity (first kind) or spatial heterogeneity in the Reynolds stresses (second kind). Because these secondary flows produce coherent mean upwelling and downwelling, they are of great importance for understanding and predicting land-atmosphere exchanges of momentum, heat, water vapor, and trace gases. Previous studies have examined turbulent secondary flows produced by spanwise heterogeneity in the aerodynamic roughness length, which creates spatial variability in the Reynolds stresses, resulting in regions of mean upwelling and downwelling of regions of high and low roughness, respectively. However, the extent to which thermal heterogeneities can create and sustain turbulent secondary flows remains an open question. Using a suite of large eddy simulations of moderately convective boundary layers with uniform aerodynamic roughness but spanwise-variable surface heat flux, we investigate the mechanisms responsible for producing and sustaining thermally-driven turbulent secondary flows. We find that shear and buoyancy production over elevated heat flux regions necessitates lateral entrainment of fluid with low turbulent kinetic energy (TKE), inducing mean counter-rotating cells aligned such that upwelling and downwelling occur over high and low heat flux regions, respectively. This result illustrates that buoyancy production of TKE alters aggregate flow response and thus is a distinctly different mechanism responsible for sustenance of secondary flows than others that have been identified previously.