Increased small impactor flux on the Moon as inferred from regolith thickness at the Chang’E-5 Region

Domain: earth_science
Published: 2025-12-22
Editorial outcome: major_revision

Abstract
Knowledge of the small impactor flux on the Moon can provide critical insights into lunar geological evolution and the dynamics of the inner solar system. However, it cannot be directly deciphered from the observed populations of small craters due to their rapid degradation. Impact-generated lunar regolith records the contribution from historical impacts of various sizes and is particularly sensitive to sub-hectometer craters, providing a unique means to estimate the small impactor flux. Here we estimate the distribution of regolith thickness in the Chang’E-5 landing region (the youngest radiometric-dated homogeneous mare unit) through both statistics of fresh crater morphology and regolith evolution simulations using different impact flux hypotheses. The large median regolith thickness of 3.6 ± 0.3 m derived from crater morphology can be attributed solely to a higher production rate of small lunar craters at sub-hectometer scales, which is 1.2−1.6 times higher than previously estimated. Such an increased small impactor flux may occur at ∼600 Ma and indicate recent asteroid breakup events, or could have persisted over the past 2 Ga, as a result of the preferential delivery of small impactors through the Yarkovsky and YORP effects.

Introduction
Knowledge of the impactor flux onto the Moon is significant for dating the lunar surface and understanding its geological evolution1,2. It can also provide critical insights into the populations of main belt asteroids and the dynamics of the inner solar system3,4. Impact craters preserve key records of the bombardment history5, serving as the most direct means of estimating the impactor flux. By linking the radiometric ages of the Apollo, Luna and Chang’E (CE) samples with the crater size-frequency distributions (CSFDs)6, the lunar chronology function was established1,7, which supports a relatively constant cratering rate over the past 3 Ga. However, due to multiple erosion processes (for example, micrometeorite bombardment, ejecta burial, and seismic shaking)8–10, impact craters undergo a gradual increase in diameter and decrease in depth since their formation11,12, only surviving for a limited time on the lunar surface13. Particularly, smaller impact craters degrade faster than larger ones11, with craters below 20–50 m in diameter surviving for <∼50−125 Myr13. Craters smaller than a few hundred meters across have typically reached an equilibrium abundance in the lunar mare14. The scale-dependent nature of crater degradation can lead to a considerable bias between the observed size-frequency distribution of small craters and their true production11. Therefore, previous estimates of the small impactor flux based on crater populations maybe inaccurate, as supported by multiple lines of evidence15–17. For example, recent observations from the Lunar Reconnaissance Orbiter Camera (LROC) Narrow Angle Camera (NAC) temporal image pairs revealed that the lunar crater production rate at decameter scales is 1.33 times higher than the widely accepted Neukum production function (NPF)15. The 40Ar/39Ar dating of impact glass beads in the Apollo 12, 14, 16, and 17 samples showed that their formation ages are primarily concentrated into the most recent ∼400 Ma16,17, which may indicate a higher small impactor flux.

Impact products offer another unique perspective for estimating the impactor flux18,19. Almost the entire lunar surface is covered with a fragmental regolith layer with thickness ranging from a few meters to tens of meters20, which is primarily generated by the impact cratering process21. Shortly after the emplacement of basaltic mare lava flows, little to no regolith is present22. Meteorite bombardment by both small and large impactors can shatter the coherent protolith and produce new regolith22. The cratering rate on the lunar surface directly controls the regolith growth rate, referred to as the impact flux trend23. As regolith thickens, small impactors overturn the existing regolith layer24, while larger impactors capable of penetrating the pre-existing regolith can fracture the underlying substrate25. As a result, regolith growth rate decreases over time and becomes lower in older units, a phenomenon known as the buffering trend or self-limiting process22,23. Considering the formation mechanism, lunar regolith is generally expected to be thicker in older surfaces, which has long been presumed in previous studies18,22,26–30. Observations from radar31, microwave radiometer32, and crater morphology20 collectively suggest that regolith thickness across lunar maria is statistically correlated with age (Supplementary Text 1 and Supplementary Figs. 1–5), though the correlation is not particularly strong due to several factors (measurement uncertainty, target property, randomness of large impacts, and the buffering effect; Supplementary Text 2). Besides, the growth of lunar regolith is dominated by and sensitive to craters at sub-hectometer scales25, making regolith thickness particularly suitable for deciphering the small impactor flux. Recent advances in crater shape models, cratering mechanisms, and topographic degradation have facilitated the development of a 3-D regolith evolution model25,33 that can well characterize the self-limiting growth process of lunar regolith resulting from continuous bombardment. By comparing the observed and modeled regolith thicknesses, new constraints on the small lunar impactor flux can be provided.

To achieve this, a geologically young, isotopically dated lunar landing area with homogeneous target properties representative of typical mare basalts is required. First, high-resolution optical imagery of the landing area enables a reliable estimate of regolith thickness. The precise radiometric age avoids the large dating uncertainty associated with crater statistics. Second, the impact flux inferred from a representative mare unit can serve as a benchmark for comparison with the standard production function (for example, the NPF1), which was established across broader mare regions. Third, the lunar regolith developed in such a unit is formed mainly by small, local impacts with less influence from exotic materials or other geological processes. Most large source craters that contributed to regolith formation are only slightly degraded and thus remain clearly identifiable. Finally, the growth of a relatively thin regolith layer is less affected by the buffering effect, enhancing the ability to detect variations in impact flux from regolith thickness. To date, a number of lunar samples with isotopic ages ranging broadly from >3.9 Ga to 25 Ma (Cone crater) have been acquired6. However, most of them are older than 3 Ga, and the corresponding geological units have undergone complex evolution processes, making it difficult to reconstruct the growth history of regional regolith. Some samples collected by the Apollo missions are inferred to source from the Copernican-aged young impact craters2,6 (for example, Copernicus and Tycho). However, their impact melt sheets show non-uniform geological features (for example, cracks and knobs; Supplementary Fig. 6), which could affect the regolith growth process.

China’s CE-5 landing site (43.058°N, 51.916°W) is located on the youngest isotopically dated mare basalt unit (P58/EM4) in the northern Oceanus Procellarum, with a precise radiometric age of 2.03 ± 0.004 Ga34. Spectroscopic observations show that the landing region exhibits uniform target characteristics several tens of kilometers across35 (Supplementary Fig. 7). Statistical analysis of rocky crater populations across 15 different mare units in Oceanus Procellarum, Mare Imbrium, Mare Nubium, Mare Tranquillitatis, Mare Humorum, and Mare Serenitatis reveals that the P58 unit exhibits no anomalous target mechanical properties36. Also, the rock abundance of the P58 unit aligns well with the temporal trend observed across 291 mare units on the Moon, indicating that the P58 unit is representative of typical mare basalts37 (Methods and Supplementary Fig. 8). Besides, high-resolution descent camera data covering the landing area, as well as the thorough investigation of source craters for the CE-5 samples38, facilitate an accurate identification of crater morphology and systematic modeling of regolith evolution. Therefore, the CE-5 landing region presents a unique opportunity to unravel the small impactor flux during the past 2 Ga.

In this study, we provide evidence for an increased flux of lunar impactors less than tens of meters, by comparing the regolith thickness in the CE-5 region estimated using fresh crater morphology with that from regolith growth modeling.

Results
Regolith thickness from crater morphology
Photographs from Ranger and Lunar Orbiter (LO) showed that the small fresh craters on the lunar surface exhibit four typical morphological types39: normal, central mound, flat-bottomed, and concentric. Extensive impact experiments revealed that the formation of these craters with different morphologies depends on the ratio of crater size to local regolith thickness40 (see Methods). Compared with in situ seismic experiments41,42 and radar remote sensing techniques43, the morphology and size-frequency distribution of small fresh craters can be used to derive the distribution of regolith thickness over a large area, independent of the geophysical model. Therefore, this approach has been widely used for regolith thickness estimation in different regions across the Moon12,18,20,27.

High-resolution imagery with appropriate solar incidence angle is key to ensuring the reliability of this method, as it determines the accuracy of crater morphology identification. Existing studies35,44,45 on the regolith thickness in the CE-5 region were conducted based on optical images with a low resolution of 1.5 m, which could lead to misidentifications of crater morphology (see Supplementary Fig. 9 for an example) or biased crater samples toward the large diameter ranges. If a central mound, flat-bottomed, or concentric crater is misidentified as a bowl-shaped normal one, then the area fraction with thick regolith deposits increases and thus the regional regolith thickness would be overestimated (see Methods for detail). Their results indicated a regolith thickness of 5 m35 or 4–6 m44,45 in the CE-5 region, which is generally thicker than that in the older Apollo landing sites (3.7−4.4 m41,42 from seismic experiments and 3.1–5.0 m46 from crater morphology; Supplementary Table 1). Considering the statistical correlation between regolith thickness and surface age, this may again suggest an overestimated thickness.

Here, we select a flat region across the CE-5 landing site with homogeneous spectral characteristics for crater morphology identification (17.3 km2; Supplementary Fig. 10), where the descent camera images are available. An obvious secondary crater field is excluded because secondary craters tend to be normal ones and may affect the statistics. The FeO/TiO2 content, optical maturity, and crater density of the study area all fall within one standard deviation of the mean value of those in the P58/EM4 unit (34,797 km2; Supplementary Fig. 11), indicating that the selected CE-5 landing region can well represent the whole basalt unit in terms of target properties and impact history. The CE-5 descent camera images (<0.5 m/pixel), the generated digital orthophotomap (DOM) and digital terrain model (DTM) (0.5 m/pixel), as well as four LROC NAC images (0.8−0.9 m/pixel; Supplementary Table 2) are used for crater measurements (see Methods). To ensure the accurate identification of crater diameter and morphology, only craters with diameters larger than 10 pixels were considered47. In total, 361 fresh impact craters with diameter ranging from 5 to 120 m were counted (Fig. 1a), including 174 normal, 118 central mound and flat-bottomed, and 69 concentric ones. Based on the fractions of these counted craters using high-resolution imagery (Supplementary Fig. 12), the regional regolith thickness is estimated to vary from ∼1 to 8 m, with a median of 3.6 ± 0.3 m (1σ uncertainty; Fig. 1b). This is in line with the observation from the CE-5 Lunar Regolith Penetrating Radar48 (>2.0−2.5 m), and generally lower than the regolith thickness in the older Apollo landing regions (Supplementary Table 1). Using a radiometric age of 2.03 ± 0.004 Ga, the mean regolith growth rate in the CE-5 region is calculated to be 1.8 ± 0.15 m/Gyr, which is larger than those in the Apollo landing regions (for example, 1–1.4 m/Gyr in the Apollo 11, 12, 15, and 17 regions; Supplementary Table 1). This could result from the self-limiting formation nature of lunar regolith22, which leads to a progressively decreasing growth rate as regolith thickens over time, or reflect an elevated small impactor flux recently. Systematic simulations on regolith evolution in the CE-5 region will provide further insights into the underlying mechanisms.

Regolith growth modeling
The production, distribution, and transportation of lunar regolith are closely linked to the impact cratering process21,22. We have developed a Lunar Topography and Regolith Evolution Model (LTREM) to simulate the formation and growth of lunar regolith resulting from continuous bombardment and topographic degradation25,33. When each individual impact crater forms, we consider the fragmentation of bedrock during the excavation stage, the formation of a breccia lens during the modification stage, and the secondary cratering during the ejecta deposition process49 (Methods). LTREM can quantitatively simulate regolith growth process by accounting for various components, including contributions from impact craters of various sizes and origins (i.e., primary and secondary), crater morphologies, the buffering effect of the existing regolith layer, and topographic degradation. This enables a systematic evaluation of the key factors that govern the present-day regolith thickness observed on the lunar surface. The uncertainty in the simulated regolith thickness arises from three main sources: the populations of primary and secondary craters, the production of regolith by individual craters, and the topographic degradation rate (see Methods).

In this study, we utilize an observation-constrained LTREM with consideration of global-scale impact events to investigate the regolith evolution over a 2×2 km area around the CE-5 landing site from 2.03 Ga to the present. The minimum crater diameter is set to be 2 m (twice the pixel resolution). Smaller craters are generally unable to penetrate the pre-impact regolith layer to produce new regolith (Supplementary Fig. 13), thus hardly affecting the results. Specifically, the input crater populations include three parts: the global source craters of the CE-5 regolith, invisible eroded craters at the local scale, and secondary craters formed by primaries (Methods). Based on our previous investigations38,50, there are eight major source craters outside the P58/EM4 unit (Supplementary Table 3 and Fig. 2a), and 2171 craters inside P58/EM4 (Fig. 2b). In particular, all the visible craters within a 2×2 km study area across the CE-5 landing site were counted (Fig. 2c). There are 1223 in total with diameters ranging from ∼4 to 408 m. Their size-frequency distribution aligns with the NPF at hectometer scales but falls below the NPF at diameters smaller than 85 m (Fig. 2d). This indicates a large number of small, invisible craters that once existed on the lunar surface but have been erased by erosion processes10. The populations of these invisible craters can be accounted for by using a crater production function (for example, the NPF). As will be shown later, they contribute significantly to lunar regolith growth and therefore their actual production can be well constrained through regolith thickness. The formation ages of the most important contributors among the CE-5 source craters are estimated based on the degradation states or the size-frequency distributions of superposed small craters (Supplementary Table 3; Methods), which are incorporated in our simulation. To ensure statistical significance, all the subsequent modeling results for each parameter set are averaged over at least 100 realizations, as further realizations yield negligible changes in the outcomes (Methods).

We first present the results by taking the NPF as the crater production function. As shown in Fig. 3, while new craters continue to form in the simulation area, older craters degrade progressively, and some fade away eventually. Large primary craters such as Philolaus (72.22°N, 32.88°W; 71.4 km) outside the area also contributed to the development of local regolith through their distal ejecta and secondary craters. The regolith layer thickens over time and is continuously reworked, as indicated by infilling within crater interiors and the redistribution of ejecta deposits. The populations of primary and secondary craters produced in the simulation area over 2 billion years are shown in Supplementary Fig. 14. Since the formation of regolith thickness is most sensitive to small craters ranging from several meters to 10s of meters (Supplementary Fig. 13)25, the presence of over 4×106 craters with D>2 m ensures the statistical significance required to infer the impactor flux. As shown in Fig. 4a, the spatial distribution of regolith thickness is highly correlated with impact craters, ranging from nearly zero on the steep crater walls to more than ten meters in the crater interior (breccia lenses) and exterior (ejecta blankets) regions. The median regolith thickness across the simulation area is 2.90 ± 0.14 m when both the breccia lenses and ejecta are considered as regolith (1σ uncertainty; see Methods), and 1.27 ± 0.06 m if only the ejecta blankets are incorporated. The ejecta-only case represents the thickness of the fine-grained part within the regolith layer, which is significantly lower than that derived from crater morphology (3.6 ± 0.3 m). Therefore, we follow previous research22,29, considering both breccia lenses and ejecta as lunar regolith in the subsequent analysis, which can be regarded as, to some extent, an upper-bound estimate.

From a perspective of source craters, we find that 90.9 ± 0.1% of the regolith is generated by primary craters, with 50.2% from the observed craters and 40.7% from eroded primary craters. The remaining 9.1% of lunar regolith is produced by secondary craters. The simulations also show that the majority of the regolith at the CE-5 landing site originated within a 1 km radius (∼93.5%), with only 5.2 ± 1.8% of material being distal ejecta from large impact craters outside the P58/EM4 unit (Fig. 4b). For the first time, we incorporate the contribution from eroded craters to the regolith thickness and provide a more accurate estimate for the fraction of exotic materials in the CE-5 samples, which is consistent with laboratory analyses (∼5% exotic materials)51,52. In terms of the growth dynamics, the median and quartiles of regolith thickness grow smoothly with time in general, with several leaps associated with the formation of Xu Guangqi and three large impact craters outside the P58/EM4 unit (labeled in Fig. 4c). The invisible craters play a dominant role in regolith accumulation during the initial stage (the red line in Fig. 4d) while their contribution decreases gradually as regolith thickens, reflecting the regolith buffering trend. The contributions of observed craters show a similar, though less pronounced, decreasing trend at first, but begin to increase at about 600 Ma (the black line in Fig. 4d). This may imply a higher recent impact flux, but further evidence is required because the number of craters with age information (listed in Supplementary Table 3) is insufficient to ensure statistical significance. The overall effect of secondary craters on the regolith growth rate is quite weak (typically accounting for 1%; the blue line in Fig. 4d), except during periods when large impact craters formed and deposited ejecta into the landing region (for example, Philolaus, Harpalus, Sharp B, and Copernicus). This is primarily due to the substantially shallower excavation depths of secondary craters compared with primaries53. As a result, secondary craters play a dominant role in mixing the surficial regolith24 but contribute little to the formation of new regolith.

Evidence for an increased impactor flux
Through a large number of simulations with varying parameter sets, including the transient crater depth, secondary crater population, and topographic diffusivity (see Methods), the model prediction (2.90 ± 0.14 m) when using the NPF as the production function is always found to be lower than the observation (3.6 ± 0.3 m). Such a difference is statistically significant at the 2σ level (>95% confidence; Methods), indicating that it cannot be attributed to the variation of other factors but likely results from an underestimated production of small, primary craters. We further test three different crater production function hypotheses other than the NPF (PF 1). These include, a 1.33 times higher impactor flux than the NPF as supported by observations from the LROC NAC temporal image pairs (PF 2)15, an increase of impact flux at 400 Ma with a factor of 3.7 as inferred from the age of impacts spherules (PF 3)17, and an increase of impact flux at 290 Ma with a factor of 2.6 as revealed by the ages of rocky craters (PF 4)54 (Supplementary Fig. 15). The probability distributions of regolith thickness across the study area modeled using these PFs are shown in Fig. 5a. The median regolith thicknesses are 3.45 ± 0.17 m (PF 2), 3.66 ± 0.18 m (PF 3), and 3.22 ± 0.15 m (PF 4), respectively (the dashed lines in Fig. 5a), all of which align with the observation within the 1σ uncertainty range. This indicates that the small crater production and the corresponding impactor flux (<10s of meters) is most probably higher than that predicted by the NPF.

Such an underestimation of the small impactor flux could probably result from the scale-dependent nature of crater degradation, which is more pronounced for small impact craters10,11. Previous studies showed that the lifetimes of impact craters with diameters smaller than 20−50 m are no more than 50−125 Ma12,13. This implies that the record of very small impact craters is still incomplete even on young geological units such as craters Aristarchus (164 ± 1.4 Ma)54, Tycho (109 ± 4 Ma)6, and North Ray (50.3 ± 0.8 Ma)6, where crater counting was conducted to extend the diameter range of the NPF down to 10 m1,55. In contrast, the lifetimes of craters larger than several hundred meters are estimated to be several billion years9,13, whose populations are therefore not significantly influenced by degradation in the mare units. Consequently, the lunar chronology function, which was recently recalibrated through the radiometric age of the CE-5 samples and size-frequency distributions of craters with diameters >200 m in the CE-5 region, is minimally affected by the small impactor flux.

