Accurate representation of heat fluxes is crucial for understanding land–atmosphere interactions and improving atmospheric simulations. However, a common issue arises with flux imbalance, where the measured turbulent heat flux tends to be underestimated due to the nonlocal effects of atmospheric secondary circulations. This study evaluated four flux imbalance prediction models by analyzing data from large eddy simulations performed over heterogeneous land surfaces. For that, a checkerboard pattern of soil moisture was used to define the lower boundary conditions for the atmosphere, across heterogeneity scales ranging from 50 m to 2.4 km. The results show that the selected models can effectively predict flux imbalance when provided with proper semi-empirical factors. The presence of two distinct secondary circulations, thermally-induced mesoscale circulation and turbulent organized structures, account for the nonlinear effect of the heterogeneity scale on the flux imbalance, but it does not affect the performance of the selected models. This study suggests that the flux imbalance prediction models are useful for improving e.g. eddy-covariance measurements. Additionally, a quadrant analysis showed an increasing difference between ejections and sweeps with height, which explains the decrease and increase of the turbulent heat flux and flux imbalance, respectively, and underscores the importance of accounting for vertical variations in turbulent fluxes to represent atmospheric processes accurately.