The early impact bombardment extensively fractured the lunar crust resulting in the formation of the so-called megaregolith. Previous estimates of megaregolith distribution vary significantly with respect to the vertical extent and the size-frequency distribution of fragments was rarely studied. We built a spatially resolved numerical model to simulate the process of cumulative impact fragmentation, aiming to backtrack the megaregolith evolution history and to constrain its fragment distribution. The results highlight the pivotal role of basin-forming events on the megaregolith formation. Especially the South-Pole Aitken (SPA) impact established the initial megaregolith structure which remained distinct after 0.5 Ga subsequent fragmentation. At 3.8 Ga, the megaregolith displays substantial lateral variation and layering: the highly fractured upper layer of ∼2.5 km is dominated by meter-scale fragments; the disturbed lower layer deeper than tens of kilometers is mainly consisting of kilometer-scale fragments; the transition zone >5 km contains fragments of various size scales.