We further conduct a series of simulations by considering different cratering rates varying from 1 to 2 times relative to the NPF over the past 2 Ga. Within the 1σ uncertainty, the modeled regolith thickness aligns with that observed when the relative cratering rate ranges between 1.2 and 1.6 (the black line in Fig. 5b). The optimal match occurs at a relative cratering rate of 1.45, which could be the most probable flux of small lunar impactors. The aforementioned PFs 2–4 correspond to equivalent cratering rates of 1.33, 1.53, and 1.23 times the NPF during 2 Ga, respectively (the colored dots in Fig. 5b), all falling within this range. Such an elevated small impactor flux is significant for dating very young craters and volcanic units (for example, irregular mare patches56) and also implies higher regolith overturn and boulder breakdown rates, which can provide critical insights into the geological evolution of the Moon.

Information from the CE-5 regolith thickness alone can only provide evidence for an overall increased impact flux over 2 billion years, without the ability to further determine the timing when it started to increase nor whether the increase factor varied across different crater size ranges. As a preliminary investigation, we explore several possibilities regarding these two aspects. Based on the estimated ages of major source craters around the CE-5 landing site through their degradation states in this study (for example, C1–C12; Supplementary Table 3), the timing of this increase is estimated to be ∼600 Ma (the black line in Fig. 4d). Several independent studies also support an increase in recent impact flux within a similar time interval. For example, the topographic diffusivity, which reflects the flux of micrometeorite bombardment, was estimated to be higher over the recent 700 Ma than during earlier periods9. The ages of impact glass beads in the CE-5 samples were found to be concentrated in the latest ∼600 Ma57. Moreover, seismic observations of impacts on Mars detected by NASA’s InSight also revealed a higher production rate of small craters (∼1.1–3.0 times) at scales of meters to tens of meters at present58,59. If this is the case, then it indicates recent catastrophic asteroid disruption events in the source regions of lunar and martian impactors, for example the breakups of the parent bodies of L chondrites (470 Ma60,61) and the asteroid Baptistina (160 Ma62) in the main belt.

Alternatively, the elevated impact flux may have persisted during the past 2 Ga, with the increase being more pronounced for smaller impactor sizes. This aligns with recent observations that the power-law slope of the crater size-frequency distribution at tens of meters is steeper than that of the NPF (Supplementary Fig. 16). We performed additional simulations assuming two plausible steeper production functions (PFs 5 and 6; Supplementary Text 3 and Supplementary Figs. 17 and 18) and found that the simulated regolith thicknesses under both scenarios can match the observation. This suggests that the elevated impact flux could be skewed towards smaller impactor sizes, which is consistent with current understanding of the source and delivery mechanism of Earth-Moon impactors. That is, the kilometer-sized and smaller Earth-Moon impactors, which are sourced directly from near-Earth objects (NEOs), were fundamentally resupplied from the main belt primarily via the Yarkovsky and YORP (Yarkovsky–O’Keefe–Radzievskii–Paddack) effects and orbital resonance from ∼3.8 Ga to the present63,64. Unlike collisions and gravitational perturbations, the Yarkovsky and YORP effects change the orbits of asteroids through thermal thrust in the long term, and are more pronounced for smaller asteroids64. The preferential delivery of smaller asteroids to the inner Solar System would result in more abundant small impactors in the NEO populations than those in the main belt3, thus leading to a steeper power-law slope for smaller impactor populations.

Nevertheless, even though the CE-5 region represents the most suitable study area currently available, the increased small impactor flux inferred from the regolith thickness at an individual site may still be affected by local stochastic variability. The broader applicability of our conclusions to the entire Moon remains to be validated in the future, for instance through seismic impact detections by the CE-7 seismometer65 and observations of lunar impact flashes during forthcoming missions66. As the accuracy of regolith thickness estimates improves and the understanding of target property diversity deepens, the framework developed in this study can be applied to a wider range of geologically representative regions, allowing for a more comprehensive assessment of the magnitude and temporal evolution of the small impactor flux on the Moon.

Methods
Regolith thickness estimation from small crater morphology
Optical imagery showed that small impact craters (D<250 m) on the lunar surface exhibit four different morphological types: normal (bowl-shaped), central mound, flat-bottomed, and concentric39. Follow-up laboratory impact experiments and numerical simulations revealed that the strength discontinuity between the surface porous layer (for example, the lunar regolith) and coherent substrate is the primary factor controlling the crater morphology40. Specifically, a normal crater forms when the ratio of crater diameter (D) to regolith thickness (T) < 4 ± 0.1 (1σ uncertainty); a concentric crater forms when D/T > 9 ± 0.6 (1σ uncertainty); a central mound or flat-bottomed crater forms when 4 ± 0.1 < D/T < 9 ± 0.640. Therefore, the size-frequency distributions of small, fresh craters with different morphologies can be used to estimate the distribution of regolith thickness in an area.

The morphological types of impact craters can be determined based on their shadow patterns in optical images27. The CE-5 digital orthophoto map (DOM; 0.5 m/pixel) generated from the descent images and four LROC NAC images (0.8−0.9 m/pixel) with different solar incidence and azimuth angles covering the study area (Supplementary Table 2) are primarily used for crater morphology identification. To enhance reliability, high-resolution descent camera images (<0.5 m/pixel) and the generated CE-5 DTM (0.5 m/pixel) are also utilized to assist in crater morphology determination. Impact craters less than ten pixels in diameter are excluded in the statistics because of their indiscernible internal structures and large relative errors in diameter measurements47. Due to topographic degradation, crater morphology may change over time (for example, from concentric to normal), and crater diameter also increases progressively12. Therefore, only fresh craters are considered in this study. They typically exhibit sharp rims, bright halos, and/or ray patterns, with very few superposed smaller craters. Impact craters formed on slopes (for example, superposed on the walls of larger craters) are also excluded because they are susceptible to the mass wasting process and tend to exhibit asymmetrical features.

We used the kernel density estimator and bootstrap method67 to estimate the size-frequency distributions (SFDs) of craters with different morphologies and the confidence intervals (Supplementary Fig. 12). The percentages of normal and concentric craters within each diameter range represent the area fractions of regolith thickness larger than D/(4 ± 0.1) and less than D/(9 ± 0.6), respectively. From all the measured craters of different sizes and morphological types, the distribution of regolith thickness in the study area can be obtained27. The 1σ uncertainty of regolith thickness was calculated by considering both the variation in range of the crater formation condition (D/T) and the SFD error. In principle, the cumulative distribution of regolith thickness can be estimated through both the normal and concentric craters and their results are generally consistent (the red and blue lines in Fig. 1b). However, because of the limited number of fresh craters at large diameters, the distribution exhibits noticeable fluctuations, particularly for concentric craters. Here, we adopt the result calculated from normal craters as the best estimate of median regolith thickness due to its smaller uncertainty.

For this approach, if a concentric, flat-bottomed, or central mound crater is misclassified as a normal crater (due to crater degradation or low-resolution imagery), the proportion of normal craters in a given diameter interval would increase. Because the proportion of normal impact craters indicates regions where the regolith thickness exceeds D=(4 ± 0:1), this could lead to a higher median thickness (for example, 5 m35). In another approach for regolith thickness estimation through the inner and outer diameters of concentric craters44, the inner structures of small concentric craters cannot be recognized when low-resolution imagery is used. This would lead to biased crater samples that include only large concentric craters, resulting in an overestimate of the regolith thickness (e.g., 4–6 m45).

Lunar topography and regolith evolution model (LTREM)
LTREM simulates the formation and evolution of lunar regolith driven by continuous bombardment over geological timescales. It considers the impact fragmentation, ejecta deposition, secondary cratering, and topographic degradation process25,33. In LTREM, the lunar target is characterized as a two-layer structure, with a progressively growing fine-grained regolith layer atop the underlying coherent substrate (mare basalt in this study). Although complex internal layering within the regolith has been revealed by the Apollo drilling cores22, such microscale structure does not influence the overall regolith thickness. The two-layer simplification adopted by LTREM is based on the widely recognized impact experiments39,40 and simulations68, which have been demonstrated to be valid for regolith thickness estimation. For each impact crater, LTREM simulates the formation of its ejecta blanket and breccia lens during the excavation and modification stages, respectively. During this process, the existing regolith is redistributed, and new regolith may be produced if the transient cavity can shatter the underlying basalt. Generally, LTREM consists of three steps that are repeated for each time interval during the simulation: (1) crater production determination; (2) regolith generation and migration; (3) topographic degradation. Through LTREM, one can obtain the time-varying growth rate and spatial distribution of regolith thickness in an area of interest.

In this study, we employ an observation-constrained LTREM by considering both the global and regional scale impact events that affected the regolith evolution in a 2×2 km area across the CE-5 sampling site. Specifically, the crater populations consist of three parts: the observed source craters delivering materials to the CE-5 region (Fig. 2), invisible craters due to degradation, and secondary craters formed by primaries. Generally, we regard the observed source craters as primary craters because their SFD at the hectometer-scale aligns well with the NPF (Fig. 2d). Taking a minimum ejecta thickness of 0.1 mm as the screening threshold, a total of 956 source craters outside the simulation area38 are considered in the simulation. Within the study area, all the 1223 visible impact craters are counted and taken as model inputs (Fig. 2c). These observed source craters have known locations and diameters, but most have no age information except for the 28 major contributors listed in Supplementary Table 3. The undated craters would form at random times during a geological period from 2 Ga to the present. The populations of invisible craters are determined according to the crater production model extrapolated down to a minimum diameter of 2 m (Supplementary Text 4 and Supplementary Figs. 14b and 19).

We note that the slope of the production function for craters <10 m may be steeper than that of the extrapolated NPF69,70. However, craters in this size range mostly garden rather than producing new regolith, and thus have a negligible effect on the simulated median regolith thickness. They would form at random times from 2 Ga to the present following the probability distribution of the selected impact flux model. Secondary craters were modeled based on the mass-frequency distribution of ejecta fragments and the scaling law of crater size49 (see Supplementary Text 5 for details), with key parameters validated through secondaries produced by the Orientale basin and several other primary craters49. This approach is more refined than using a generic global secondary crater flux (for example, ref. 71), as it accounts for the specific source craters that contributed ejecta to the CE-5 landing region. The key parameter is the power-law slope of the secondary crater SFDs, which is typically −4 as supported by multiple observations72. We also tested other potential values provided by different studies (−4.1415 and −3.573), which will be discussed later in the parameter sensitivity analysis section. Overall, the population of secondary craters used in this study is comparable with that in ref. 74.

Each impact crater is assigned a size-dependent topographic profile (for example, the decreasing depth/diameter ratio versus crater diameter) before being added to the lunar surface. The primary crater profiles are taken from the newly developed shape models of small fresh lunar craters ranging from 36 m to ∼1 km12,75, which include various morphological types and could capture any potential compression during the impact cratering process (Supplementary Text 6 and Supplementary Figs. 20 and 21). The shape models for sub-hectometer craters exhibit a volume-neutral characteristic (Fig. 6 in Zhang et al.25), and therefore small craters incapable of penetrating the existing regolith layer neither increase nor reduce the regolith. Secondary craters are assumed to be parabolic in shape with lower depth/diameter ratios than primaries due to their lower impact velocities (<2.4 km/s), which can be determined according to the distance from their parent primary craters (Supplementary Text 5). When an impact crater forms, we first simulate the opening of the transient crater during the excavation stage and then the collapse during the modification stage. The fragmentation and uplift of bedrock, redistribution of pre-existing regolith, and deposition of impact ejecta are implemented by updating the elevation of two-layer target according to the crater geometry. The crater addition with consideration of the complex pre-impact topographic relief is treated according to cratering dynamics, where hypervelocity impacts essentially erase the original terrain within the crater, whereas the lower-energy ejecta deposition allows the pre-impact surface beyond the rim to be partially preserved (“Crater formation on a two-layer target” in Zhang et al.33). The results of crater superposition using this method have been validated against observations based on LROC NAC DTM data (Fig. 2 in ref. 33).

During the crater formation process, breccia lenses and ejecta blankets on the crater interior and exterior regions are major sources of lunar regolith22, whose thicknesses can be constrained by the topographic profiles of transient and final craters. Since the shape models of final craters were established from high-resolution topographic data (NAC DTMs and SLDEMs), the main uncertainty arises from the shape of the transient crater. Previous laboratory experiments and impact modeling showed that the transient crater can be characterized as a parabolic cavity with a depth/diameter ratio of ∼1/321, which is adopted in our simulation. Different potential depth/diameter ratios are also tested to modify the transient crater shape and will be analyzed later.

The topographic degradation process is simulated using a diffusion model9, expressed as:
∂h/∂t = κ∇2h, (1)
where h is the surface elevation, t is time, and κ is the scale-dependent diffusivity, linked via:
κ/Dφ, (2)
where D is the fresh crater diameter and φ is the scaling factor11, which is typically 0.9. In LTREM, an equivalent diffusivity for the entire simulation area is employed because craters of different diameters are superimposed on each other. However, each location should really have its own diffusivity. According to the median crater size in this study (∼2.5 m), the diffusivity is selected to be 0.031 m2/Myr as determined from the scaling relationship of Eq. (2). We note that applying a uniform equivalent diffusivity across the entire simulation area could lead to an underestimation of the degradation rate for large craters and slightly alter the spatial distribution pattern of regolith thickness. However, different diffusivities tested show that this parameter has a negligible influence on the median regolith thickness (see Parameter Sensitivity Analysis in Methods), as the degradation process primarily redistributes rather than produces lunar regolith.

Model convergence
In our simulation, most of the craters are random in size, location, and/or formation age following the production function. However, although the modeled regolith thickness at a single pixel exhibits strong fluctuations across different realizations, the median thickness across the entire simulation area remains highly stable (Supplementary Fig. 22). This is because the diameters, locations, and ages of major source craters are already known and have been incorporated in our model. For example, Fig. 2 shows 2171 craters with known size and location, among which the 28 major source craters further have age information (Supplementary Table 3). Our observation-constrained Monte Carlo model effectively mitigates the randomness introduced by a large number of impact craters, resulting in outputs that converge within a narrow range. The simulation results show that the median regolith thickness across the study area stabilizes rapidly varying less than 0.01 m after ∼20 realizations (Supplementary Fig. 23a). The relative change induced by each new realization drops to less than 0.1% after ∼30 realizations (Supplementary Fig. 23b). Based on 100,000 bootstrapping samples, the standard error of the median regolith thickness averaged over 100 realizations was estimated to be 0.003 m, and the 95% confidence interval was 2.902–2.912 m (Supplementary Fig. 24). This uncertainty is negligible compared with the error introduced by model parameters (see Parameter Sensitivity Analysis). Other statistical parameters, such as the quartiles, also converge within a narrow range after several tens of realizations. To ensure the statistical significance, all the results provided in this study are the average of at least 100 realizations.

Parameter sensitivity analysis
LTREM includes a series of parameters describing the grid settings, impact flux, crater morphology, and topographic features. Essentially, lunar regolith is produced by repeated impact events, so the regolith thickness is primarily affected by the crater population and regolith production through individual impact craters. The former includes both the primary and secondary crater populations, where the primary impact flux is the parameter to be constrained in this study, and the secondary impact flux is governed mainly by the power-law slope of the secondary populations. The lunar regolith generated by an individual impact can be constrained by the topographic profiles of the transient and final craters. Since the shapes of final craters were directly extracted from the elevation data of fresh craters12,75, the transient crater depth is the main factor influencing regolith production. Topographic diffusivity controls the rate of lateral transport of lunar regolith and hence its spatial distribution, and therefore could also make a difference to the median regolith thickness. In summary, the SFD slope of secondary craters, transient crater depth, and topographic diffusivity are the three main factors influencing the simulated regolith thickness.

We note that target properties (for example, bulk density and strength) can also affect the size of each impact crater and thus the volume of regolith produced (that is, the volume of shattered bedrock, see ref. 25 for detail). Previous studies showed that the rocky crater populations and rock abundance can act as probes for the mechanical properties of the protolith36,37. Both indicate that the P58 unit does not exhibit any anomalies in target mechanical properties. Specifically, the frequency of rocky craters in P58 unit (3.6 × 10−3 km−2) is close to the median value (4.0 × 10−3 km−2) and falls within the interquartile range (3.2 × 10−3 – 4.8 × 10−3 km−2) among 15 different mare units across Oceanus Procellarum, Mare Imbrium, Mare Nubium, Mare Tranquillitatis, Mare Humorum, and Mare Serenitatis36. From our analysis, the median rock abundance of the P58 unit (0.34%) also aligns well with the temporal trend derived from 291 mare units across the Moon (0.33 ± 0.08%; Supplementary Fig. 8). Moreover, the size-frequency distributions of craters are directly taken as model inputs in this study and they are not derived from impactor populations and target properties using the scaling law76. Given that the crater populations (D>200 m) in the CE-5 region and the radiometric age of the CE-5 samples align well with the lunar chronology function7, the target properties of the P58/EM4 unit are expected to be broadly consistent with those of other mare basalts, for example, the Apollo and Luna landing regions.

We first allow a single parameter to vary within its potential range while fixing other parameters to analyze its influence on the regolith thickness. As shown in Supplementary Table 4, a steeper (−4.1415) secondary SFD slope can increase the regolith thickness by 1.0%, and a shallower slope reduces the regolith thickness by 2.2%. When the transient crater depth/diameter ratio ranges from 3/8 to 3/1021, the median regolith thickness varies by up to 8.8%. The effect of topographic diffusivity is relatively minor, as setting it to zero or increasing it by a factor of 10 only leads to a variation of no more than 1.1% in the median regolith thickness.

We further estimate the overall uncertainty of the LTREM using the Monte Carlo method by simultaneously allowing the three parameters to take random values throughout the potential parameter space. The median regolith thicknesses from 1000 simulations are shown in Supplementary Fig. 25 with a minimum of 2.56 m and a maximum of 3.20 m, respectively. However, most results are concentrated within a narrow range with one standard deviation of 0.14 m, corresponding to a relative model uncertainty of 4.8%.

Statistical significance of the difference between observation and simulation
The observed and modeled regolith thicknesses are 3.6 ± 0.3 m and 2.9 ± 0.14 m (1σ uncertainty), respectively. Here, we employed the effect size and bootstrap method to evaluate the statistical significance of the difference. Effect size is a quantitative measure of the magnitude of the difference between two independent groups, which is defined as:
d = (x1 - x2) / sqrt((σ1^2 + σ2^2) / 2), (3)
where x1 and x2, and σ1 and σ2 are the means and standard deviations of two independent groups. According to commonly used guidelines for interpretation77, an effect size greater than 0.2, 0.5, and 0.8 is considered small, medium, and large, and values above 2 are regarded as huge. Based on Eq. (3), the effect size of the difference between the observed and modeled regolith thicknesses is 2.11, exceeding the uncertainties by more than two standard deviations. Besides, from 100,000 bootstrapping samples, the 95% confidence interval (2σ level) of the difference between the observation and simulation is estimated to be 0.06–1.35 m (Supplementary Fig. 26), which is completely above zero. In conclusion, the difference between the observed and modeled regolith thicknesses is statistically significant at 2σ level.

