The heterogeneity of sea ice in the marginal ice zone (MIZ) causes multiscale secondary circulations that are challenging to model or observe. Their absent or inadequate representation in ocean-atmosphere interactions in climate models is partially responsible for the underestimation of Arctic sea ice loss. Observationally, such circulations obfuscate the interpretation of polar atmospheric chemistry measurements, among others. To address this open challenge, large-eddy simulations are conducted over real-world satellite-sensed sea ice maps, and idealizations of these maps that alter the ice pattern but conserve its fraction, showing that the ice fraction of a surface is not sufficient to predict surface heat flux. In a second suite of simulations, three other heterogeneity metrics (representing the surface fragmentation, tortuosity, and patch size variability) are then introduced; varying them among simulations suggests that surface-air interactions will be significantly impacted. Based on these results, we then show, using a multi-linear regression of map features, that MIZ surface parameterizations in Earth Systems Models can be improved if they account for these various features of the surface. Finally, the simulations also illustrate that the contribution of roughness changes to the resulting atmospheric circulations are minimal compared to thermal heterogeneity in the MIZ.