The first principle nonlinear gyrokinetic numerical simulation successfully reproduces the experimental observations of non-diffusive large scale avalanching events in a KSTAR MHD-quiescent L-mode plasma. Power law scaling of electron temperature fluctuation δTe and Hurst exponent factor H from simulation and experiment are in good agreement. In addition, the simulation verifies that the global pattern of mean zonal flow is corrugated with staircase-like structure and responsible for the creased profile of δTe, which is also observed in the experiment. We report on a novel finding that the zonal flow staircase constrains the radial extent of electron heat avalanches through shearing the electron temperature gradient fluctuation, while previous studies were mostly on ion heat transport.