Dating of individual impact craters
For craters larger than 1 km in diameter, their ages are taken from previous studies45 through CSFD measurements on the ejecta blanket regions. For impact craters smaller than several hundred meters, their ages were estimated based on the degradation states. According to the morphological prominence (for example, the sharpness of the crater rim, depth/diameter ratio, inner wall slope), craters with different degradation states can be classified into three major types13: A (fresh), B (moderately degraded), and C (heavily degraded), and two transitional types: AB and BC. The ages of craters with different degradation states can be estimated according to their types: A and AB (<3% of lifetime), B (3−20% of lifetime), BC (20−50% of lifetime), and C (>50% of lifetime). The lifetime (t_life, in Myr) of a crater refers to its survival time on the lunar surface (from formation to erasure), which can be calculated as13:
t_life = 2.5D (D<160 m); 8D^0.900 (160 m < D < 1.2 km) (4)
where D is the crater diameter in meters.

For validation, we conducted a detailed degradation modeling for the 13 small impact craters listed in Supplementary Table 3 (C1−C12 and the unnamed crater). The topographic evolution of each crater is simulated based on the diffusion model (Eq. 1) using a scale-dependent diffusivity (Eq. 2). The optimal estimate of the crater age can be obtained when the modeled depth/diameter ratio matches the observed value extracted from the NAC DTM. We tested a wider range of potential initial diameters ranging from 0.625 D to 1 D, with both the initial crater topographic profiles of Yang et al.12 and Cai and Fa75 considered. As shown in Supplementary Table 5, the ages derived from modeling are generally consistent with those estimated through the degradation states.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Please see attached review.", "reviewer_2": "Summary of review: The manuscript presents the application of a regolith depth evolution model to the Chang’E 5 landing site. While the model application has merit pending validation, the key conclusion concerning the impact rate relies on an unsupported correlation between regolith thickness and surface age. Because of this, I do not believe that it satisfies the goal of Nature Communications “to publish rigorous and comprehensive papers that represent advances of significance to specialists within each field.”\n\nReview:\n1) The state of scientific knowledge does not support the fundamental assumption made in this manuscript that regolith depth is clearly and directly linked to the impact flux. The manuscript claims a “strong correlation between regolith thickness and age across different geological units 18,23” (Lines 64-65); however, this statement is not supported by the literature and is not supported by the provided citations.\nRecent studies examining the abundance of excavated bedrock on mare surfaces indicate a poor relationship between mare regolith depth and time (Chertok et al., 2023; Elder et al., 2023). For example, the abundance of rocky craters suggests that Mare Humorum has the shallowest regolith despite being the oldest mare in the study by Chertok et al. (2023). In general, the depth of the rock-free regolith shows significant variability, which is not correlated to the age of the surface (e.g., Elder et al. 2019; Venkatraman et al. 2023). Rather than developing due to fluctuations in the impact rate over time, regolith depth in the maria may be dominantly controlled by the material properties of the bedrock, as was proposed by Chertok et al. (2023), Elder et al. (2023) and also suggested by Zhang et al. (2022) who stated, “Target properties greatly influence the crater size and hence the volume of regolith produced by impact cratering. For subhectometer impact craters, the strength of bedrock plays a prominent role.”\nRemote sensing radar also illuminates the poor correlation between observed regolith depth and surface age. For example, Shkuratov & Bondarenko (2000) reported that Sinus Rorus and Mare Nectarus have similar regolith depths despite having drastically different surface ages (Sinus Rorus is < 2 Ga, while Nectarus is ~ 4 Ga; Hiesinger et al. (2011)). Further, regolith depths derived from crater morphology (Bart et al. 2011) from various maria suggest no clear trend in mare regolith depth with time. The regolith depths reported in this manuscript (Lines 158-160) fit within the existing picture of a poor correlation between surface age (and, by analogy, impact flux) and regolith depth.\nThe 3-D regolith evolution model utilized in this manuscript is not validated against observations of regolith depth as a function of time through the lens of morphology (e.g., Bart et al. 2011), radar (e.g., Shkuratov & Bondarenko 2000) or rock abundance (e.g., Elder et al. 2023; Chertok et al. 2023; Elder et al. 2019; Venkatraman et al. 2023); the evident validation was a reference to the model by Hirabayashi et al. (2018) shown in Zhang et al. (2022); however, although Hirabayashi et al. (2018) validated their model against the Apollo 15 landing site, Hirabayashi et al. did not validate the relationship of regolith growth with time. Comparison to another theoretical model is not validation.\nIn conclusion, because there is no unambiguous statistically significant link between regolith depth and surface age in the lunar maria — the literature suggests surface age and mare regolith depths are, in fact, completely un-correlated — the extrapolation of flux from regolith depth and surface age is fundamentally flawed and requires validation not yet present in the literature nor provided in the manuscript.\n2) It is unclear whether the results presented in this manuscript are statistically significant. In Figure 4b, all points lie within error, and only a 2σ error is shown, rather than the errors inherited from uncertainty in flux, crater morphology, and material properties. The treatment of secondary impacts may have led to the finding of a high cratering rate relative to the NPF. The manuscript utilized an absolute secondary cratering rate significantly lower than that used by Costello et al. (2021), who showed secondary impacts dominate regolith mixing at sub-hectometer scales and validated their model against the Apollo cores (Morris 1978) and the observed flux of ‘splotches’ (Speyerer et al. 2016). Secondary impacts are more abundant than primaries at hectometer scales. An underestimate in their contribution would lead to an overestimation of the relative contribution of primaries and the false conclusion that there are more primaries (an uptick in the flux). Including secondaries in a crater count is a known cause of offset from the NPF. The manuscript did not sufficiently rule out the likely possibility that secondary craters are the culprit in driving offsets from the NPF.", "reviewer_3": "The study, titled “Increased Small Impactor Flux on the Moon as Inferred from Regolith Thickness at the Chang’E-5 Region” by Mingwei Zhang et al., is very well-written and brought a novel approach to understand the Moon’s impact flux. The research goal is clear; it aims to determine the main driver of unusual regolith thickness constrained from the observed morphology of small craters from Chang’E 5 landing site. The research method volumetrically calculates the contribution of regolith during cratering under different target properties, allowing the authors to separate out the influences of small primary craters, invisible primary craters, secondary craters, and topographic diffusion. By using a conventional chronology functon(Neukum) where a constant impact flux is assumed as constant in the past 2 billion years, the modeled regolith thickness on average with parameter tuning (secondary population and diffusivity) cannot match the regolith thickness constrained by observation data, therefore impact flux in <2 Ga is required to increase for thickening regolith and unusual growth. However, the observed drilling core and seismic data shows a more complicated layer structure of regolith, instead of the underlying assumption that uses two layers model to classify the crater morphology. How the Zhang et al. model tackles this complexity (regolith dynamics and interaction with cratering) is unclear in the current manuscript. In addition, the crater population in Chang’E 5 site is fitted with Neukum function and seems relatively consistent with isotopically dated ages of mare basalt (2 Ga). In the current manuscript, it is not clear that how the crater population with increased impact flux in the model differs from the observed crater population. Some concerns in terms of model capability and stochastic nature of cratering have to be addressed for the finding: the increased impact flux is recorded in the regolith thickness data at the Chang’E 5 landing site. Please see the annotated document for more details and suggestions to improve the clarity of the manuscript.\nWarm regards,\nYa Huei Huang\n2025/3/5"}, "author_response": "Answers to Reviewer 1\nWe appreciate the reviewer for their positive review of our manuscript. We carefully revised\nour manuscript according to the reviewer’s comments. Specifically, we discussed the effect of\ntarget mechanical properties, updated the references, and corrected the typos. A point-by-point\nresponse is provided below.\nComments from the reviewer are in black, and line numbers in the reviewer’s comments refer\nto those in our original manuscript.\nReplies from the authors are in red, and line numbers in our replies refer to those in the track\nchanges revision. References are provided at the end of the file.\nQ1. This manuscript briefly discusses variation in mechanical properties at the Chang’E-5 site,\nand concludes that the site is relatively homogenous. However, it does not discuss other ideas\nabout mechanical property variations and their potential effects on the conclusions presented\nhere. For example, how would the ideas presented in Chertok et al. 2023 (Journal of\nGeophysical Research: Planets, 128, e2023JE007767. https://doi. org/10.1029/2023JE007767) \nrelate to the assumptions made here about the scaling of primary craters? What if the\nmechanical properties of the mare at the Chang’E-5 landing site differ significantly from those\nat other maria? This could potentially weaken the conclusion that the Moon experienced a\nvariation in the flux of small impactors. This should at least be discussed as a potential source\nof error or fodder for future examination.\nThanks for raising this issue. We agree that the target mechanical properties are a non-negligible\nfactor that deserves careful consideration.\nBased on the analysis of Chertok et al. (2023), the frequency (3.6×10-3 km-2) of rocky craters in\nthe P58 unit (CE-5 region) is quite close to the median value (4.0×10-3 km-2) and falls within\nthe interquartile range (3.2×10-3–4.8×10-3 km-2) among 15 different mare units across Oceanus\nProcellarum, Mare Imbrium, Mare Nubium, Mare Tranquillitatis, Mare Humorum, and Mare\nSerenitatis (see their Table 1). This indicates that the P58 unit exhibits no anomalies in target\nmechanical properties. We describe this in L137–141 and L589–593 in the revision.\nElder et al. (2023) showed that rock abundance can act as a probe of the target mechanical\nproperties. The median rock abundance of the P58 unit (0.34%) aligns well with the temporal\ntrend derived from 291 mare units across the Moon (0.33±0.08%), indicating that the P58 unit\nis representative of typical mare basalts. We added Supplementary Fig. 2 and describe this in\nL141–144 and L593–596.\nAlso, the crater population for D>200 m in the P58 unit and the radiometric age of the CE-5\nsamples align well with the lunar chronology function (Yue et al. 2022), which was derived\nusing data from much broader regions across the Moon (Neukum et al. 2001). This agreement\nalso indicates that the mechanical properties of the P58 unit are broadly consistent with other\nmare basalts, for example the Apollo and Luna landing regions. This has been described in\nQ2. In Extended Data Table 2, reference 45 is incorrect for the ages of Copernicus and\nAristarchus. In reference 45, those ages were taken from other references, and the primary\nreferences should be used here.\nWe now provide the primary references for these two craters.\n• Copernicus crater:\nStöffler, D. & Ryder, G. Stratigraphy and isotope ages of lunar geologic units: Chronological\nstandard for the inner solar system. Space Sci. Rev. 1, 9–54, (2001).\n• Aristarchus crater:\nZanetti, M. et al. Evidence for self-secondary cratering of Copernican-age continuous\nejecta deposits on the Moon. Icarus 298, 64–77, (2017).\nQ3. Abstract: missing words “higher than” or similar in this sentence: “The large median\nregolith thickness of 3.6±0.3 m derived from crater morphology can be attributed solely to a\nhigher production rate of small lunar craters at sub-hectometer scales, which is 1.2−1.7 times\nthan previously estimated.”\nWe added higher as suggested (See L27).\nQ4. P4: The word “without” in this sentence is confusing; do you mean, “because they don’t\nshow degradation?” “Second, most large impact craters are identifiable without heavy\ndegradation, whose population can directly serve as model input in regolith evolution\nsimulations”\nHere we mean the large source craters are only slightly degraded. In the revision, we changed\nthe sentence to “Most large source craters that contributed to regolith formation are only slightly\ndegraded and thus remain clearly identifiable.” for clarification. See L107–108.\nQ5. P14: “Because the proportion of normal impact craters indicates regions where the regolith\nthickness exceeds (D/40.1)” should, I think, be D/(4±0.1).\nWe double-checked the sentence and it was correct in the original manuscript (4±0.1). See L453.\nQ6. P17: “stabilizes” is misspelled.\nWe corrected this typo in the revision.\nQ7. Caption for Extended data Fig. 6: “centered” is misspelled.\n“Centred” is changed to “centered”.\nWe appreciate the referee for evaluating our manuscript. We carefully revised our manuscript\naccording to the reviewer’s comments. Specifically, we provided strong evidence for the\ncorrelation between regolith thickness and age, discussed the effect of target properties, and\nperformed statistical analyses to support our conclusions. A point-by-point response is provided\nbelow.\nComments from the reviewer are in black, and line numbers in the reviewer’s comments refer\nto those in our original manuscript.\nReplies from the authors are in red, and line numbers in our replies refer to those in the track\nchange revision. References are provided at the end of the file.\nQ1. The state of scientific knowledge does not support the fundamental assumption made in\nthis manuscript that regolith depth is clearly and directly linked to the impact flux. The\nmanuscript claims a “strong correlation between regolith thickness and age across different\ngeological units 18, 23” (Lines 64-65); however, this statement is not supported by the literature\nand is not supported by the provided citations.\nRecent studies examining the abundance of excavated bedrock on mare surfaces indicate a poor\nrelationship between mare regolith depth and time (Chertok et al., 2023; Elder et al., 2023). For\nexample, the abundance of rocky craters suggests that Mare Humorum has the shallowest\nregolith despite being the oldest mare in the study by Chertok et al. (2023). In general, the depth\nof the rock-free regolith shows significant variability, which is not correlated to the age of the\nsurface (e.g., Elder et al. 2019; Venkatraman et al. 2023). Rather than developing due to\nfluctuations in the impact rate over time, regolith depth in the maria may be dominantly\ncontrolled by the material properties of the bedrock, as was proposed by Chertok et al. (2023),\nElder et al. (2023) and also suggested by Zhang et al. (2022) who stated, “Target properties\ngreatly influence the crater size and hence the volume of regolith produced by impact cratering.\nFor subhectometer impact craters, the strength of bedrock plays a prominent role.”\nRemote sensing radar also illuminates the poor correlation between observed regolith depth and\nsurface age. For example, Shkuratov & Bondarenko (2000) reported that Sinus Rorus and Mare\nNectarus have similar regolith depths despite having drastically different surface ages (Sinus\nRorus is < 2 Ga, while Nectarus is ~ 4 Ga; Hiesinger et al. (2011)). Further, regolith depths\nderived from crater morphology (Bart et al. 2011) from various maria suggest no clear trend in\nmare regolith depth with time. The regolith depths reported in this manuscript (Lines 158-160)\nfit within the existing picture of a poor correlation between surface age (and, by analogy, impact\nflux) and regolith depth.\nThe 3-D regolith evolution model utilized in this manuscript is not validated against\nobservations of regolith depth as a function of time through the lens of morphology (e.g., Bart\net al. 2011), radar (e.g., Shkuratov & Bondarenko 2000) or rock abundance (e.g., Elder et al.\n3\n2023; Chertok et al. 2023; Elder et al. 2019; Venkatraman et al. 2023); the evident validation\nwas a reference to the model by Hirabayashi et al. (2018) shown in Zhang et al. (2022); however,\nalthough Hirabayashi et al. (2018) validated their model against the Apollo 15 landing site,\nHirabayashi et al. did not validate the relationship of regolith growth with time. Comparison to\nanother theoretical model is not validation.\nIn conclusion, because there is no unambiguous statistically significant link between regolith\ndepth and surface age in the lunar maria — the literature suggests surface age and mare regolith\ndepths are, in fact, completely un-correlated — the extrapolation of flux from regolith depth\nand surface age is fundamentally flawed and requires validation not yet present in the literature\nnor provided in the manuscript.\nWe respectfully disagree with the reviewer’s assertion that “regolith thickness is completely\nuncorrelated with surface age”. However, we acknowledge that the physical properties of target\nmaterial are non-negligible factors that can influence the regolith growth process and should be\nconsidered in our analysis. In principle, lunar regolith is primarily formed through mechanical\nfragmentation caused by continuous impact events, and thus its thickness is directly linked to\nthe impact flux. Other factors, including the physical properties of the target material, could\ninfluence regolith thickness to some extent (as indicated by certain outliers). However, these\noutliers do not negate the overall trend of increasing regolith thickness with time.\nOur study uses the crater size-frequency distribution derived from the NPF for both the input\nmodel and the benchmark comparison. Provided that the P58 unit has material properties that\nare representative of mare regions used to construct the NPF, then our approach will by\nconstruction avoid any dependence on material properties. In this response, we demonstrate\nthat our study area (P58) is highly representative and does not exhibit any anomalous target\nmaterial properties. Below, we provide a detailed explanation from the following perspectives:\nthe mechanism of regolith formation, factors influencing regolith growth, supporting\nobservational evidence, and most importantly, the representativeness of our study area in terms\nof target properties.\n1. Mechanism of regolith growth. It is widely accepted that lunar regolith is primarily formed\nby meteorite impacts on the lunar surface (e.g., Oberbeck et al. 1973; Melosh 1989; Mckay et\nal. 1991; Hirabayashi et al. 2018; Xie et al. 2021). Shortly after the emplacement of basaltic\nmare lava flows, little to no regolith is present. Meteorite bombardment by both small and large\nimpactors can shatter the protolith and produce regolith. As regolith thickens, small impactors\nprimarily overturn the existing regolith layer, while only larger impactors capable of penetrating\nthe pre-existing regolith can fracture the underlying substrate and create new regolith.\nConsequently, for a given geological unit on the lunar surface, regolith thickness increases over\ntime, with a gradually decreasing growth rate (Melosh 1989; Mckay et al. 1991; Head and\nWilson 2020; Hirabayashi et al. 2018). Based on the formation mechanism, a number of studies\n4\nmodeled the regolith growth process and their results supported the correlation between regolith\nthickness and time. Here are several examples:\n• Based on the observation that regolith thickness correlates with crater abundance at the\nSurveyor landing sites, Shoemaker et al. (1969) characterized the regolith thickness as a\nfunction of the equilibrium diameter of impact craters (a proxy of regional impact flux).\n• Oberbeck et al. (1973) modeled regolith growth using different crater populations as inputs,\nand their results showed that the regolith thickness increases roughly as (time)0.7 (see their\nFig. 8).\n• Hirabayashi et al. (2018) developed an analytical model describing the growth of regolith\nthickness with time and impact flux (Fig. 11b), and calibrated their model using the result\nof the seismic experiment at the Apollo 15 landing site.\n• Vanga et al. (2022) demonstrated that the rock abundance of intercrater regions in the\nmaria decreases with older ages, which is consistent with regolith thickening over time.\nThey further verified this using a Monte Carlo simulation (their Fig. 3).\n• Zhang et al. (2023a, b) investigated the formation of regolith through individual impact\ncraters and the growth of regolith thickness due to continuous bombardment, which also\nshowed the correlation of regolith thickness with age (Fig. 10 in Zhang et al. 2023a).\n2. Other factors affecting the regolith growth process. When comparing regolith thicknesses\nacross different geological units on the Moon, several additional factors may influence the\nmeasurements and obscure the fundamental relationship between regolith thickness and surface\nage. These factors include: (1) the target properties as proposed by the reviewer; (2) the\ncontribution from individual large craters, which are highly stochastic in both space and time;\nand (3) the buffering trend in the regolith growth process. For these reasons, there could be\nsome outliers that do not strictly follow the correlation between regolith thickness and age. For\nexample, the regolith thickness is abnormally thin in Mare Australe, with a median thickness\nof 2 m despite of its older age of 3.88 Ga (Atang and Bart 2024), which is likely due to its\nanomalous mechanical properties (Elder et al. 2023). Moreover, the regolith growth rate\nbecomes slower in older units due to the buffering trend. As a result, the contribution from\nimpact flux diminishes over time, while the influence of other factors becomes more\npronounced. Therefore, to capture the variation in impact flux through regolith thickness, it is\nessential to focus on young surfaces, such as the CE-5 landing region investigated in this study.\nFinally, we want to emphasize that both factors (2) and (3) have been explicitly incorporated\ninto our model.\n3. Observational evidence from literature. Multiple lines of observational evidence support\na correlation between regolith thickness and surface age, which remains discernible despite the\ninfluence of the aforementioned factors.\n• Based on the crater morphology, Oberbeck and Quaide (1968) estimated regolith\n5\nthicknesses in 12 Lunar Orbiter sites, finding that “it has been demonstrated that\ndifferences in thickness of the regolith from place to place reflect differences in time since\nthe formation of new rock surfaces in those areas.” This finding was summarized in\nMelosh (1989) (Table 10.1), Melosh (2011) (p. 284), and also documented by McKay et\nal. (1991) (Lunar Sourcebook), where the authors stated “regolith thickness is related to\nage” (p.343 in Chapter 7).\n• Using Earth-based 70-cm Arecibo radar data, Shkuratov and Bondarenko (2001) estimated\nregolith thickness across the lunar nearside. Their results demonstrated a strong correlation\nbetween regolith thickness and surface age (see Fig. 11). As the referee noted, the regolith\nthicknesses at Sinus Roris (<2 Ga) and Mare Nectaris (~4 Ga) are 7.5 m and 9.6 m,\nrespectively, and the difference between them may appear small. However, this individual\ncase does not undermine the overall correlation between regolith thickness and age, as\nregolith growth tends to slow down with increasing thickness known as the buffering trend.\n• Since 2009, the high-resolution optical images acquired by LROC NAC enabled a more\nreliable estimate on regolith thickness. Bart et al. (2011) estimated the regolith thicknesses\nat 30 different regions across the Moon, finding that “regolith thickness correlates quite\nwell with crater density”. A subsequent reanalysis by Bart (2025) of these data (Fig. 4)\nalso showed a clear trend of increasing regolith thickness with surface age. To\nquantitatively assess the monotonic relationship between regolith thickness and surface\nage reported by Bart et al. (2011), we applied the Mann–Kendall trend test to their data\n(blue dots in Fig. 3 in Bart 2025). The Mann–Kendall trend test is a non-parametric method\nwidely used to detect statistically significant trends in geophysical data, which does not\nassume any particular data distribution (Hussain and Mahmud 2019). The test result\n(p=0.029<0.05) indicates a statistically significant increasing trend in regolith thickness\nwith increasing surface age, suggesting that older surfaces tend to develop thicker regolith\nlayers with a confidence level exceeding 95%. In addition to this study, Fa et al. (2014)\nfound that the regolith thickness in four geological units in Sinus Iridum increases with\nage (Fig. 12), though the increase rate is small (note that their surface ages are older than\n2.9 Ga, and therefore the buffering trend is more pronounced). Xie et al. (2021)\nsystematically investigated the regolith thicknesses at 13 different geological units on the\nMoon and revealed that regolith thickness correlated well with surface age (Fig. 2).\nWe have carefully reviewed the references provided by the reviewer (i.e., Chertok et al., 2023;\nElder et al., 2023; Elder et al. 2019; Venkatraman et al. 2023). However, these studies focused\non the variations in target properties and did not investigate or refute the correlation between\nregolith thickness and age.\n• Chertok et al. (2023) showed that the difference in rocky crater populations may act as a\nprobe for the mechanical strength of the protolith and the thickness of buried lava flows.\nbut merely suggested that the target property constitutes an additional factor influencing\nregolith development.\n• Elder et al. (2023) found that surface rock abundance decreases with unit age as expected\nfor a thickening regolith and stated “this implies that regolith thickens with time, resulting\nin fewer impacts capable of excavating rocks on older surfaces”. This finding supports the\ngeneral trend that regolith thickness increases over time. The regional variability was\nattributed to variations in target mechanical properties of the mare basalts — a factor that\ncould affect the regolith growth process. For example, that Mare Humorum and southern\nOceanus Procellarum are anomalously rocky, and Mare Australe is rock-poor for its age.\n• Elder et al. (2019) analyzed the subsurface rock content based on the rock abundance in\nthe ejecta blankets of cold-spot craters. They did not investigate the relation between\nregolith thickness and age.\n• Venkatraman et al. (2023) estimated the regolith thickness using rocky craters, but they\ndid not investigate the relation between regolith thickness and age.\n4. Representativeness of the P58 unit. As the youngest isotopically dated mare basalt unit on\nthe Moon, the P58 unit (CE-5 region) is particularly suitable for estimating impact flux through\nregolith thickness. In particular, P58 is geologically representative and does not show any\nanomalous target material properties.\n• Based on the analysis of Chertok et al. (2023), the frequency (3.6×10-3 km-2) of rocky\ncraters in the P58 unit is quite close to the median value (4.0×10-3 km-2) and falls within\nthe interquartile range (3.2×10-3–4.8×10-3 km-2) among 15 different mare units across\nOceanus Procellarum, Mare Imbrium, Mare Nubium, Mare Tranquillitatis, Mare\nHumorum, and Mare Serenitatis (see their Table 1). This indicates that the P58 unit\nexhibits no anomalies in target mechanical properties. We describe this in L137–141 and\nL589–593 in the revision.\n• Elder et al. (2023) showed that rock abundance can act as a probe for target mechanical\nproperties. In the revision, we perform a statistical analysis of rock abundance across 291\nmare units using the updated dataset of Powell et al. (2023). The median rock abundance\nof the P58 unit (0.34%) aligns well with the temporal trend derived from 291 mare units\nacross the Moon (0.33±0.08%), indicating that the P58 unit is representative of typical\nmare basalts. We added Supplementary Fig. 2 and described this in L141–144 and L593–\n596.\n• Additionally, the crater population for D>200 m in the P58 unit and the radiometric age of\nthe CE-5 samples align well with the lunar chronology function (Yue et al. 2022), which\nwas derived using data from much broader regions across the Moon (Neukum et al. 2001).\nThis agreement also indicates that the mechanical properties of the P58 unit are broadly\nconsistent with other mare basalts, for example the Apollo and Luna landing regions. This\n7\nhas been described in\nL598–602 in Methods.\nFinally, we want to emphasize that in this study, the crater size-frequency distribution derived\nfrom the NPF is directly used as model input, rather than being modeled through scaling laws\nthat incorporate the impactor and target properties. This approach inherently avoids the\npotential influence of target material characteristics.\nIn the revision, we clarify the mechanism of the regolith growth process in the introduction,\nprovide supporting observational evidence for the correlation between regolith thickness and\nage from the literature (L64–76), mention other factors that could affect regolith production\n(L76–82), and emphasize the representativeness of the P58 unit in terms of target properties\n(L137–143, L584–602 and Supplementary Fig. 2).\nQ2. It is unclear whether the results presented in this manuscript are statistically significant. In\nFigure 4b, all points lie within error, and only a 2σ error is shown, rather than the errors inherited\nfrom uncertainty in flux, crater morphology, and material properties. The treatment of\nsecondary impacts may have led to the finding of a high cratering rate relative to the NPF. The\nmanuscript utilized an absolute secondary cratering rate significantly lower than that used by\nCostello et al. (2021), who showed secondary impacts dominate regolith mixing at sub-\nhectometer scales and validated their model against the Apollo cores (Morris 1978) and the\nobserved flux of ‘splotches’ (Speyerer et al. 2016). Secondary impacts are more abundant than\nprimaries at hectometer scales. An underestimate in their contribution would lead to an\noverestimation of the relative contribution of primaries and the false conclusion that there are\nmore primaries (an uptick in the flux). Including secondaries in a crater count is a known cause\nof offset from the NPF. The manuscript did not sufficiently rule out the likely possibility that\nsecondary craters are the culprit in driving offsets from the NPF.\n• In our original manuscript, we estimated the uncertainty in regolith thickness by\naccounting for variations in primary and secondary crater populations, crater morphology\n(specifically transient crater depth), and topographic diffusivity, as detailed in the\nparameter sensitivity analysis section of the Methods. We conducted a comprehensive\nsensitivity analysis using 1,000 Monte Carlo simulations.\nIn the revision, we perform additional statistical analyses. Using the widely adopted\nbootstrap method, we estimate the 95% confidence interval for the difference between the\nobserved and modeled regolith thicknesses from 100,000 bootstrap samples, yielding a\nrange of 0.06–1.35 m. This result indicates that the difference is statistically significant\n(Supplementary Fig. 4). We also estimate the standard error and 95% confidence interval\nof the mean value over 100 simulations, which was found to be less than 0.17%\n(Supplementary Fig. 6). See L316–318 and L555–559.\n8\n• We agree with the reviewer that secondary craters play a dominant role in regolith mixing.\nHowever, our study focuses on the growth of regolith thickness rather than the mixing\nof the existing regolith layer. Secondary craters typically exhibit much lower depth-to-\ndiameter ratios (1/50), primarily due to their lower impact velocities (Schultz and Gault\n1985). As a result, they can overturn the surficial regolith layer but rarely generate new\nregolith. We explained this in L306–309.\nThe secondary cratering rate used in our study is comparable with that in O’Brien and\nByrne (2021) but lower than that in Costello et al. (2021). There are two reasons for this:\n(1) The rate employed by Costello et al. (2021) is similar to the population of splotches\nreported by Speyerer et al. (2016). However, the origin of these splotches remains a subject\nof debate. As noted by Speyerer et al. (2016), “these changes in reflectance lacked rims,\nor any other resolvable topographic signature; these splotches merely churn the target\nsurface and do not form classic impact craters”. In Speyerer et al. (2016), the detected\nprimary craters are mostly <20 m in diameter, with a maximum of 43 m, while the\nsplotches range from <10 m up to about 40 m. Given that the maximum diameter of\nsecondary craters is typically no more than 5% of the primaries (supported by numerous\nobservations, e.g., Bierhaus et al. 2018 and also adopted by Costello et al. 2021), the\nobserved splotches with diameters of 10–40 m would require their primaries having\ndiameters ranging from ~200–800 m. Such large and recent impacts would be easily\ndetectable on the lunar surface, yet no such craters have been reported. Therefore, it is\ninappropriate to interpret splotches as typical secondary craters.\nThe primary crater flux adopted by Costello et al. (2021) is only approximately 33%–50%\nof the NPF within 10–100 m diameter range (Fig. 1). Because primary craters excavate to\nsignificantly greater depths than secondary craters, an unusually high flux of secondary\ncraters would be required to achieve an overturn depth comparable with that observed in\nthe Apollo core samples.\ncrater populations used in Costello et al. (2021).\n(2) In addition to the aforementioned issues, Costello et al. (2021) adopted a generalized,\nglobally averaged secondary crater flux. However, the formation of secondary\ncraters—particularly large ones capable of excavating bedrock—is closely linked to\nspecific large primary impacts. In this study, we accounted for the specific source\ncraters that contributed to the CE-5 regolith, and estimated the corresponding\nsecondary crater flux based on their sizes, locations, and ages. The parameters used in\nour secondary crater modeling were adopted from Xie et al. (2020), which have been\nvalidated through observations of secondaries produced by the Orientale basin and\nseveral other primary craters (see their Figs. 2 and 3). Therefore, we think that the\nparameters for secondary craters in our study are more accurate than those used in\nCostello et al. (2021). We described this in L500–503.\nWe appreciate the reviewer’s positive evaluation and constructive comments on our manuscript.\nIn the revision, we carefully revised our manuscript according to the reviewer’s comments.\nSpecifically, we explained the validity of using a two-layer model and performed additional\nstatistical analyses to support our conclusions. A point-by-point response is provided below.\nComments from the reviewer are in black, and line numbers in the reviewer’s comments refer\nto those in our original manuscript.\nReplies from the authors are in red, and line numbers in our replies refer to those in the track\nchange revision References are provided at the end of the file.\nMajor comments\nQ1. However, the observed drilling core and seismic data shows a more complicated layer\nstructure of regolith, instead of the underlying assumption that uses two layers model to classify\nthe crater morphology. How the Zhang et al. model tackles this complexity (regolith dynamics\nand interaction with cratering) is unclear in the current manuscript.\n• We acknowledge that the both the Apollo and Luna core samples revealed a complex\nlayered structure within the regolith, characterized primarily by variations in grain size,\ncomposition, maturity, and noble gas concentrations. However, this microscale internal\nstructure does not impact the method we use to estimate regolith thickness, which is based\non the mechanical properties of the regolith and the underlying bedrock.\n• We do not agree with the referee’s comment that seismic data indicates a more complicated\nlayer structure of regolith. Actually, the Apollo seismic experiments showed a relatively\nhomogeneous regolith structure: “P-wave velocities measured in regolith at the Apollo\nlanding sites are very low and uniform and range from 92 to 114 m/sec” (Carrier et al.\n1991, Lunar Sourcebook, p.516 in chapter 9). The uniform P-wave velocity within the\nregolith, along with its increase at the base of the regolith layer, suggests a homogeneous\nregolith overlying bedrock. This seismic characteristic is the basis for estimating regolith\nthickness (see Nakamura et al., 1975, The Moon, pp. 57–66).\n• The two-layer target model (i.e., a loose granular layer atop a cohesive substrate) has been\nvalidated through both impact experiments and numerical simulations. The most widely\nused method for estimating regolith thickness—small crater morphology—was developed\nthrough impact experiments on a regolith-bedrock two-layer target (Quaide and Oberbeck,\n1968; Oberbeck and Quaide, 1967, 1968). The four observed morphological types of small\nlunar craters were successfully reproduced in numerical impact cratering simulations using\na two-layer structure (e.g., Senft and Stewart, 2007). In addition, comparisons with\nregolith thickness measurements from lunar seismic experiments have demonstrated that\nthis method is reliable and largely unaffected by internal layering within the regolith (e.g.,\napproximation of the lunar surface (i.e., regolith-bedrock) and has long been employed in\nestimating regolith thickness. We explained this in L468–472 in the revision.\nQ2. In addition, the crater population in Chang’E 5 site is fitted with Neukum function and\nseems relatively consistent with isotopically dated ages of mare basalt (2 Ga). In the current\nmanuscript, it is not clear that how the crater population with increased impact flux in the model\ndiffers from the observed crater population.\nWe discussed this in L333–346 in our original manuscript.\nThe diameter ranges of impact craters in Yue et al. (2022) and in our study are different. In Yue\net al. (2022), craters with diameters >200 m were used to fit the chronology function. In contrast,\nour study revealed an increased impact flux for sub-hectometer small craters (<100 m), a\ndiameter range to which regolith thickness is most sensitive (see Extended Data Fig. 5). Within\nthis diameter range, the rapid degradation of impact craters produces a bias between the\nobserved (e.g., the Neukum production function) and the actual crater production.\nMinor comments\nQ3. L63–66. This is perhaps the most critical observation from the manuscript. Please consider\nto add a visual in the main text.\nWe added several sentences describing the growth process of lunar regolith to clarify the\ncorrelation between regolith thickness and age. See L64–76.\nQ4. L79. Please add a brief description on this term (first use). If the “self-limiting growth\nprocess” does occur, why the regolith thickness data can tell us about the impact flux is a bit\nodd and not clear in the current manuscript. For example, one expects to see the presence of\n“self-limiting growth process” at the lunar surface.\nThank you for the suggestion. In the revision, we defined the self-limiting process and the\nbuffering trend in regolith growth (See L68–73). Due to the self-limiting nature of regolith\nformation, its growth is sensitive to impact flux primarily on geologically young lunar surfaces.\nThis is why we selected the CE-5 landing region for this study, as it is the youngest isotopically\ndated mare basalt unit.\nQ5. L94–95. Please clarify the sentence and provide a more appropriate reference. At least,\nreference 29 is the book, so it would be great to know the page numbers that for an interested\nreader to read it.\nThis sentence described the self-limiting growth process. See our response to Q4. This is\ndescribed on p.285–286 of Mckay et al. (1991). We provided the chapter information in the\nrevision. See L119.\ngrowth rates, 1.8 and 1.2, are on the same order regardless the uncertainties in deriving older\nages. Consider to delete or revise the following sentences, “This could result from the self-\nlimiting…” It actually weakens the hypothesis proposed in the study.\nThis can be found on p. 198 of chapter 10 in ref. 21. The mean growth rates at the Apollo\nlanding sites can be calculated based on regolith thicknesses and surface ages, which were\nsummarized in Extended Data Table 1. We cited this table in the revision.\nThis represents a 50% difference, which we consider to be significant. Given that the Apollo\n11 landing region is 3.8 Ga old, the low regolith growth rate reflects the self-limiting process.\nThis supports our choice of the CE-5 landing site, because the relatively young surface is less\naffected by the self-limiting process, making it more suitable for assessing the impact flux.\nQ7. L229–231. Please remind of regolith deposition being static. If the reviewer’s\nunderstanding is correct, this static regolith growth mechanism cannot dynamically capture the\nevolution of regolith. If so, please have this discussion for parameter uncertainties.\nThe percentages reported here are derived from Fig. 3d, where the contributions to regolith\nformation from the observed, invisible, and secondary craters were recorded. As a large number\nof craters form in LTREM, the bedrock is shattered and excavated, and the regolith is deposited\nand transported, which is a dynamic process. Different parameters primarily affect the total\nregolith thickness, while the relative contribution percentages remain largely stable. The 1σ\nuncertainty is 0.1%. We added this in the revision (L282).\nQ8. L225–257. The modeled regolith thickness using NPF, 2.48–3.32 meter (3 sigma\nsignificance) is on the same order of 2.7–4.5 meter though. I suggest to have more statistical\nanalysis using non-biased test on the bootstrapped samples of the 100 realization data (please\nrefer to Huang et al., (2022) supplementary test).\nThanks for the suggestion. In the revision, we performed additional statistical analyses of our\nresults as suggested by the reviewer. We estimated the 95% confidence interval of the difference\nbetween the observed and modeled regolith thicknesses from 100,000 bootstrapping samples,\nwhich is 0.06–1.35 m. This indicates that the larger regolith thickness is statistically significant\n(Supplementary Fig. 4). We also estimated the standard error and 95% confidence interval\nof the mean regolith thickness over 100 realizations, which is less than 0.17% (Supplementary Fig.\n6). These analyses show that our conclusion regarding an increased small impactor flux is\nreliable. See L316–318 and L554–558.\nQ9. L416. The heterogeneity of ray distribution, as well as its ballistic sedimentation,\nhaven’t explored yet, and how its sporadic nature affects SFD is not discussed.\ndistribution of distal ejecta (e.g., Elliott et al., 2018; Minton et al., 2019), our approach provides\na first-order, average approximation for the cumulative effect of secondaries generated by a\nlarge number of primary impacts. Because secondary craters play a dominant role in mixing\nthe surficial regolith rather than producing new regolith, this would not affect our results in any\nsignificant manner. We discussed this in Supplementary Text 1 in the revision. See L47–52 in\nSupplementary Information.\nQ10. L496–497. This sentence sounds another way around; for example, from Minton et al.,\n(2019) the distal ejecta (lateral transport of lunar regolith) drives the change of surface\ntopography. Please consider to revise this sentence and/or to add a reference to support this.\nHere we mean this parameter (diffusivity) controls the rate of regolith transport. This sentence\nis changed to “Topographic diffusivity controls the rate of lateral transport of lunar regolith”.\nSee L575.\nQ11. L499–500. A side note as “the volume of regolith produced” is first appeared here. How\nthe volume of regolith produced in the manuscript is not stated in the main text and the\nsupplementary. Zhang et al. (2023) has a very detailed description. I suggest to have at least the\ndefinition of “the volume of regolith produced” and their components (shattered rock, the\nvolume of regolith changed …) in the main text or the supplementary and refer details to\nZhang et al., (2023).\nWe described this briefly and added the reference in the revision. See L584–585.\nQ12. L506–517. This reviewer’s impression is that the LTREM allows us to break down the\ncomponents of regolith production. With this model capability, the authors can test each\ncomponent that potentially drive the regolith production in Chang’E 5 landing site, and\nfurthermore the authors can test whether impact flux in <2 Gy is increased. If this chain of\nthought is right, please have similar sentences in the main text as it was not clear and strong.\nThanks for the suggestion. We emphasized this chain in the main text. See L223–228.\n• Atang, E. F. M., & Bart, G. D. Lunar Regolith Depth Vs. Surface Age: New Measurements\nin the Lunar Maria. in 55th Annual Lunar and Planetary Science Conference. 1549, (2024).\n• Bart, G. D. Evaluating the Depth–Age Hypothesis for the Evolution of the Lunar Regolith.\nThe Planetary Science Journal, 6(2), 36, (2025).\n• Bierhaus, E. B. et al. Secondary craters and ejecta across the solar system: Populations and\neffects on impact‐crater–based chronologies. Meteorit. Planet. Sci. 53, 638–671, (2018).\n• Carrier III, W. D. et al. in Lunar Source-Book: A User's Guide to the Moon. Physical\nproperties of the lunar surface. (eds G. H. Heiken, D. T. Vaniman, & B. M. French) Ch. 9,\n475–594 Cambridge University Press, (1991).\n• Chertok, M. A., Lucey, P. G., Costello, E. S. & Ireland, S. M. The rock abundance of crater\npopulations as a probe of mare protolith properties. J. Geophys. Res. Planets 128,\ne2023JE007767, (2023).\n• Costello, E. S., Ghent, R. R. & Lucey, P. G. Secondary impact burial and excavation\ngardening on the Moon and the depth to ice in permanent shadow. J. Geophys. Res. Planets\n126, e2021JE006933, (2021).\n• Di, K., et al. Lunar regolith thickness determination from 3D morphology of small fresh\ncraters. Icarus, 267, 12–23, (2016).\n• Elder, C. M. et al. The subsurface coherent rock content of the moon as revealed by cold‐\nspot craters. J. Geophys. Res. Planets 124(12), 3373–3384, (2019).\n• Elder, C. M. et al. The variability of lunar mare basalt properties from surface rock\nabundance. Planet. Sci. J. 4, 244, (2023).\n• Elliott, J. R., Huang, Y.-H., Minton, D. A. & Freed, A. M. The length of lunar crater rays\nexplained using secondary crater scaling. Icarus 312, 231–246, (2018).\n• Fa, W. Z., Liu, T. T., Zhu, M. H. & Haruyama, J. Regolith thickness over Sinus Iridum:\nResults from morphology and size-frequency distribution of small impact craters. Journal\nof Geophysical Research-Planets 119, 1914–1935, (2014).\n• Head, J. W. & Wilson, L. Rethinking lunar mare basalt regolith formation: New concepts\nof lava flow protolith and evolution of regolith thickness and internal structure. Geophys.\nRes. Lett. 47, e2020GL088334, (2020).\n• Hirabayashi, M. et al. The role of breccia lenses in regolith generation from the formation\nof small, simple craters: Application to the Apollo 15 landing site. J. Geophys. Res. Planets\n123, 527–543, (2018).\n• Hussain, M M., & Mahmud, I. PyMannKendall: A python package for non parametric\nMann Kendall family of trend tests. Journal of Open Source Software, 4(39), 1556, (2019).\nhttps://doi.org/10.21105/joss.01556\n• McKay, D. S. et al. in Lunar Source-Book: A User's Guide to the Moon. The Lunar\nRegolith. (eds G. H. Heiken, D. T. Vaniman, & B. M. French) Ch. 7, 285–356 Cambridge\nUniversity Press, (1991).\n15\n• Melosh, H. J. Impact Cratering: A Geologic Process. (Oxford University Press, 1989).\n• Melosh, H. J. Planetary surface processes. (Cambridge University Press, 2011).\n• Minton, D. A., Fassett, C. I., Hirabayashi, M., Howl, B. A. & Richardson, J. E. The\nequilibrium size-frequency distribution of small craters reveals the effects of distal ejecta\non lunar landscape morphology. Icarus 326, 63–87, (2019).\n• Nakamura, Y., Dorman, J., Duennebier, F., Lammlein, D. & Latham, G. Shallow lunar\nstructure determined from the passive seismic experiment. The Moon 13, 57–66, (1975).\n• Neukum, G., Ivanov, B. A. & Hartmann, W. K. Cratering records in the inner solar system\nin relation to the lunar reference system. Space Sci. Rev. 96, 55–86, (2001).\n• Oberbeck, V. R. & Quaide, W. L. Estimated thickness of a fragmental surface layer of\nOceanus Procellarum. J. Geophys. Res. 72, 4697–4704, (1967).\n• Oberbeck, V. R. & Quaide, W. L. Genetic implications of lunar regolith thickness\nvariations. Icarus 9, 446-465, (1968).\n• Oberbeck, V. R., Quaide, W. L., Mahan, M. & Paulson, J. Monte Carlo calculations of\nlunar Regolith thickness distributions. Icarus 19, 87–107, (1973).\n• O'Brien, P., & Byrne, S. Physical and chemical evolution of lunar mare regolith. J.\nGeophys. Res. Planets 126(2), e2020JE006634, (2021)\n• Quaide, W. L. & Oberbeck, V. R. Thickness determinations of the lunar surface layer from\nlunar impact craters. J. Geophys. Res. 73, 5247–5270, (1968).\n• Schultz, P. H. & Gault, D. E. Clustered impacts: Experiments and implications. Journal of\nGeophysical Research: Solid Earth 90, 3701–3732, (1985).\n• Senft, L. E., & Stewart, S. T. Modeling impact cratering in layered surfaces. J. Geophys.\nRes. Planets 112(E11), (2007).\n• Shkuratov, Y. G. & Bondarenko, N. V. Regolith layer thickness mapping of the Moon by\nradar and optical data. Icarus 149, 329–338, (2001).\n• Shoemaker, E. M., et al. Observations of the lunar regolith and the Earth from the\ntelevision camera on Surveyor 7. J. Geophys. Res. Planets, 74(25), 6081–6119, (1969).\n• Speyerer, E. J., Povilaitis, R. Z., Robinson, M. S., Thomas, P. C. & Wagner, R. V.\nQuantifying crater production and regolith overturn on the Moon with temporal imaging.\nNature 538, 215–218, (2016).\n• Stöffler, D. & Ryder, G. Stratigraphy and isotope ages of lunar geologic units:\nChronological standard for the inner solar system. Space Sci. Rev. 1, 9–54, (2001).\n• Vanga, S. et al. Rock abundance on the lunar mare on surfaces of different age:\nImplications for regolith evolution and thickness. Geophys. Res. Lett. 49,\ne2021GL096710, (2022).\n• Venkatraman, J. et al. Statistical estimates of rock-free lunar regolith thickness from\ndiviner. Planetary and Space Science 229, 105662, (2023).\n• Xie, M., Liu, T. & Xu, A. Ballistic sedimentation of impact crater ejecta: Implications for\nthe provenance of lunar samples and the resurfacing effect of ejecta on the lunar surface.\n• Xie, M., Xiao, Z., Xu, L., Fa, W. & Xu, A. Change in the Earth–Moon impactor population\nat about 3.5 billion years ago. Nat. Astron. 5, 128–133, (2021).\n• Yang, X., Fa, W. Z., Du, J., Xie, M. G. & Liu, T. T. Effect of topographic degradation on\nsmall lunar craters: Implications for regolith thickness estimation. Geophys. Res. Lett. 48,\ne2021GL095537, (2021).\n• Yue, Z. et al. Updated lunar cratering chronology model with the radiometric age of\nChang’e-5 samples. Nat. Astron. 6, 541–545, (2022).\n• Zanetti, M. et al. Evidence for self-secondary cratering of Copernican-age continuous\nejecta deposits on the Moon. Icarus 298, 64–77, (2017).\n• Zhang, M., Fa, W. & Eke, V. R. Modeling the evolution of lunar regolith: 1. Formation\nmechanism through individual simple impact craters. J. Geophys. Res. Planets 128,\ne2023JE007850, (2023).\n• Zhang, M., Fa, W., Barnard, E. M. & Eke, V. R. Modeling the evolution of lunar regolith:\n2. Growth rate and spatial distribution. J. Geophys. Res. Planets 128, e2023JE008035,\n(2023).\n17"}, {"round": 2, "reviewer_comments": {"reviewer_2": "[Editorial Note: Please see the final page of this file for the attached figure]\nThank you for addressing my review and for your thorough response. I am satisfied with the manuscript revisions regarding statistical significance and the treatment of secondaries.\nConcern remains regarding a fundamental assumption of a direct and obvious correlation between regolith depth and age. Since my review, Bart (2025) has published a review of observational evidence of regolith depth from diverse sources collected over the last several decades (including those cited in the rebuttal to my previous review). Bart (2025) demonstrates that the observations of regolith depth do not support the depth-age hypothesis. To illustrate this point, in an attached figure, I plot regolith depth versus mare age for the CE'5 regolith depth estimates from this manuscript, two of the data sources included in Bart (2025) (Bart et al. (2011) and Shkuratov & Bondarenko (2001)), and a recently published assessment of regolith depth evolution on crater ejecta (Thanker & Neish 2025). This figure includes error bars, which compound the poor statistical correlation between depth and time identified by Bart (2025). The data show no positive correlation between age and regolith depth. A straight line could fit through all of the data in the figure, including the CE'5 data. Observations are relatively sparse for surfaces younger than 2.5 Ga, and the 70-cm radar scattering model underpinning the Shkuratov & Bondarenko (2001) regolith depth measurements has unresolved issues (Campbell 2002), but the available data do not support a clear relationship between regolith depth and age, even for these young surfaces.\nThe apparent lack of correlation between observed regolith depth and mare age must be reconciled with the statement in the revised manuscript: \"Therefore, regolith thickness is statistically linked to surface age and can serve as a robust indicator of lunar impact history, particularly for a young lunar surface, as supported by both observational and modeling results.\" (L72-73)\nI agree that impacts likely dominate regolith generation, and I agree that the properties of the study landing region are relatively uniform. In light of these agreements and having read the same assertions that were cited in the rebuttal, like the direct quote from the Lunar Sourcebook, I was surprised by the lack of depth-age correlation when assessing the sum of regolith depth observations collected over the last few decades; however, given the apparent poor correlation between observed regolith depth and age in the mare, it is clear that a direct relationship between flux, surface age, and regolith depth cannot be taken as a given.\nConsider discussing the issues with the depth-age hypothesis, focusing on the unique case/location presented in this manuscript, and including more of the material written in the thorough rebuttal to my previous review, particularly concerning the representativeness of the P58 unit. An explicit argument would strengthen the manuscript that while the story of regolith depth and its evolution is complicated and poorly correlated with surface age when viewed across multipule mare (e.g., Bart 2025, Chertok et al. 2023), we can use models to assess individual locations, like Hirabayashi et al. (2018) did for Apollo 15 and this manuscript does for the CE'5 site, because we expect the effects of factors like protolith properties, local large crater ejecta, development of a regolith buffer that is subsequently mixed etc., to be within the scope of model parameterization.\nReferences:\nBart, G. D. (2025). Evaluating the Depth–Age Hypothesis for the Evolution of the Lunar Regolith. The Planetary Science Journal, 6(2), 36.\nThaker, A. D., & Neish, C. D. (2025). Assessing the depth of the regolith around radar-dark halo craters on the moon. Icarus, 116623.\nCampbell, B. A. (2002). Comment on “Regolith layer thickness mapping of the Moon by radar and optical data,” by YG Shkuratov and NV Bondarenko. Icarus, 158(2), 560-561.", "reviewer_3": "The manuscript, titled “Increased Small Impactor Flux on the Moon as Inferred from Regolith Thickness at the Chang’E-5 Region” by Mingwei Zhang et al., maintains a clear and exciting research goal: use the regolith thickness distribution at the Chang'E-5 region to infer the impact flux rate at the recent two billion years on the Moon. The reivewer thanks the Zhang et al. for making effort to explain the LTREM model, however, the model itself still requires clarifications. First of all, the assumption that the relationship between ages and regolith thickness can be relaxed as there remains uncertainties regarding each individual site (e.g., Mare Australe's unusual regolith thickness distribution). The authors ensured the target properties of the Chang’E-5 region with the comparison studies. Ultimately, we would like to have a model that determines the main driver of regolith thickness from the observed morphology of small craters and the crater size frequency distribution at a site. Where the model breaks down, it may uncover an intrinsic properties of the underlying surface or impactor populations.\nThe most concerning question to this reviewer is the quantification of complex interaction in multiple cratering processes. This is an important task as the modeled regolith thickness from LTREM with increased impact flux scenarios can be explained by two standard deviation from the observed regolith thickness distribution at the Chang’E-5 site. For example, it will be interesting to see how two layers' thickness flucaute over time. Is it truly sufficient to have two layers model (LTREM) to explain the observed regolith thickness distribution? How the two layers regolith model can dynamically capture the random nature of cratering process and complexity of interaction between regolith growth and compression, as well as target properties, has to be tested an quantified against more sites.\nPlease see few minor comments in the annotated document regarding the Neukum flux in the model.\nWarm regards,\nYa Huei Huang\n2025/5/26"}, "author_response": "Answers to Reviewer 2\nWe appreciate the reviewer’s careful review and constructive comments. In the revision, we systematically analyzed the correlation between regolith thickness and age based on multi-source remote sensing observations (Supplementary Figs. 1–5), and provided an in-depth discussion on the factors that may affect the thickness–age hypothesis. We agree with the reviewer that this hypothesis may not be universally applicable to all mare basalts across the Moon, and thus we focused on demonstrating the homogeneity and representativeness of the P58 unit.\nComments from the reviewer are in black, and line numbers in the reviewer’s comments refer to those in our original manuscript.\nReplies from the authors are in red, and line numbers in our replies refer to those in the track change revision. References are provided at the end of the file.\nQ1. Concern remains regarding a fundamental assumption of a direct and obvious correlation between regolith depth and age. Since my review, Bart (2025) has published a review of observational evidence of regolith depth from diverse sources collected over the last several decades (including those cited in the rebuttal to my previous review). Bart (2025) demonstrates that the observations of regolith depth do not support the depth-age hypothesis. To illustrate this point, in an attached figure, I plot regolith depth versus mare age for the CE'5 regolith depth estimates from this manuscript, two of the data sources included in Bart (2025) (Bart et al. (2011) and Shkuratov & Bondarenko (2001)), and a recently published assessment of regolith depth evolution on crater ejecta (Thanker & Neish 2025). This figure includes error bars, which compound the poor statistical correlation between depth and time identified by Bart (2025). The data show no positive correlation between age and regolith depth. A straight line could fit through all of the data in the figure, including the CE'5 data. Observations are relatively sparse for surfaces younger than 2.5 Ga, and the 70-cm radar scattering model underpinning the Shkuratov & Bondarenko (2001) regolith depth measurements has unresolved issues (Campbell 2002), but the available data do not support a clear relationship between regolith depth and age, even for these young surfaces.\nThe apparent lack of correlation between observed regolith depth and mare age must be reconciled with the statement in the revised manuscript: \"Therefore, regolith thickness is statistically linked to surface age and can serve as a robust indicator of lunar impact history, particularly for a young lunar surface, as supported by both observational and modeling results.\" (L72-73)\nI agree that impacts likely dominate regolith generation, and I agree that the properties of the study landing region are relatively uniform. In light of these agreements and having read the same assertions that were cited in the rebuttal, like the direct quote from the Lunar Sourcebook, I was surprised by the lack of depth-age correlation when assessing the sum of regolith depth observations collected over the last few decades; however, given the apparent poor correlation between observed regolith depth and age in the mare, it is clear that a direct relationship between flux, surface age, and regolith depth cannot be taken as a given.\nConsider discussing the issues with the depth-age hypothesis, focusing on the unique case/location presented in this manuscript, and including more of the material written in the thorough rebuttal to my previous review, particularly concerning the representativeness of the P58 unit. An explicit argument would strengthen the manuscript that while the story of regolith depth and its evolution is complicated and poorly correlated with surface age when viewed across multiple mare (e.g., Bart 2025, Chertok et al. 2023), we can use models to assess individual locations, like Hirabayashi et al. (2018) did for Apollo 15 and this manuscript does for the CE-5 site, because we expect the effects of factors like protolith properties, local large crater ejecta, development of a regolith buffer that is subsequently mixed etc., to be within the scope of model parameterization.\nThanks for your constructive comment. We agree with the reviewer that “while the story of regolith depth and its evolution is complicated and poorly correlated with surface age when viewed across multiple mare, we can use models to assess individual locations”. Our response is structured into four parts: (1) the attached figure, (2) evaluation of the thickness–age relationship based on observations, (3) factors influencing this relationship, and (4) the representativeness and uniqueness of the CE-5 region.\n1. On the attached figure\nThanks for providing the figure and we appreciate this effort to integrate existing results. In this figure, the regolith thickness does not show a clear correlation with age. However, we want to emphasize a few limitations of this compilation. First, the statistical metrics are not consistent: mean and standard deviation for Shkuratov & Bondarenko (2001) and Thanker & Neish (2025), whereas median and min/max values for Bart et al. (2011). Second, there are inherent systematic differences among different studies due to varying data sources, methodologies, and underlying assumptions. For example, the crater morphology method is based on the strength discontinuity between regolith and bedrock (Bart et al. 2011), whereas radar observations rely on the differences in their dielectrical properties (Shkuratov & Bondarenko 2001). Besides, regolith growth on ejecta deposits (Thanker & Neish 2025) may not be directly comparable with that on mare basalts, because fine-grained materials within ejecta blankets may themselves already be part of the regolith. All these issues could affect the interpretation of the plotted trends.\n2. Regolith thickness vs. age from multi-source observations\nTo systematically evaluate the correlation between regolith thickness and age, we analyzed regolith thickness estimated from radar observations, microwave radiometer measurements, and the crater morphology method, respectively. Then, two statistical tests—the Spearman rank correlation and the Mann-Kendall test—were applied to these results. The Spearman rank correlation is a non-parametric test that assesses the strength and direction of the monotonic relationship between two variables, without assuming a linear relationship or normal distribution. The correlation coefficient (ρ: –1 to 1) indicates how strongly two sets of ranks are correlated (Spearman 1987). The Mann-Kendall test is a non-parametric statistical test used to detect the presence of a monotonic trend in a time series without requiring any particular distribution for the data (Hussain & Mahmud 2019). For both methods, the p value represents the probability of observing the data under the null hypothesis, with a small p value indicating a statistically significant monotonic correlation.\n2.1 Radar observations\nUsing the Earth-based 70-cm Arecibo radar data, Shkuratov and Bondarenko (2001) inverted the regolith thickness over the lunar nearside based on a multiple scattering model of radar waves. We georeferenced the regolith thickness map of Shkuratov and Bondarenko (2001) with the mare basalts dated by Hiesinger et al. (Supplementary Fig. 1). We then calculated the mean regolith thickness within each basalt unit containing more than 100 pixels (400–2500 km2), as shown in Supplementary Fig. 2. The blue line shows the prediction of the analytical regolith growth model (Zhang et al. 2023a), which is applicable for surface younger than 3.5 Ga due to different earlier impactor populations (Xie et al. 2021).\nFor all the dots, the Spearman rank correlation yielded a ρ value of 0.12 with a p-value of 0.09, indicating no significant monotonic relationship between regolith thickness and age. Consistently, the Mann-Kendall test detected no significant trend with a p-value of 0.08. We further performed an analysis on basalts whose target properties showed no significant anomalies (the red dots), defined as rock abundance values falling within one standard deviation in Supplementary Fig. 7. The results also indicate no significant trend.\nIn conclusion, the regolith thickness derived from Shkuratov and Bondarenko (2001) shows no correlation with age. However, this lack of correlation is likely due to problems with the radar data calibration and inversion method (Campbell 2002; Fa and Wieczorek 2012). A detailed discussion including the data calibration, radar scattering model, and polarization properties can be found in the Introduction of Fa and Wieczorek (2012) (second-to-last paragraph of page 772), and is not repeated here.\n20\n15\n10\n2.2 Microwave radiometer measurements\nThe global distribution of lunar regolith thickness was derived from the 3 GHz brightness temperature observed by CE-1 microwave radiometer (Fa and Jin 2010; Supplementary Fig. 3). Similarly, we calculated the mean regolith thickness within each basalt unit containing more than 10 pixels (~9000 km2 at the equator), as shown in Supplementary Fig. 4.\nFor all the dots, the Spearman rank correlation yielded a ρ value of 0.40 with a p-value of 1.2´10-6, indicating a statistically significant moderate correlation between regolith thickness and age. Consistently, the Mann-Kendall test detected an increasing trend with a p-value of 1.0´10-6. The red dots represent the basalts whose target properties showed no significant anomalies and the statistical results for these red dots also indicate a moderate correlation between regolith thickness and age.\nSupplementary Fig. 3 Supplementary Fig. 4\n2.3 Crater morphology method\nBart et al. (2011) performed a global survey of regolith thickness using small crater morphology. We analyzed the relationship between the median regolith thickness and age based on their results, using the data reported in Figure 3 of Bart (2025) (Supplementary Fig. 5). The Spearman rank correlation yielded a ρ value of 0.57 with a p-value of 0.02, indicating a statistically significant correlation between regolith thickness and age. Consistently, the Mann-Kendall test detected an increasing trend with a p-value of 0.03.\n)m(\nssenkciht\nhtiloger\nnaeM\n(a)\n1 2 3 4\nAge (Ga)\n10\n5\n)m(\nssenkciht\nhtiloger\nnaeM\n(a)\n1 2 3 4\nAge (Ga)\n8\n6\n4\nIn summary, the regolith thickness estimated from microwave radiometer measurements (Fa and Jin, 2010) and the crater morphology method (Bart et al., 2011) shows a moderate correlation with age, whereas the thickness estimated from radar observations (Shkuratov and Bondarenko, 2001) shows no such correlation. As already mentioned, there are several issues with the radar-based method of Shkuratov and Bondarenko (2001); therefore, we consider regolith thickness and age to show an overall moderate statistical correlation. Our statement that “regolith thickness is statistically linked to surface age” is based on the statistical tests described above, and does not imply that older surfaces necessarily have thicker regolith than younger surfaces. In fact, we agree with the reviewer’s perspective and the difference lies only in the wording of our description.\nIn the revision, we discussed this correlation thoroughly in the main text and point out that, while there is a moderate statistical correlation, some regions clearly deviate from this trend, in order to avoid any misinterpretation (see L81–91). We also added Supplementary Text 1 and Figs. 1–5 providing the regolith thickness obtained from radar observations, microwave radiometer measurements and crater morphology.\n3. Reasons for the variability in regolith thickness\nSince lunar regolith is primarily formed by impacts (e.g., Mckay et al. 1991), its thickness within a given geological unit is expected to increase over time. However, the regolith thicknesses in mare units with different ages across the Moon exhibit considerable variation, and their correlation with age is not particularly strong (Supplementary Figs. 2, 4 and 5). This may be attributed to the following four factors.\n3.1 Limitations associated with the data and methods\nFor regolith thickness, radar-derived estimates are strongly influenced by the choice of radar wave scattering model and radar data radiometric calibration error (Shkuratov and Bondarenko 2001), which lead to large errors for regolith thickness inversion (Campbell 2002; Fa and Wieczorek 2012). In contrast, the crater morphology method is based on impact experiments (Quaide and Oberbeck 1968; Oberbeck and Quaide 1967, 1968) and further supported by numerical simulations (Senft and Stewart 2007), making it relatively more reliable.\nFor ages of mare basalts, crater-counting results reported by Hiesinger et al. can differ substantially from sample-based ages. For example, the P58 unit was previously estimated to be ~1.3 Ga by crater counting (Hiesinger et al. 2011), whereas isotopic dating of the CE-5 samples yields an age of ~2.0 Ga (Li et al. 2021). For the CE-3 landing region, Hiesinger et al. (2011) estimated an age of 2.96 Ga, whereas Qian et al. (2023) reported an age of 1.4 Ga, both of which used crater counting method. Without isotopic dating, it remains difficult to assess its accurate age. Similar discrepancies are expected in other mare units as well.\nThe errors in both regolith thickness and age estimates can weaken the apparent correlation between them.\n3.2 Variation of target property\nRecent observations of both rocky crater populations and surface rock abundances indicate that mechanical properties vary substantially across different mare basalt units (e.g., Elder et al., 2023; Chertok et al., 2023). The strength and bulk density of the target material can influence crater sizes and thus affect regolith production rates (e.g., Zhang et al., 2023a). This issue has been discussed in the reviewer’s comments and our previous response.\n3.3 Random nature of large impacts\nLunar regolith is primarily formed by local small impacts but is also influenced by exotic ejecta deposits (e.g., Mckay et al. 1991; Xie et al. 2020a). Numerous small craters can be considered spatially uniform, whereas large craters are relatively rare and exhibit highly stochastic spatial patterns. Areas affected by ejecta and secondary craters from large impacts tend to have thicker regolith. The spatial randomness and temporal episodicity of large impacts result in regolith thicknesses that do not increase steadily with time across different geological units (e.g., crater Xu Guangqi in the simulation area and crater Copernicus outside the area).\n3.4 Buffering trend in regolith growth\nAs regolith thickens, small impactors merely overturn the existing regolith layer, while only larger impactors capable of penetrating the pre-existing regolith can fracture the underlying substrate and produce new regolith. As a result, regolith growth rate decreases over time and is limited in older surfaces.\nThese factors are mentioned in the introduction section (see L92–95) and detailed in the Supplementary Text 2.\n4. Regolith evolution in the CE-5 region\nAs suggested by the reviewer, though regolith thickness and age lack a strong correlation across the Moon, it is feasible to conduct a detailed modeling of the regolith evolution at individual appropriate regions, where the influences of the above factors can be minimized as much as possible. The CE-5 landing region is such a unique location, providing an unprecedented opportunity to estimate the impact flux.\n4.1 The high-resolution imagery acquired by the descent camera enabled a reliable estimate of the regional distribution of regolith thickness. The precise isotopic age of 2.03±0.004 Ga avoids the large dating uncertainties associated with crater counting (addressing factor 3.1).\n4.2 Spectroscopic observations show that the CE-5 landing region of several tens of kilometers across exhibits uniform target characteristics (Extended Data Figs. 1 and 4). Both the frequency of rocky craters and statistics of rock abundance indicate that the P58 unit is representative of typical mare basalts used to construct the NPF (Supplementary Fig. 7). This was explained in detail in our prior response (addressing factor 3.2).\n4.3 As the youngest radiometrically dated mare unit, the CE-5 region has a relatively simple geological evolution history. The source craters contributing to regolith generation at the CE-5 region have been thoroughly documented. By taking these source craters as model inputs, the contribution from large impacts has been explicitly incorporated (addressing factor 3.3).\nBesides, the regolith growth in such a young surface is less influenced by the buffering effect (addressing factor 3.4).\nIn the revision, we emphasized the feasibility of inferring the impact flux from regolith evolution modeling at the CE-5 landing site. See L151–166.\nIn conclusion, we agree with the reviewer that the existing observations do not support a strong correlation between regolith thickness and age. We rephrased our statement regarding the correlation, discussed the potential reasons for the variation of regolith thickness across the Moon, and explained why it is feasible to estimate the impact flux from the regolith thickness at the CE-5 landing region. We explained these issues in the introduction (see L71–65, L120–139, L151–166) and provided more detailed descriptions in Supplementary Texts 1 and 2 and Supplementary Figs. 1–5.\nAnswers to Reviewer 3\nWe appreciate the reviewer’s evaluation on our manuscript. In the revision, we carefully revised our manuscript according to the reviewer’s comments. Specifically, we explained the validity of using a two-layer model and performed additional statistical analyses to support our conclusions. A point-by-point response is provided below.\nComments from the reviewer are in black, and line numbers in the reviewer’s comments refer to those in our original manuscript.\nReplies from the authors are in red, and line numbers in our replies refer to those in the track change revision References are provided at the end of the file.\nMajor comments\nQ1. However, the observed drilling core and seismic data shows a more complicated layer structure of regolith, instead of the underlying assumption that uses two layers model to classify the crater morphology. How the Zhang et al. model tackles this complexity (regolith dynamics and interaction with cratering) is unclear in the current manuscript.\n• We acknowledge that the both the Apollo and Luna core samples revealed a complex layered structure within the regolith, characterized primarily by variations in grain size, composition, maturity, and noble gas concentrations. However, this microscale internal structure does not impact the method we use to estimate regolith thickness, which is based on the mechanical properties of the regolith and the underlying bedrock.\n• We do not agree with the referee’s comment that seismic data indicates a more complicated layer structure of regolith. Actually, the Apollo seismic experiments showed a relatively homogeneous regolith structure: “P-wave velocities measured in regolith at the Apollo landing sites are very low and uniform and range from 92 to 114 m/sec” (Carrier et al. 1991, Lunar Sourcebook, p.516 in chapter 9). The uniform P-wave velocity within the regolith, along with its increase at the base of the regolith layer, suggests a homogeneous regolith overlying bedrock. This seismic characteristic is the basis for estimating regolith thickness (see Nakamura et al., 1975, The Moon, pp. 57–66).\n• The two-layer target model (i.e., a loose granular layer atop a cohesive substrate) has been validated through both impact experiments and numerical simulations. The most widely used method for estimating regolith thickness—small crater morphology—was developed through impact experiments on a regolith-bedrock two-layer target (Quaide and Oberbeck, 1968; Oberbeck and Quaide, 1967, 1968). The four observed morphological types of small lunar craters were successfully reproduced in numerical impact cratering simulations using a two-layer structure (e.g., Senft and Stewart, 2007). In addition, comparisons with regolith thickness measurements from lunar seismic experiments have demonstrated that this method is reliable and largely unaffected by internal layering within the regolith (e.g., 11\nDi et al., 2016; Yang et al., 2021). Therefore, the two-layer model can serve as an effective approximation of the lunar surface (i.e., regolith-bedrock) and has long been employed in estimating regolith thickness. We explained this in L468–472 in the revision.\nQ2. In addition, the crater population in Chang’E 5 site is fitted with Neukum function and seems relatively consistent with isotopically dated ages of mare basalt (2 Ga). In the current manuscript, it is not clear that how the crater population with increased impact flux in the model differs from the observed crater population.\nWe discussed this in L333–346 in our original manuscript.\nThe diameter ranges of impact craters in Yue et al. (2022) and in our study are different. In Yue et al. (2022), craters with diameters >200 m were used to fit the chronology function. In contrast, our study revealed an increased impact flux for sub-hectometer small craters (<100 m), a diameter range to which regolith thickness is most sensitive (see Extended Data Fig. 5). Within this diameter range, the rapid degradation of impact craters produces a bias between the observed (e.g., the Neukum production function) and the actual crater production.\nMinor comments\nQ3. L63–66. This is perhaps the most critical observation from the manuscript. Please consider to add a visual in the main text.\nWe added several sentences describing the growth process of lunar regolith to clarify the correlation between regolith thickness and age. See L64–76.\nQ4. L79. Please add a brief description on this term (first use). If the “self-limiting growth process” does occur, why the regolith thickness data can tell us about the impact flux is a bit odd and not clear in the current manuscript. For example, one expects to see the presence of “self-limiting growth process” at the lunar surface.\nThank you for the suggestion. In the revision, we defined the self-limiting process and the buffering trend in regolith growth (See L68–73). Due to the self-limiting nature of regolith formation, its growth is sensitive to impact flux primarily on geologically young lunar surfaces. This is why we selected the CE-5 landing region for this study, as it is the youngest isotopically dated mare basalt unit.\nQ5. L94–95. Please clarify the sentence and provide a more appropriate reference. At least, reference 29 is the book, so it would be great to know the page numbers that for an interested reader to read it.\nThis sentence described the self-limiting growth process. See our response to Q4. This is described on p.285–286 of Mckay et al. (1991). We provided the chapter information in the revision. See L119.\nThis can be found on p. 198 of chapter 10 in ref. 21. The mean growth rates at the Apollo landing sites can be calculated based on regolith thicknesses and surface ages, which were summarized in Extended Data Table 1. We cited this table in the revision.\nThis represents a 50% difference, which we consider to be significant. Given that the Apollo 11 landing region is 3.8 Ga old, the low regolith growth rate reflects the self-limiting process. This supports our choice of the CE-5 landing site, because the relatively young surface is less affected by the self-limiting process, making it more suitable for assessing the impact flux.\nQ7. L229–231. Please remind of regolith deposition being static. If the reviewer’s understanding is correct, this static regolith growth mechanism cannot dynamically capture the evolution of regolith. If so, please have this discussion for parameter uncertainties.\nThe percentages reported here are derived from Fig. 3d, where the contributions to regolith formation from the observed, invisible, and secondary craters were recorded. As a large number of craters form in LTREM, the bedrock is shattered and excavated, and the regolith is deposited and transported, which is a dynamic process. Different parameters primarily affect the total regolith thickness, while the relative contribution percentages remain largely stable. The 1σ uncertainty is 0.1%. We added this in the revision (L282).\nQ8. L225–257. The modeled regolith thickness using NPF, 2.48–3.32 meter (3 sigma significance) is on the same order of 2.7–4.5 meter though. I suggest to have more statistical analysis using non-biased test on the bootstrapped samples of the 100 realization data (please refer to Huang et al., (2022) supplementary test).\nThanks for the suggestion. In the revision, we performed additional statistical analyses of our results as suggested by the reviewer. We estimated the 95% confidence interval of the difference between the observed and modeled regolith thicknesses from 100,000 bootstrapping samples, which is 0.06–1.35 m. This indicates that the larger regolith thickness is statistically significant (Supplementary Fig. 4). We also estimated the standard error and 95% confidence interval of the mean regolith thickness over 100 realizations, which is less than 0.17% (Supplementary Fig. 6). These analyses show that our conclusion regarding an increased small impactor flux is reliable. See L316–318 and L554–558.\nQ9. L416. The heterogeneity of ray distribution, as well as its ballistic sedimentation, haven’t explored yet, and how its sporadic nature affects SFD is not discussed.\nQ10. L496–497. This sentence sounds another way around; for example, from Minton et al., (2019) the distal ejecta (lateral transport of lunar regolith) drives the change of surface topography. Please consider to revise this sentence and/or to add a reference to support this.\nHere we mean this parameter (diffusivity) controls the rate of regolith transport. This sentence is changed to “Topographic diffusivity controls the rate of lateral transport of lunar regolith”. See L575.\nQ11. L499–500. A side note as “the volume of regolith produced” is first appeared here. How the volume of regolith produced in the manuscript is not stated in the main text and the supplementary. Zhang et al. (2023) has a very detailed description. I suggest to have at least the definition of “the volume of regolith produced” and their components (shattered rock, the volume of regolith changed …) in the main text or the supplementary and refer details to Zhang et al., (2023).\nWe described this briefly and added the reference in the revision. See L584–585.\nQ12. L506–517. This reviewer’s impression is that the LTREM allows us to break down the components of regolith production. With this model capability, the authors can test each component that potentially drive the regolith production in Chang’E 5 landing site, and furthermore the authors can test whether impact flux in <2 Gy is increased. If this chain of thought is right, please have similar sentences in the main text as it was not clear and strong.\nThanks for the suggestion. We emphasized this chain in the main text. See L223–228.\n• Atang, E. F. M., & Bart, G. D. Lunar Regolith Depth Vs. Surface Age: New Measurements in the Lunar Maria. in 55th Annual Lunar and Planetary Science Conference. 1549, (2024).\n• Bart, G. D. Evaluating the Depth–Age Hypothesis for the Evolution of the Lunar Regolith. The Planetary Science Journal, 6(2), 36, (2025).\n• Bierhaus, E. B. et al. Secondary craters and ejecta across the solar system: Populations and effects on impact‐crater–based chronologies. Meteorit. Planet. Sci. 53, 638–671, (2018).\n• Carrier III, W. D. et al. in Lunar Source-Book: A User's Guide to the Moon. Physical properties of the lunar surface. (eds G. H. Heiken, D. T. Vaniman, & B. M. French) Ch. 9, 475–594 Cambridge University Press, (1991).\n• Chertok, M. A., Lucey, P. G., Costello, E. S. & Ireland, S. M. The rock abundance of crater populations as a probe of mare protolith properties. J. Geophys. Res. Planets 128, e2023JE007767, (2023).\n• Costello, E. S., Ghent, R. R. & Lucey, P. G. Secondary impact burial and excavation gardening on the Moon and the depth to ice in permanent shadow. J. Geophys. Res. Planets 126, e2021JE006933, (2021).\n• Di, K., et al. Lunar regolith thickness determination from 3D morphology of small fresh craters. Icarus, 267, 12–23, (2016).\n• Elder, C. M. et al. The subsurface coherent rock content of the moon as revealed by cold‐spot craters. J. Geophys. Res. Planets 124(12), 3373–3384, (2019).\n• Elder, C. M. et al. The variability of lunar mare basalt properties from surface rock abundance. Planet. Sci. J. 4, 244, (2023).\n• Elliott, J. R., Huang, Y.-H., Minton, D. A. & Freed, A. M. The length of lunar crater rays explained using secondary crater scaling. Icarus 312, 231–246, (2018).\n• Fa, W. Z., Liu, T. T., Zhu, M. H. & Haruyama, J. Regolith thickness over Sinus Iridum: Results from morphology and size-frequency distribution of small impact craters. Journal of Geophysical Research-Planets 119, 1914–1935, (2014).\n• Head, J. W. & Wilson, L. Rethinking lunar mare basalt regolith formation: New concepts of lava flow protolith and evolution of regolith thickness and internal structure. Geophys. Res. Lett. 47, e2020GL088334, (2020).\n• Hirabayashi, M. et al. The role of breccia lenses in regolith generation from the formation of small, simple craters: Application to the Apollo 15 landing site. J. Geophys. Res. Planets 123, 527–543, (2018).\n• Hussain, M M., & Mahmud, I. PyMannKendall: A python package for non parametric Mann Kendall family of trend tests. Journal of Open Source Software, 4(39), 1556, (2019).\nhttps://doi.org/10.21105/joss.01556\n• McKay, D. S. et al. in Lunar Source-Book: A User's Guide to the Moon. The Lunar Regolith. (eds G. H. Heiken, D. T. Vaniman, & B. M. French) Ch. 7, 285–356 Cambridge University Press, (1991).\n15\n• Melosh, H. J. Impact Cratering: A Geologic Process. (Oxford University Press, 1989).\n• Melosh, H. J. Planetary surface processes. (Cambridge University Press, 2011).\n• Minton, D. A., Fassett, C. I., Hirabayashi, M., Howl, B. A. & Richardson, J. E. The equilibrium size-frequency distribution of small craters reveals the effects of distal ejecta on lunar landscape morphology. Icarus 326, 63–87, (2019).\n• Nakamura, Y., Dorman, J., Duennebier, F., Lammlein, D. & Latham, G. Shallow lunar structure determined from the passive seismic experiment. The Moon 13, 57–66, (1975).\n• Neukum, G., Ivanov, B. A. & Hartmann, W. K. Cratering records in the inner solar system in relation to the lunar reference system. Space Sci. Rev. 96, 55–86, (2001).\n• Oberbeck, V. R. & Quaide, W. L. Estimated thickness of a fragmental surface layer of Oceanus Procellarum. J. Geophys. Res. 72, 4697–4704, (1967).\n• Oberbeck, V. R. & Quaide, W. L. Genetic implications of lunar regolith thickness variations. Icarus 9, 446-465, (1968).\n• Oberbeck, V. R., Quaide, W. L., Mahan, M. & Paulson, J. Monte Carlo calculations of lunar Regolith thickness distributions. Icarus 19, 87–107, (1973).\n• O'Brien, P., & Byrne, S. Physical and chemical evolution of lunar mare regolith. J. Geophys. Res. Planets 126(2), e2020JE006634, (2021)\n• Quaide, W. L. & Oberbeck, V. R. Thickness determinations of the lunar surface layer from lunar impact craters. J. Geophys. Res. 73, 5247–5270, (1968).\n• Schultz, P. H. & Gault, D. E. Clustered impacts: Experiments and implications. Journal of Geophysical Research: Solid Earth 90, 3701–3732, (1985).\n• Senft, L. E., & Stewart, S. T. Modeling impact cratering in layered surfaces. J. Geophys. Res. Planets 112(E11), (2007).\n• Shkuratov, Y. G. & Bondarenko, N. V. Regolith layer thickness mapping of the Moon by radar and optical data. Icarus 149, 329–338, (2001).\n• Shoemaker, E. M., et al. Observations of the lunar regolith and the Earth from the television camera on Surveyor 7. J. Geophys. Res. Planets, 74(25), 6081–6119, (1969).\n• Speyerer, E. J., Povilaitis, R. Z., Robinson, M. S., Thomas, P. C. & Wagner, R. V. Quantifying crater production and regolith overturn on the Moon with temporal imaging. Nature 538, 215–218, (2016).\n• Stöffler, D. & Ryder, G. Stratigraphy and isotope ages of lunar geologic units: Chronological standard for the inner solar system. Space Sci. Rev. 1, 9–54, (2001).\n• Vanga, S. et al. Rock abundance on the lunar mare on surfaces of different age: Implications for regolith evolution and thickness. Geophys. Res. Lett. 49, e2021GL096710, (2022).\n• Venkatraman, J. et al. Statistical estimates of rock-free lunar regolith thickness from diviner. Planetary and Space Science 229, 105662, (2023).\n• Xie, M., Liu, T. & Xu, A. Ballistic sedimentation of impact crater ejecta: Implications for the provenance of lunar samples and the resurfacing effect of ejecta on the lunar surface. 16\nJ. Geophys. Res. Planets 125, e2019JE006113, (2020).\n• Xie, M., Xiao, Z., Xu, L., Fa, W. & Xu, A. Change in the Earth–Moon impactor population at about 3.5 billion years ago. Nat. Astron. 5, 128–133, (2021).\n• Yang, X., Fa, W. Z., Du, J., Xie, M. G. & Liu, T. T. Effect of topographic degradation on small lunar craters: Implications for regolith thickness estimation. Geophys. Res. Lett. 48, e2021GL095537, (2021).\n• Yue, Z. et al. Updated lunar cratering chronology model with the radiometric age of Chang’e-5 samples. Nat. Astron. 6, 541–545, (2022).\n• Zanetti, M. et al. Evidence for self-secondary cratering of Copernican-age continuous ejecta deposits on the Moon. Icarus 298, 64–77, (2017).\n• Zhang, M., Fa, W. & Eke, V. R. Modeling the evolution of lunar regolith: 1. Formation mechanism through individual simple impact craters. J. Geophys. Res. Planets 128, e2023JE007850, (2023).\n• Zhang, M., Fa, W., Barnard, E. M. & Eke, V. R. Modeling the evolution of lunar regolith: 2. Growth rate and spatial distribution. J. Geophys. Res. Planets 128, e2023JE008035, (2023).\n17"}, {"round": 3, "reviewer_comments": {"reviewer_2": "I appreciate that responses to my reviews have been comprehensive. However, a core problem remains. The manuscript does not present a convincing argument that it is currently feasible to predict flux based on regolith depth, as significant and incorrect assumptions are made regarding the relationship between regolith depth and time, leaving significant uncertainties unaddressed.\nThe timescales examined by statistical methods to correlate age and regolith depth using microwave and crater morphology (>3 Ga) are not relevant to the subtle and recent variations in flux discussed in this manuscript. The data analyzed for the microwave correlation spans ages from approximately 1 Ga to 4 Ga, with an apparent uptick after 3.5 Ga. Similarly, the crater morphology method (Bart et al., 2011) also spans the 2.5 to 3.5 Ga period. I suspect that the statistical correlation results from including regolith depths for surfaces older than 3 Ga, where the impact rate increases significantly. The correlation between regolith depth and age (and recent impact flux) will only be convincing if the data examined are within the relatively flat 1 to 3 Ga timescale.\nConsider performing these statistical analyses of regolith depth versus surface age over the limited timescale of 1-3 Ga to assess whether a definitive depth-age-flux correlation exists for relatively young surfaces (like the CE'5 landing site).\nThe statistical uncertainty in the link between the model and observed regolith depth must also be reflected in the uncertainty in modeled regolith depht and the flux implied by the difference between the model and observations (e.g., Fig. 4). In the response to reviewers, section \"1. Two-sigma confidence interval\", a modeled regolith depth of 2.9±0.14 m is reported; however, the reported error does not account for the uncertainty we find between the observed regolith depth and model predictions versus time at relevant timesclaes. A larger uncertainty in the model-predicted regolith depth will result in significantly different statistics and confidence when compared to the observed depth at the CE'5 landing site. Without an analysis of model uncertainty, it is unclear whether the possibly significant uncertainty invalidates predictions regarding an uptick in flux.", "reviewer_3": "Dear Zhang et al.,\nThis reviewer is pleased by the overall context and revised edits provided by Zhang et al titled “Increased Small Impactor Flux on the Moon as Inferred from Regolith Thickness at the Chang’E-5 Region”. This reviewer would like to thank authors for making their regolith evolution model more understandable and thinking about statistical interpretation from variety sources of the observational constrained regolith thickness dataset. The hypothesis regarding the relationship between the surface age and the regolith growth remain testable, requiring more research work on this topic. Zhang et al., also have provided sufficient content on addressing the effects of primary small impacts on the regolith growth model. As the authors proposed the primary small impacts are dominant source of regolith growth at the Chang’E 5 landing site, that is, this population is currently invisible due to the crater degradation and crater equilibrium, as a result, uncertainties concerning the size frequency distribution (~<85 meter in diameter) are high and have to be constrained from the additional information such as the observed glass spherule ages, boulder distribution, and LROC impact flash observations. The authors have shown a vast room of improving understanding of primary small crater population as their modeling under five types of production functions (in combination of magnitude of increased impact rate and corresponding crater size range) can match the observed regolith thickness. However, this reviewer still feels unclear about the fragmentation of bedrock part in the regolith evolution model though the authors have provided quite some work regarding the working principles of the model. For example, the regolith is defined as breakdown of bedrock and created via either ejecta blanket or breccias lens. This should be fine until the surface is covered by the craters and more subsequent craters in next time step. How does the model handle the surface covered by loosely regolith materials? This is back to the original point from this reviewer regarding the compaction part of the model. Also, this reviewer strongly suggests the authors, if available, have some general discussion on how types of crater morphology reflects the underlying regolith layer structure (the experimental study in Quaide and Oberbeck, 1968). As the backbone of the manuscript is based on the relationship between regolith thickness and crater morphology, it would be great for supplementing this plot if it is available or the analysis is not too challenging. Finally, this reviewer has minor comments for the authors to improve the manuscript context regarding the geological history at Chang’E 5 landing site from Evans et al., (2018) and Fu et al., (2021). Please see the annotated PDF of the author’s revised manuscript (highlighted in purple color with comment box).\nWarm regards,\nYa Huei Huang\n2025/8/22\nReference:\nQuaide and Oberbeck (1968)\nhttps://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/JB073i016p05247\nEvans et al., (2018)\nhttps://doi.org/10.1029/2017JE005421\nFu et al., (2021)\nhttps://doi.org/10.1029/2020JE006797", "reviewer_5": "Review of “Increased Small Impactor 1 Flux on the Moon as Inferred from Regolith Thickness at the Chang’E-5 Region”\nSummary and evaluation\nThis paper makes the case that the regolith thickness in the Chang’E-5 landing region demands a higher impact flux than is found in the standard chronology model. The argument for this is based on a coupled regolith and terrain evolution model. If correct, this enhanced impact flux is important.\nA [very] major strength of this paper is the care with which the regolith evolution modeling has been considered – this is absolutely state-of-the-art and will be a valuable addition to the literature. The main [minor] weak element I’d raise is that this paper is essentially trying to make a fairly significant claim based on this model by comparing it to a kind-of-messy single data point (i.e., the regolith thickness on C’E-5 surface e.g. Fig. 1b). This could lead us astray if either the thickness data has unexpected systematics (i.e., I think the traditional regolith-thickness derivation using morphology were developed on larger craters), or if the model functional behavior or parameterizations are incorrect.\nOverall, however, I think the real methodological innovations of this work outweigh any concerns I have about the applications of the model, so I would argue for acceptance after minor revisions at the author’s discretion. I hope the authors might explore additional applications of their model later and give a few ideas for this in my comments (e.g. #3).\nGeneral comments:\n1) I have been added as a reviewer in a later round. Reading the back and forth in the rebuttal document, I see that the authors ran into an earlier reviewer skepticism regarding the notion that regolith thickens with time and that older units have thicker regolith. I agree with the authors here: it’s a noisy relationship, but one that both (a) is likely to hold from first principles, and (b) looks empirically supported, from a variety of observations. In fact, I think we would be surprised if the relationship was stronger (in agreement with Line 88-95) given (a) the inherent variability of the relevant processes, (b) the relatively narrow age of most mare surface area, and (c) the buffering effect.\n2) I see that, in at least the original iteration of the analysis, the Neukum production function was extrapolated to sizes <10m in diameter (e.g. ExData Fig 6), which is formally below where it was defined. I think it is extremely plausible that the slope of the PF steepens below some size because of the known steeper PF of small grains (Grun 1985). Regardless of the details here (whether steepening happens at 2 m or 20 m), I do not think it alters the results much. Craters in this regime mostly garden (like secondaries) rather than increase regolith thickness, but it is something worth considering.\n3) I think everywhere this paper talks about regolith thickness, it is talking about the median, but there is potentially more of a handle on the relevant parameters if it was possible to invert for the full regolith thickness distribution. I would admit that the available data to do this might not be confident enough, but I would encourage further thinking on that path.\nSpecific comments:\nLine 46: The survival time has some implicit ‘visibility threshold’ built in, see Minton et al 2019; Fassett et al. 2022. Fassett et al 2022 also calculate slightly longer survival times than given here. Note that in looking at Fig. 2, I think the authors’ visibility threshold (ie how degraded a crater has to be to be recognized) is lower than mine.\nLine 63. Little to no regolith is present. The cited Head and Wilson paper disagrees\n(https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2020GL088334). I tend to think the authors are probably right (although regardless the first 50cm-1m of growth isn’t that important because it grows quickly, due to the lack of a significant buffering effect).\nLine 81. I know this paragraph was largely added at the behest of the earlier reviewer, but it feels distracting here. Is there a place later in the paper where it logically could flow better? I think I agree with everything here.\nLine 400. Was this supposed to be a reference to supplementary figure 10-11?\nLine 570-574. I’m not surprised that the diffusivity has negligible effect on the median, though the use of a constant, lower-than-applicable effective diffusivity for big craters (~50-100m+) potentially causes an under-thickening of regolith in their interiors? I think the model is showing this behavior by looking rougher at the sub-saturated sizes in Fig. 3, 1500 Myr, than the real Moon.\nLine 618. Is it transient crater depth, or the whole shape? Curious: Does you model intentionally seek to conserve mass/volume?\nLine 702. It would be nice if the code for LTREM were made available for reproducibility reasons.\nLine 953. I don’t think I understand Fig. 1b. Why is it allowable for this curve to not be monotonically increasing? As I read this, the series for concentric craters would predict less area with thickness <6.5m than <5m, which doesn’t seem logical. For reference, I think this figure is a little like observational data akin to Fig. 6 of Hirabayashi et al. 2018 (though perhaps with axes switched). I would find it useful to have this translated into a best estimate for the regolith cfd/sfd (areal fraction with thickness <X or =X or >X or whatever).\nLine 991. This figure is pretty cool…nice work.\nSuppl. Figure 2, 4. I’ll note that these look sort of like Vanga et al.’s Fig 3 median regolith thicknesses, though probably with some difference in the values."}, "author_response": "Answers to Reviewer 2\nWe appreciate the reviewer’s comments on our manuscript. A point-by-point response is provided below.\nComments from the reviewer are in black, and line numbers in the reviewer’s comments refer to those in our original manuscript.\nReplies from the authors are in red, and line numbers in our replies refer to those in the track change revision. \nQ1. I appreciate that responses to my reviews have been comprehensive. However, a core problem remains. The manuscript does not present a convincing argument that it is currently feasible to predict flux based on regolith depth, as significant and incorrect assumptions are made regarding the relationship between regolith depth and time, leaving significant uncertainties unaddressed.\nThe timescales examined by statistical methods to correlate age and regolith depth using microwave and crater morphology (>3 Ga) are not relevant to the subtle and recent variations in flux discussed in this manuscript. The data analyzed for the microwave correlation spans ages from approximately 1 Ga to 4 Ga, with an apparent uptick after 3.5 Ga. Similarly, the crater morphology method (Bart et al., 2011) also spans the 2.5 to 3.5 Ga period. I suspect that the statistical correlation results from including regolith depths for surfaces older than 3 Ga, where the impact rate increases significantly. The correlation between regolith depth and age (and recent impact flux) will only be convincing if the data examined are within the relatively flat 1 to 3 Ga timescale.\nConsider performing these statistical analyses of regolith depth versus surface age over the limited timescale of 1-3 Ga to assess whether a definitive depth-age-flux correlation exists for relatively young surfaces (like the CE-5 landing site).\nAs noted by Reviewer #5, the correlation between regolith thickness and age holds from first principles, and it is also supported by observations as shown in Supplementary Figs. 2, 4, and 5. We understand the reviewer’s concern that this correlation may be obscured by other factors, especially for surfaces younger than 3 Ga. This has been discussed in the manuscript (L76–78 and Supplementary Text 2), and we agree that the statistical correlation is not particularly strong due to (1) uncertainties in the estimates of age and regolith thickness, (2) variations in target properties, and (3) the stochastic nature of large impact events. However, our goal is not to demonstrate a universal correlation between regolith thickness and age across all the mare units (<3 Ga). It is exactly because of the noise present in the global observations that we chose the CE-5 region for this study. Owing to the unique geological setting, the returned samples, and the availability of high-resolution data at the CE-5 landing site, the uncertainties in both age and regolith thickness, the effects of target property variations, and the spatiotemporal randomness of large impacts, are substantially reduced. Therefore, it is feasible to infer the impact flux from the regolith thickness in the CE-5 region. Please see our prior response for more detail regarding the uniqueness of the CE-5 region.\nQ2. The statistical uncertainty in the link between the model and observed regolith depth must also be reflected in the uncertainty in modeled regolith depth and the flux implied by the difference between the model and observations (e.g., Fig. 4). In the response to reviewers, section \"1. Two-sigma confidence interval\", a modeled regolith depth of 2.9±0.14 m is reported; however, the reported error does not account for the uncertainty we find between the observed regolith depth and model predictions versus time at relevant time scales. A larger uncertainty in the model-predicted regolith depth will result in significantly different statistics and confidence when compared to the observed depth at the CE'5 landing site. Without an analysis of model uncertainty, it is unclear whether the possibly significant uncertainty invalidates predictions regarding an uptick in flux.\nThe uncertainty we find between the observed regolith thickness and predictions from an analytical model versus time represents arises from a series of factors described in Supplementary Text 2 (see also our reply to Q1), which differs from the model/observation uncertainty at the CE-5 region in this study. In such a carefully selected study area, the uncertainties associated with age and regolith thickness measurements are substantially reduced, the target exhibits no anomalous mechanical properties, the buffering effect is relatively weak, and the major source craters have been thoroughly identified and incorporated into our model. Therefore, the overall uncertainty at the CE-5 landing site is relatively small and has been adequately reported (3.6±0.3 m for observation and 2.9±0.14 m for simulation).\nAnswers to Reviewer 3\nWe appreciate the reviewer’s comments on our manuscript. In the revision, we carefully revised our manuscript. A point-by-point response is provided below.\nComments from the reviewer are in black, and line numbers in the reviewer’s comments refer to those in our original manuscript.\nReplies from the authors are in red, and line numbers in our replies refer to those in the track change revision.\nQ1. However, this reviewer still feels unclear about the fragmentation of bedrock part in the regolith evolution model though the authors have provided quite some work regarding the working principles of the model. For example, the regolith is defined as breakdown of bedrock and created via either ejecta blanket or breccias lens. This should be fine until the surface is covered by the craters and more subsequent craters in next time step. How does the model handle the surface covered by loosely regolith materials? This is back to the original point from this reviewer regarding the compaction part of the model.\nThese methodological details have been described in our previous studies (Section 4.1 in Zhang et al., 2023a, JGR-Planets, and Section 2.2 in Zhang et al., 2023b, JGR-Planets). In brief, we estimate the volume of fractured bedrock based on the geometry of the final and transient crater profiles, rather than modeling the physical fracturing process. This approach is applicable to various preimpact target conditions (including regolith-covered surfaces), and any compaction or expansion is inherently incorporated through the crater topographic profile. The treatment adopted here differs from the porosity compaction approach used in the reviewer’s recent study (Huang et al., Nat. Geosci. 15, 531–535, 2022). In the revision, we clarified the principle underlying the simulation of bedrock fragmentation and provided a more precise reference to our previous work (Section 2.2 in Zhang et al., 2023b). See L522–526.\nQ2. Also, this reviewer strongly suggests the authors, if available, have some general discussion on how types of crater morphology reflects the underlying regolith layer structure (the experimental study in Quaide and Oberbeck, 1968). As the backbone of the manuscript is based on the relationship between regolith thickness and crater morphology, it would be great for supplementing this plot if it is available or the analysis is not too challenging.\nThanks for the suggestion. We described the relationship between regolith thickness and crater morphology in the section “Regolith Thickness Estimation from Small Crater Morphology” in Methods. The underlying principle is straightforward: a normal crater forms when the ratio of crater diameter (D) to regolith thickness (T) <4±0.1; a concentric crater forms when D/T >9±0.6; a central mound or flat-bottomed crater forms when 4±0.1< D/T <9±0.6. See the diagram below.\nOn this basis, the percentage of fresh craters with different morphological types (Supplementary Fig. 8) could be converted to the area fraction of regolith thickness (Fig. 1b). For example, in the diameter range D_i, the percentage of concentric craters is 90%, then 90% of the area has a regolith thickness less than D_i/(9±0.6).\nThe crater morphology method has been widely used for regolith thickness estimation in the literature (e.g., Fa et al., 2014, JGR-Planets, 10.1002/2013JE004604; Yang et al., 2021, GRL, 10.1029/2021GL095537; Xie et al., 2021 NA, 10.1038/s41550-020-01241-8; Joulaud et al., 2024, Space Sci Rev, 10.1007/s11214-024-01101-1), and a diagram can be found in Fig. 2 of Zhang et al., 2023a. Therefore, we provided the relevant references and explained the method within the text without adding a diagram in the manuscript. See L426–441. For clarification, we added the reference to Supplementary Fig. 8 in the revision.\nMinor comments\nQ3. L161. May want to be careful with this. Fu et al., (2022) proposed possible non-mare basalt in Chang'E 5 landing site, as well as few quasi-circular mass anomaly from gravity data analysis (Evans et al., (2018). After all, the homogeneous layer of the topmost surface might not be accurate for underlying geological setting. Reference:\nEvans et al., (2018) https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2017JE005421\nFu et al., (2021) https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2020JE006797\nThanks for the references. The proportion of non-mare materials in the CE-5 regolith is only ~5%, as supported by our modeling results (Fig. 4b) and laboratory measurements (Wu et al., 2024, Annual Review of Earth and Planetary Sciences, 52, 159–194), which align with Fu et al. (2021). The contamination of exotic material and the mass anomaly do not influence the basalt mechanical properties in this study. Therefore, we prefer not to change here.\nQ4. L173. “that are widespread on the lunar surface” Consider to delete it.\nWe deleted “that are widespread”. See L146.\nQ5. L233–235. Method section seems lacking of description for fragmentation of bedrocks.\nPlease see our response to Q1.\nQ6. L433–442. Quaide and Oberbeck (1968) was the first to extrapolate the regolith thickness using two layers model, in which is based on the correlation between features of craters versus the apparent diameter of a crater and the regolith thickness. This reviewer has two suggestions:\n1.It might not be straightforward for readers to estimate regolith thickness. Though in the following text the authors explained pros and cons from the types of crater morphology. Please consider to add simple text regarding the regolith thickness estimate. For example, Figure 10 in Quaide and Oberbeck (1968) is an illustrative example.\n2. Please consider to provide your data and plot similar to Figure 8, D_F/D_A vs. D_A/t in Quaide and Oberbeck (1968) study if it is available.\nThanks for the suggestions. There are two approaches for regolith thickness estimation proposed by Quaide & Oberbeck (1968)\n1. Based on the diameter of the crater interior feature and the crater rim-to-rim diameter (D_F/D_A vs. D_A/t, as mentioned in suggestion #2).\n2. Based on the size-frequency distribution of craters with four morphological types.\nPlease note that the reviewer’s suggestion applies to Approach 1. In this study, we employed the Approach 2, which utilizes the relative proportions of crater types and does not require D_F. For an illustrative example, please see the diagram in our response to Q2. The text regarding the regolith thickness estimate can be found at L426–441. We also present our data in a format similar to Fig. 2 in Oberbeck and Quaide, 1967, JGR, 10.1029/JZ072i018p04697 (Supplementary Fig. 8). In the revision, we provided the reference to our previous work for a more detailed description of the workflow. See L432.\nAnswers to Reviewer 5\nWe appreciate the reviewer for the positive review and constructive comments on our manuscript. We carefully revised our manuscript according to the reviewer’s comments.\nParticularly, we discussed the limitation of our conclusion based on a single data point and the future applications of our model. A point-by-point response is provided below.\nComments from the reviewer are in black, and line numbers in the reviewer’s comments refer to those in our original manuscript.\nReplies from the authors are in red, and line numbers in our replies refer to those in the track changes revision.\nGeneral comments:\nQ1. A [very] major strength of this paper is the care with which the regolith evolution modeling has been considered – this is absolutely state-of-the-art and will be a valuable addition to the literature. The main [minor] weak element I’d raise is that this paper is essentially trying to make a fairly significant claim based on this model by comparing it to a kind-of-messy single data point (i.e., the regolith thickness on C’E-5 surface e.g. Fig. 1b). This could lead us astray if either the thickness data has unexpected systematics (i.e., I think the traditional regolith-thickness derivation using morphology were developed on larger craters), or if the model functional behavior or parameterizations are incorrect.\nOverall, however, I think the real methodological innovations of this work outweigh any concerns I have about the applications of the model, so I would argue for acceptance after minor revisions at the author’s discretion. I hope the authors might explore additional applications of their model later and give a few ideas for this in my comments (e.g. #3).\nWe thank the reviewer for the positive comments on the methodological innovations of our work. We agree that our conclusion based on a single data point at the CE-5 site may be limited in generality due to local stochastic variability or potential model biases. In the final paragraph of the revision, we described the limitation in inferring an increased small impactor flux to avoid overinterpretation and discussed avenues for future improvements within the framework of our model (see also our response to Q4). See L382–393.\nQ2. I have been added as a reviewer in a later round. Reading the back and forth in the rebuttal document, I see that the authors ran into an earlier reviewer skepticism regarding the notion that regolith thickens with time and that older units have thicker regolith. I agree with the authors here: it’s a noisy relationship, but one that both (a) is likely to hold from first principles, and (b) looks empirically supported, from a variety of observations. In fact, I think we would be surprised if the relationship was stronger (in agreement with Line 88-95) given (a) the inherent variability of the relevant processes, (b) the relatively narrow age of most mare surface area, and (c) the buffering effect.\nWe totally agree and appreciate the reviewer for the positive and supportive comment.\nQ3. I see that, in at least the original iteration of the analysis, the Neukum production function was extrapolated to sizes <10m in diameter (e.g. ExData Fig 6), which is formally below where it was defined. I think it is extremely plausible that the slope of the PF steepens below some size because of the known steeper PF of small grains (Grun 1985). Regardless of the details here (whether steepening happens at 2 m or 20 m), I do not think it alters the results much. Craters in this regime mostly garden (like secondaries) rather than increase regolith thickness, but it is something worth considering.\nThanks for providing the reference. Based on the impactor population of Grün et al. (1985), Fassett et al. (2022, JGR-Planets, 127, e2022JE007510) estimated crater production rates for diameters <10 m using the crater size scaling law (their Fig. S1, shown below). The PF in the smaller diameter range (<~2.7 m) is indeed steeper than the NPF. As the reviewer noted, craters in this regime primarily contribute to the gardening process and rarely produce new regolith due to the shielding effect (Extended Data Fig. 5). In the revision, we discussed this issue in the Methods section. See L489–494.\nQ4. I think everywhere this paper talks about regolith thickness, it is talking about the median, but there is potentially more of a handle on the relevant parameters if it was possible to invert for the full regolith thickness distribution. I would admit that the available data to do this might not be confident enough, but I would encourage further thinking on that path.\nThanks for the suggestion. We agree with the reviewer that the full regolith thickness distribution carries more information. We provide the spatial distribution of the modeled regolith thickness in Fig. 4a and the probability density functions (PDFs) in Fig. 5a. To better characterize the overall features of the long-tailed distribution pattern and to facilitate comparison with observations, we adopted the median as a representative measure.\nFor future work, the regolith thickness corresponding to the peaks of the PDFs (mode), as well as the quartiles and the interquartile range (IQR) of the thickness, could serve as useful statistical measures to more comprehensively characterize the regolith evolution process. For example, for the same median thickness, a large IQR indicates a highly heterogeneous regolith thickness distribution. This may be related to the slope of the crater production function: a steeper PF, dominated by small craters, may lead to a more uniform regolith distribution and thus a smaller IQR, whereas a shallower PF with more large craters may produce a wider thickness spread. This may provide a worthwhile path for future exploration, especially as more accurate observations of regolith thickness become available to support such analysis. In the revision, we briefly discussed the potential improvement to offer insights for future studies. See L388–393.\nSpecific comments:\nQ5. Line 46: The survival time has some implicit ‘visibility threshold’ built in, see Minton et al 2019; Fassett et al. 2022. Fassett et al 2022 also calculate slightly longer survival times than given here. Note that in looking at Fig. 2, I think the authors’ visibility threshold (ie how degraded a crater has to be to be recognized) is lower than mine.\nYes, the survival time here is estimated based on Basilevsky (1976), which did not provide an explicit visibility threshold (a value <1/12). Considering the low-resolution orbital imagery in the 1970s, it is likely higher than 0.05 used in Minton et al. (2019) or 0.04 in Fassett et al. (2022). For Fig. 2, we checked the data and found that the visibility threshold is indeed lower than 0.05, since we want to count as many craters as possible that contribute to the regolith at the landing site.\nQ6. Line 63. Little to no regolith is present. The cited Head and Wilson paper disagrees\n(https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2020GL088334). I tend to think the authors are probably right (although regardless the first 50cm-1m of growth isn’t that important because it grows quickly, due to the lack of a significant buffering effect).\nWe agree with the reviewer. We cited Head and Wilson paper mainly because of their introduction of the concepts of impact flux trend and buffering trend, which are relevant for our modeling framework. In the revision, we cited McKay et al. 1991, Lunar Sourcebook here.\nQ7. Line 81. I know this paragraph was largely added at the behest of the earlier reviewer, but it feels distracting here. Is there a place later in the paper where it logically could flow better? I think I agree with everything here.\nWe agree with the reviewer. In the revision, we simplified this text and moved it forward to improve the overall logic. See L73–78. We also added a brief summary at the end of Supplementary Text 2.\nQ8. Line 400. Was this supposed to be a reference to supplementary figure 10-11?\nYes. We added the reference in the revision. See L367.\nQ9. Line 570-574. I’m not surprised that the diffusivity has negligible effect on the median, though the use of a constant, lower-than-applicable effective diffusivity for big craters (~50-100m+) potentially causes an under-thickening of regolith in their interiors? I think the model is showing this behavior by looking rougher at the sub-saturated sizes in Fig. 3, 1500 Myr, than the real Moon.\nYes, this is an insightful comment. Using a uniform equivalent diffusivity across the entire simulation area indeed leads to an underestimation of the degradation rate for large craters. This could make a difference to the spatial distribution pattern of regolith thickness. However, since degradation mainly redistributes regolith rather than producing new regolith, it does not affect the median regolith thickness in any significant manner, as the reviewer noted.\nIn future work, we plan to update the model by allowing the degradation rate to vary across regions according to local slope to better simulate the degradation of craters with different sizes (larger craters generally contain more steeply sloped areas). We described the limitation of the equivalent diffusivity in the revision. See L550–556.\nQ10. Line 618. Is it transient crater depth, or the whole shape? Curious: Does you model intentionally seek to conserve mass/volume?\nIt is the whole shape of the transient crater that influences regolith thickness. However, since a bowl-shaped cavity is commonly used to represent transient crater geometry (e.g., Melosh 1989; Hirabayashi et al., 2018, JGR-Planets), the crater depth becomes the key parameter. In this study, we varied the transient crater shape by assigning different transient crater depths (L537–538). For accuracy, we refer to the transient crater depth here, though both descriptions are correct.\nThe volume change during impact cratering process is considered based on the observed topographic profiles of 99 fresh lunar craters (extracted by Cai and Fa 2020, doi.org/10.1029/2020JE006429; Yang et al. 2021, doi.org/10.1029/2021GL095537). As shown in Fig. 6d in Zhang et al. 2023a (JGR-Planets, 128, e2023JE007850), kilometer-sized craters exhibit a negative volume change, while the volume for craters smaller than ~100 m is nearly conserved. The volume-neutral nature of sub-hectometer craters ensures that numerous small craters that do not penetrate the existing regolith layer neither add nor remove regolith in an area. We clarified this in the revision. See L514–517.\nQ11. Line 702. It would be nice if the code for LTREM were made available for reproducibility reasons.\nWe plan to release the model as software in the future. It takes some time because of user interaction design, docstrings, and code standardization. At present, interested readers may obtain the code from the corresponding author upon request.\nQ12. Line 953. I don’t think I understand Fig. 1b. Why is it allowable for this curve to not be monotonically increasing? As I read this, the series for concentric craters would predict less area with thickness <6.5m than <5m, which doesn’t seem logical. For reference, I think this figure is a little like observational data akin to Fig. 6 of Hirabayashi et al. 2018 (though perhaps with axes switched). I would find it useful to have this translated into a best estimate for the regolith cfd/sfd (areal fraction with thickness <X or =X or >X or whatever).\nThis is a good question. In principle, the cumulative distribution of regolith thickness should be monotonic. If the spatial distribution of regolith thickness were known, the calculated cumulative distribution would necessarily be monotonic. For example, the figure below shows the monotonic cumulative distributions of regolith thickness derived from the models of Oberbeck et al. (1973), Hirabayashi et al. (2018), and Zhang et al. (2023b).\nHowever, the cumulative regolith thickness inferred from crater morphology is obtained based on the proportion of craters with different morphologies within each diameter bin. It can be regarded as a random sampling of the underlying distribution, and therefore exhibits fluctuations depending on the number of crater samples. Ideally, with a sufficient number of fresh craters, the percentage of normal craters should decrease with diameter, whereas that of concentric craters should increase, resulting in a monotonic regolith thickness distribution (as an example, see Figs. 9 and 11 in Fa et al., 2014, JGR-Planets, 119, 1914–1935). In practice, however, the limited number of observed fresh craters leads to statistical fluctuations. For example, in a small diameter range D to D , if 90% of craters are concentric, the regolith thickness in 90% of the area is estimated to be <D/(9±0.6). In the next larger bin D to D , there may be only 80% concentric craters, implying only 80% of the area having regolith thickness <D /(9±0.6), which appears somewhat counterintuitive.\nSimilar behavior has been reported in Fig. 10 of Oberbeck & Quaide, (1968), Icarus, 9(1–3), 446–465, and Extended Data Fig. 4 in Xie et al., (2021), Nature Astronomy, 5(2), 128–133.\nHere Fig. 1b is actually a cumulative frequency distribution that includes sampling uncertainty, which explains why the areal fraction with thickness <X is not strictly monotonic. Considering this, we adopt the result calculated from normal craters as the best estimate of median regolith thickness due to its smaller uncertainty. In the revision, we briefly explained the non-monotonic behavior. See L437–439.\nQ13. Line 991. This figure is pretty cool…nice work.\nThank you.\nQ14. Suppl. Figure 2, 4. I’ll note that these look sort of like Vanga et al.’s Fig 3 median regolith thicknesses, though probably with some difference in the values.\nYes. The median regolith thickness in Vange’s Fig. 3 was calculated based on the regolith growth model of Hirabayashi et al. 2018 JGR-Planets 123, 527–543. In our previous study, we compared the results of our analytical regolith growth model and Hirabayashi et al. (2018), and the general trends are similar (Fig. 10 in Zhang et al. 2023a JGR-Planets, 128, e2023JE007850)."}]

Decision letter