Turbulent fluxes make a substantial and growing contribution to the energy balance of ice surfaces globally, but are poorly constrained owing to challenges in estimating the aerodynamic roughness length (z0). Here, we used structure from motion (SfM) photogrammetry and terrestrial laser scanning (TLS) surveys to make plot-scale 2-D and 3-D microtopographic estimations of z0 and upscale these to map z0 across an ablating mountain glacier. At plot scales, we found spatial variability in z0 estimates of over two orders of magnitude with unpredictable z0 trajectories, even when classified into ice surface types. TLS-derived surface roughness exhibited strong relationships with plot-scale SfM z0 estimates. At the glacier scale, a consistent increase in z0 of ~0.1 mm d−1 was observed. Space-for-time substitution based on time since surface ice was exposed by snow melt confirmed this gradual increase in z0 over 60 d. These measurements permit us to propose a scale-dependent temporal z0 evolution model where unpredictable variability at the plot scale gives way to more predictable changes of z0 at the glacier scale. This model provides a critical step towards deriving spatially and temporally distributed representations of z0 that are currently lacking in the parameterisation of distributed glacier surface energy balance models.