Formation of Small Craters in the Lunar Regolith: How Do They Influence the Preservation of Ancient Melt at the Surface?

The impact melt that records the formation time of basins is essential for the understanding of the lunar bombardment history. To better understand melt distribution on the Moon, this study investigates mixing of melt by small impacts using a Monte Carlo numerical model. The obtained mixing behavior is then integrated into a larger scale model developed in previous work. While large impacts produce most of the melt volume in both the regolith and megaregolith, we find that the dominant source of melt near the surface is small impacts. Material in the top meter is affected mainly by impacts that form craters <5 km in diameter. In the uppermost 10 cm, melt with age <0.5 Ga is abundant; while as depth increases older melt is increasingly present. This may indicate that the excess of impact melt <0.5 Ga in lunar samples from the near surface is caused by the cumulative mixing of small impacts. A comparison of the age distribution of melt derived from craters of different sizes with that of impact glass constrains the size of spherule‐forming impacts. Our model is consistent with observations if most impact glass spherules from the near surface are produced by <100 m craters and >100 m craters do not contribute abundant spherules. The distribution of the datable melt with depth is also analyzed, which is essential for future sampling missions. Excavated materials of young and large craters (>100 m on highlands; >10 km on maria) appear to be the most fruitful targets.


Model
The melt diffusion in megaregolith caused by large impacts was investigated by Liu et al. (2020) based on their developed numerical model. In their model, with the knowledge of lunar impact flux, the formation sequence of differently sized craters was calculated. Impact events occurred in chronological sequence. For each impact event, the process of excavation, ejecta deposition, and melting were all simulated based on the scaling laws. After the formation of each impact crater, its influence on the material composition both over the surface and with depth was traced. In their model, the global lunar surface was considered, and the material composition was recorded in the uniformly distributed point sets on a sphere with the radius of the Moon.
This study is focused on the regolith mixing by small impacts, with the aim to integrate the net behavior of the process into the model of Liu et al. (2020) to improve our understanding of the presence of melt at the surface. In Liu et al. (2020), only the formation of craters larger than 5 km in diameter was simulated. To complement these simulations, this study focuses on the cratering events below this size. The newly designed Monte Carlo model runs in a slightly different way due to the considered regional area (see Section 2.1). The simulation results will be combined with those from the larger-scale simulation (Liu et al., 2020). To estimate the influence of small cratering events on the melt distribution in the megaregolith, melt products of small craters are taken to be the portion of megaregolith that is remelted by small impacts, and the initial melt composition of the megaregolith is derived from Liu et al. (2020).
Although the interpretation of impact flux of the Moon over some periods remain controversial (e.g., at ∼3.9 Ga, Boehnke & Harrison, 2016;Neukum et al., 2001;Ryder, 1990;Turner et al., 1977;Wetherill, 1975;recent 0.5 Ga, Culler et al., 2000;Mazrouei et al., 2015;McEwen et al., 1997;Robbins, 2014;Vokrouhlický et al., 2017;Hartmann, 2019, and references therein), it is generally accepted that, in the early period of the lunar bombardment history, the impact flux was relatively high and a great number of craters larger than hundreds of kilometers were formed; while in the later period, the impact flux was much smaller, both in amounts and size of the impactors, resulting in fewer large craters (Fassett et al., 2012;Neukum, 1983). The thick ejecta of large impact events, particularly of giant basins, buried the product of smaller impact events to a great depth, erasing their signature on the surface (Liu et al., 2019(Liu et al., , 2020. The occurrence of small impact events during the early period is, therefore, not involved in this study. The production and chronology functions of Neukum (1983) are applied to calculate the impact rate (Figure 2a), based on the occurrence time and the crater size of impact events using the Monte Carlo method. Due to the lower limit of the production function, the minimum crater diameter (D min ) considered is 10 m.
While the small craters are forming, impactors that form >5-km craters are not considered. However, the obtained diffusion features of small craters in this study should still be representative of most of the present-day lunar surface, since most of the near-surface has been merely gardened by small craters due to the low frequency of large craters during the late period. Based on the Neukum impact flux, it can be calculated that in the recent 3.5 Ga, ∼1,500 craters >5 km would have occurred on the Moon. Taking the disturbance region of each impact to be the areas where the majority of ejected materials are deposited (2 radii from the crater center), it can be calculated that only 0.7% (Ga −1 ) of the lunar surface is affected by large craters. This indicates that most of the lunar surface should preserve the trace of impact mixing by small craters since 3.5 Ga ago. The start time of this model is thus taken to be 3.5 Ga. For the even more recent period when the impact flux is almost constant (<3 Ga), there would be ∼800 craters >5 km and their affected region is even smaller (0.4% Ga −1 ). To understand the melt distribution over some regional areas that were distinctly affected by the larger crater during the late period, the additional simulations where the thick ejecta cover the study region at various times during the formation of small craters have been performed (Section 3.2).

Simulation Area
Due to the great number of small craters, the statistical distribution of their melts on a region can be representative of the consequence of the cumulative gardening process of the whole lunar surface. In addition, the simulation of a small region, rather than a global surface, reduces the number of impact events, increasing the model efficiency. But the area should be large enough for the occurrence of craters 1-5 km in diameter. Taking the occurrence of 1-km craters as an example, in the recent 3.5 Ga, there is a low chance they are present over an area smaller than ∼10 km 2 . However, given the frequent occurrence of craters with the considered D min (10 m), the number of impact events that need to be simulated significantly increases with the area and hence would strongly impair the model efficiency. Overlapping areas with different sizes are used to reach a compromise between the occurrence of the craters >1 km and the limit of the number of smaller craters that need to be simulated. Three square areas of different sizes ( Figure 1) are modeled: region A of 500 × 500 km 2 ; region B of 2 × 2 km 2 ; and region C of 1 × 1 km 2 . Region C is the target area and D min is 10 m. Few craters >1 km would occur in 3.5 Ga over this small area. To account for the influence of craters >1 km on the target region C, it is embedded in a larger region A where D min is taken to be 1 km. Due to the difference of D min between the regions A and C, some materials that are moved outside the target region C will not be transported back to the target region C (compared to if a D min of 10 m would be used for A as well). The transition region B where D min is the same as in target region C is thus designed as well. With this approach, the mixing process of the predominant small craters is simulated without neglecting the influence of larger craters in the vicinity.
To investigate the effect of the size of transient region B on the melt abundance in the target region C, the simulation with a larger region B was run. The results show that when the area of region B is 25 times larger than the above setting, due to the formation of more craters, the total volume of the generated melt is greater, and hence more melt can be transported to the target region C. However, it is found that the relative abundance of melt with different ages appears to be similar. In addition, it shows that the bulk fraction of melt in the top meter over the 4 km 2 region B and that over the 25 km 2 region B is almost identical with a difference <0.5%. This indicates that the average melt abundance over an area larger than 4 km is representative. If the volume of melt is of interest, it can be calculated by combining the bulk fraction and the relative abundance. To improve the model efficiency, bulk fraction and relative abundance of the smaller transient region B is used in the model. Focusing on the mixing behavior of small craters, in the following discussions, the relative melt abundance over the target region is presented where the relative melt fraction of the zone with the most abundant melt is equal to 1.0.
The study region is divided into cells to record the material composition ( Figure 1b). The target region C has the highest resolution where the cell is 5 × 5 m in size. Regions far away from the target region have relatively coarse resolutions: the cells over the region B (excluding the contained target region C) are 50 × 50 m in size; the cells over the further regions (excluding region B and target region C) is 500 × 500 m in size. The increasing resolutions with the decreasing D min and the increasing distance from the target region help to better preserve the composition of the melt generated by the very small craters in the target region.

Melt Production
Given a crater with a final diameter of D, its transient crater diameter (D t ) is 0.8D (Melosh, 1989), and the excavation depth, d exc , is approximately 0.1D t . The excavation volume (V exc ) is estimated to be 1/3 of a disk with d exc in thickness and D t in diameter. The excavated material is ejected outside the crater cavity and is deposited near the crater rim as ejecta layers. Since the majority (∼75%, Liu et al., 2019Liu et al., , 2020 of the ejected material is distributed within five radii from the crater center, only the ejecta in this range is traced and all the ejected material is taken to be displaced in this range to conserve the mass. The thickness of ejecta layers, δ, decreases with distance from crater center, r (Pike, 1974) where A (km 4 ) varies for differently sized craters to conserve mass and is calculated by taking the integrated volume within five radii exactly to be V exc .
The scaling law of Cintala and Grieve (1998) that directly relates D t with the melt production was applied to calculate the total melt volume (V melt = 1.4e −4 D t

3.85
) while investigating the impact mixing within the megaregolith (Liu et al., 2020). However, the impact energy of the small craters, which mainly relates to the impact velocity and the projectile size, may not be great enough to melt the target material. Other scaling laws that take both factors into account, therefore, are used to calculate the melt production in this study (Abramov et al., 2012). Given D t , the diameter of the projectile, D p , is where  p is projectile density (dunite, 3,320 kg/m 3 );  t is target density (anorthosite, 2,940 kg/m 3 ); g is the lunar gravity (1.61 m/s 2 ). The impact velocity (v i ) is taken to be the most common value of 18 km/s. In principle, such impacts possess enough energy to melt the target material. Then, the volume of impact melt is calculated LIU ET AL.
10.1029/2020JE006708 4 of 17 where E m is the specific internal energy of the target at shock pressure that results in melting upon release (highland anorthosite, 3.42 × 10 6 J/kg);  is impact angle and is taken to be the most common value of 45°. The majority of the generated impact melt remains inside the crater cavity, and about 25% the melt is entrained in the ejecta and distributed outside the crater. The thickness of melt, δ m , is assumed to decrease with r as an inverse square relationship, based on the results of unpublished impact simulations (pers. comm. Kai Wünnemann) Note, melt in the ejecta does not actually occur as a distinct layer, but is distributed throughout the ejecta layer.
      m / r r represents the melt content in the ejecta layer at a given distance r. Similar to A, A m (km 3 ) also varies for differently sized craters and is estimated by taking the integrated volume within five radii exactly to be the melt volume that is ejected from the cavity (i.e., 0.25V melt ). The interaction between Equations 1 and 4 implies that, according to our model, beyond a certain radial distance the emplaced ejecta are completely composed of melt. However, this is inconsistent with terrestrial observations French, 1998;Hofmann & Gnos, 2006;Stöffler et al., 2013) and suggests that either or both (1) and (4) are not valid at large distances from the crater. This is another reason that we only apply these scaling relationships to the ejecta within five crater radii of the crater rim. Impact events occur one by one starting with the oldest. For each impact event, we first calculate its gardening region (i.e., all the cells within 5 radii). The thickness of emplaced ejecta and melt in each cell is then LIU ET AL.   (Neukum, 1983). (b) Cumulative distribution of craters over the regions A and B. (c) Distance from the center of the target region of (c) >100-m craters and (d) >1-km craters, where the bigger points indicate larger crater size. In (c) and (d), the craters whose ejecta reach the center of the target region while their formations are highlighted in red. As can be seen from Figure 3a Crater 1 (120 m), Crater 2 (740 m), and Crater 3 (300 m) in (c) contribute a significant amount of melt to the target region C; Crater 4 2 km in diameter in (d) leaves a detectable trace on the target region. calculated using the scaling relationships (1)-(4) and the distance between the center of the crater and the cell. As each crater forms, the number of ejecta layers grows. To avoid the need to track the ejecta layer of each crater we therefore periodically redefine the layer sequence in each cell. The depth of each layer in the sequence is defined as: 10 −3 × 1.9 a m, where a = 0, 1, …, 20. In this way, shallower layers are thinner preserving a good resolution of the melt component in the near-surface. Each time the layer sequence is redefined, we assume that the material within each new layer is a uniform mixture of the individual ejecta layers in that depth range.

Crater Distribution
Over the central target region C and the transient region B, craters larger than 10 m are simulated. The distribution of these craters is shown in Figure 2b. There are 87,889 craters in total over the past 3.5 Ga. The maximum crater diameter is 740 m; 90% of the craters are smaller than 20 m. In region A, 1,200 craters occur. The maximum crater is 4.9 km in diameter; half of the craters are larger than 1.2 km in diameter.
Since large craters that are close to the target region would significantly alter the melt distribution of the target region, their relative distance to the target region is also presented. Figure 2c displays the distance of craters >100 m in diameter to the center of the target region: there are 82 craters in total; 46% of these craters are >1 km away from the target region. Figure 2d displays the distance of craters >1 km in diameter to the center of the target region: 90% of those craters are >100 km away from the target region; only one crater (Crater 4) is located within 10 km.

Melt Distribution Resulting From Small Craters
Each impact event produces melt with the age of its occurrence time. With the long-term impact gardening, differently aged impact melt would present in the near-surface. The calculated melt distribution over the target region C (Figure 3a) shows that melts are most abundant in the top 10 cm. At the lower depth (10 cm-∼1 m), very little melt younger than ∼500 Ma is present. Some melts are present at even lower depth (>1 m), but the concentration is small. Below 10 m almost no melt can be found. Even when the impact flux was relatively strong (>3 Ga, Figure 2a), the melt generated from small cratering events still accumulates in the top surface. The concentration of impact melt from small craters within the first meter of regolith indicates that such small events merely affect the melt composition of shallow material. Thus, small impact events, despite their frequent occurrence, do not significantly change the distribution of the melt generated by large craters that occur lower in the stratigraphy column. In other words, the melt of the large craters remains in place until the present day if there is no disturbance by other late-forming large craters.
As can be seen from Figures 2c and 2d, many large craters occur near the central target region. The ejecta of 27 craters >100 m in diameter reach the central target region when they form (highlighted in red in Figure 2c). These cratering events disperse a significant amount of melt over the target region, such as the 120-m Crater 1, 740-m Crater 2, and 300-m Crater 3 that occurred 0.4, 1.1, and 0.3 km away from the target regions at 0.4, 0.6, and 1.1 Ga, respectively. The ejecta of only a single crater >1 km (2-km Crater 4 at 1.3 Ga highlighted in red in Figure 2d) affected the central target region. However, despite being larger, due to the greater distance to the target region the deposited melt volume in the target region is less than that of Craters 1, 2, and 3. Without the disturbance of the subsequent large craters, the majority of the deposited ejecta remain around the crater rim, but some could be present further away, as a result of transport by subsequent impact gardening ( Figure S2 in Liu et al. (2020)). Therefore, the melt of the other >100-m craters could also be delivered to the target region by the subsequent impact events. However, as can be seen from Figure 3a, their volumes are too small to be detected.
The distribution of the differently aged melt varies with depth. Impact melts of small craters are abundantly deposited in the very top surface (about the top 10 cm). Except for the large volumes of melt generated by the relatively large impact cratering events (e.g., Crater 1, Crater 2, and Crater 3), at the surface, the abundance of the melt generally decreases with time, so that melt younger than ∼0.5 Ga accounts for about half of the total melt products. Deeper in the regolith (∼10 cm-1 m), the amount of the melt is lower and the melt of different ages is similarly abundant (Figure 3a). Hardly any melt generated by small impact events can be found deeper than ∼1 m, and only some melt from larger craters (hundreds of meters in diameter) is present.
No crater larger than 5 km is considered in this model. If such a crater occurs near the target region, its thick ejecta would likely change the melt distribution shown in Figure 3a. To investigate such an influence, simulations with layers of various thicknesses and ages being deposited in between the <5-km cratering events were performed. This thick layer is supposed to simulate the ejecta of large craters. These ejecta generally have a higher velocity. During emplacement mixing with the local material in the target region could, therefore, occur. The degree of local mixing is related to the distance to the crater center. If the target region is closer to the crater less of the local materials will be disturbed (Oberbeck et al., 1975). The weakest and the most intense local mixing occur close to and far away from the crater rim, respectively. Two extreme cases with no local mixing (Figures 3b1-3e1) and thorough mixing (Figures 3b2-3e2) are thus considered. Due to the constant subsequent impactor flux, the melt distribution of small craters after the deposition of thick ejecta is always identical to the scenario in Figure 3a. The distribution of the melt generated earlier than the emplacement of thick ejecta is related to the thickness and the age of the ejecta and the degree of the local mixing. As can be seen from Figures 3b1 and 3c1, when the local mixing is negligible and the deposit of the ejecta is relatively thin (Figures 3b1 and 3c1), some of the buried melt is reexcavated to the surface by subsequent impact events, and the earlier the coverage occurs the more melt can be reexcavated to the surface. When the deposit of the ejecta is thicker (Figures 3d1 and 3e1), the buried melt becomes difficult to deliver to the surface and there is almost no melt older than the ejecta coverage presented in the top surface. However, if the materials of the target region are thoroughly mixed with the foreign ejecta, the early formed melt is mixed into the high-volume foreign ejecta diluting its concentration. As a result, the LIU ET AL.
10.1029/2020JE006708 7 of 17 To present the distribution features, the relative melt fraction is calculated where the relative melt fraction of the zone with the most abundant melt is equal to 1.0 (Section 2.1). (b1-e1) and (b2-e2) display the same scenarios as (a) but with the interruption of ejecta of different thicknesses (1 and 10 m) at various times (0.5 and 2.0 Ga indicated by the vertical lines). Different extreme case scenarios are shown, in which the emplacement of the ejecta causes no local mixing and thorough local mixing. The strength of the color indicates how much of that material would be brought to the surface. For the extreme case of "no local mixing," when the ejecta coverage is thin, the abundance of the buried old melt that could be found at the surface is directly related to the occurrence time of the ejecta emplacement b1, c1; when the emplaced ejecta is rather thick, little buried old melt could be found at the surface (d1, e1). For the extreme case of "thorough local mixing," the great volume of the emplaced ejecta significantly dilutes the local material leading to a rather small fraction of the affected melt product. In such a case, the abundance of the affected old melt at the surface would be small and has no clear relation to both the occurrence time and the thickness of the emplaced ejecta.
abundance of the early formed melt after the emplacement of the ejecta is too small to be observed at the surface (Figures 3b2-3e2). Both extreme cases of mixing do not erase any information on materials but alter the spatial distribution of early generated melt products by either making them deeply buried or diluting their concentration in the near-surface. Due to the various degrees of local mixing at the different distances from large craters, the material of the target region in the greater depth might be shielded from such process leading to the surface material with the composition in between the presented two extreme estimates. Some melt of small craters could be entrained in the ejecta of large craters and add to the target surface. However, given the great excavation volume for craters >10s km and the small amount of melt derived from small impacts, the melt of small craters is itself significantly diluted while the excavating processes. Taking a 10km crater as an example, its excavation depth is about 1 km. Assuming all the components are thoroughly mixed during excavating, the fraction of the top surface (top 1 m, the seating depth of most small craters' melt) materials over the total excavated materials is ∼0.001 (1 m/1 km). Due to the low concentration, the added melt of small craters would not significantly change the melt distribution features shown in Figure 3.

How Does This Process Influence the Preservation of Ancient Melt Products Near the Lunar Surface?
Liu et al. (2020) calculated the melt abundance in the megaregolith caused by the cumulative impact mixing of craters >5 km. To better constrain the distribution of the melt close to the surface, the influence of impact gardening by smaller craters on the preservation of the melt in megaregolith needs to be considered. In this study, melt products of small craters are taken to be the portion of megaregolith that is remelted by small impact. The distribution of melt obtained in this study results from the cumulative impact mixing of small craters without the influence of craters >5 km. In an area where a large crater formed recently, the abundance of earlier-formed melts will be lower, but the abundance of the later produced melts of small craters should be comparable to an area where no larger crater formed (Figures 3a-3e).
In the following, we will focus on the melt distribution in the top meter where the material is extensively mixed by small impacts. Simulations of regions with different sizes (see Section 2.1) show that the bulk fraction of <3.5-Ga melt in the top meter is 3.5 × 10 −3 km −3 . Based on the relative distribution of melts (Figure 3a) and the estimate of the bulk fraction, the volume of melt over the global surface is calculated (black histograms in Figure 4b). As can be seen from Figures 3a and 4b, the small craters only garden the material close to the surface and the fraction of produced melt is small. The distribution of melt from craters >5 km in diameter in both the top of the surface and the underlying layers therefore almost remains unchanged and should be consistent with the predicted melt abundance in the megaregolith of Liu et al. (2020).  period is heavily affected by the formation of the large-sized craters, creating notable peaks (e.g., at 0.3 and 1.7 Ga; gray histograms in Figure 4b). There is an increase in melt products created by small craters in the last 0.5 Ga, but not in those created by larger impacts. Between 3 and 3.8 Ga, the impact flux is increasing but no basin-forming events occurred (the youngest basin, Orientale basin, by current estimates formed at ∼3.8 Ga; Fassett et al., 2012). The relative abundance of melt is higher for older ages, which is a consequence of the increasing impact rate going back in time. During the basin-forming period (>3.8 Ga), the impact rate was even higher. Unlike in the 3-3.8 Ga period, the abundance of melt formed in this period (>3.8 Ga) in the present top surface layer is relatively constant. This is likely due to wide-spread ejecta of large-scaled late-forming basins significantly altering or erasing the traces of older basins. For example, older melts that were excavated during the formation of the Imbrium basin were extensively mixed and entrained in the Imbrium ejecta dispersing them over large areas onto the surface. It is possible that all Apollo and Luna landing sites are dominated by Imbrium ejecta (e.g., Haskin et al., 1998). The melt products of the basins that were slightly older than Imbrium and were far away enough to have not been significantly affected by the Imbrium ejecta would be in great abundance on the surface and would contribute to the modeled peaks observed in Figure 4b. For example, given the formation sequence of basins in the model of Liu et al. (2020), the peaks at ∼4.0 and ∼4.2 Ga are derived from the occurrence of Crisium and Serenitatis, respectively.
In the period with the higher impact rate, the higher frequency and the larger size of impactors lead to the volume of the generated melt being orders of magnitude more abundant. The simulation results suggest that the small craters will only melt a small amount of megaregolith, although they occur in large numbers. As can be seen from Figure 4b, even in the top meter the melt generated by larger craters and basins is the predominant impact melt component. In the recent 3 Ga, the total volume of melt generated by the craters >5 km is about 10 times larger than that generated by the smaller craters. That is, on average, <10% of the megaregolith material that is exposed on the lunar surface has been remelted by small craters. Therefore, even on the surface, melt generated by large cratering events dominates.
Knowledge about the distribution of basin melt is important as it is crucial for calibrating the lunar chronology function (Hiesinger et al., 2011;Neukum & Ivanov, 1994;Neukum et al., 2001;Stöffler & Ryder, 2001;Wilhelms, 1987). Liu et al. (2020) predicted that the majority of the basin melt, especially from the younger basins, should be well-preserved in the near-surface. After the basin-forming epoch, only large impact events could significantly alter the local impact melt distribution. Due to the small amounts of basin melt that is remelted by small impacts, we expect that the distribution of basin melt should be basically identical to that shown in Liu et al. (2020).
The seating depth of the melt with different ages is presented in Figure 5. Most of the melt from small impacts is concentrated in the top layer whereas most of the melt generated by large impacts is concentrated at greater depth. In the megaregolith, the majority of the melt generated after the basin-forming period is deposited at ∼10 m ( Figure 5a); most of the melt formed during the basin-forming period was buried at LIU ET AL.
10.1029/2020JE006708 9 of 17 In (a), the red points represent the melt generated by craters <5 km. In both (a) and (b), the black points represent the melt generated by craters >5 km (Liu et al., 2020). The points indicate the depth to where 50% of the melt is seated; the upper and the lower boundaries of the gray bars are the depths to where 25%, and 75% of the melt is seated, respectively; the upper and the lower caps indicate the depths where 10% and 90% of the melt is seated, respectively. When the upper or lower caps are beyond the range of y axis, they are indicated by the arrows.
∼100 m and deeper (Figure 5b). Still, some basin melt (<10% of the total basin melt) can be found close to the surface, volumetrically eclipsing the melt from small craters (Figure 4b).

Impact Flux Over the Recent 3 Ga
Impact melt products (melt rocks and glasses) provide an important record for constraining the lunar impact flux. Returned impact melt rocks generally have ages older than 3 Ga, and therefore rarely provide insights into the impact flux during the past 3 Ga. Most of the well-constrained ages of impact melt glasses fall into the recent 3 Ga (Zellner, 2019), and thus provide a record of the impact flux during this time. Moreover, impact glass spherules are ubiquitously distributed Melosh & Vickery, 1991;Reid et al., 1972), and a large number of impact glasses have been found in the lunar regolith soil samples collected from the Apollo landing sites and have been extensively dated. Previous studies found an apparent excess of ages <0.5 Ga of impact glass spherules (Culler et al., 2000;Levine et al., 2005;Zellner & Delano, 2015). As an example, Figure 6b shows the relative abundance of the 40 Ar-39 Ar ages of 128 lunar impact glass (70 spheres and 58 shards) assessed by Zellner (2019). One explanation for the excess of impact glass in this period is an increase in the impact flux during the late Copernican (Culler, 2000). However, it may also be caused by a preservation bias (Zellner & Delano, 2015) or by the sampling strategy since all the samples were collected from the shallow surface (Huang et al., 2018;McKay et al., 1991).
As the dated lunar impact glass spheres and shards come from lunar regolith samples, whose typical sampling depth is ∼10 cm (McKay et al., 1991), their age data are best compared to the modeling results of the upper 10 cm. In our model, we find that an increased amount of impact melts younger than 500 Ma will be present in the upper 10 cm (Figure 6a) without requiring an increased impactor flux during this time. The observed excess of young lunar glass ages is, therefore, likely due to an inherent bias toward young glasses in the upper 10 cm of the lunar regolith, and hence is related to the sampling strategy undertaken by the Apollo astronauts.
LIU ET AL.

10.1029/2020JE006708
10 of 17 Based on the Cratered Terrain Evolution Model (CTEM), Huang et al. (2018) investigated the distribution of impact glass spherules with impact gardening. They also suggested that an increased abundance of younger spherule ages could be explained by sampling bias. Aiming to trace the diffusion of glass spherules, for each impact event, the newly generated impact melt that stayed inside the cavity and that was ejected was regarded as "melt body" and "glass spherules," respectively. The distal ejecta with various patterns were considered. Their results suggested that distal ejecta could be the main source of glass spherules, and the ejecta close to the crater rim may not abundantly mix spherules. In our model, we do not consider the distal ejecta and all the ejected material is taken to be emplaced within five radii from the crater center. In addition, our model does not distinguish the various types of impact melt (i.e., rocks, spherules, or shards). If, taking the same assumption as the model of Huang et al., only the ejected melt product can result in spherules, the spherule abundance should be a fraction of the total generated impact melt. Due to the nearly constant ratio of the ejected melt to the total melt for craters with a diameter of 10s-100s m (∼0.75, Cintala & Grieve, 1998), the relative abundance of the differently aged melt when the total melt products are considered should be comparable with that when only the ejected part is considered. Therefore, for the statistical melt distribution, both our model and the estimate of Huang et al. present a similar trend of the younger age bias. However, since no distal ejecta is considered in our model, the regional spatial distribution may be different from the result of Huang et al.
Although our model does not distinguish whether the impact melt products are glass spheres or shards, our results show clearly different distribution features of impact melt derived from differently sized craters, which is found to be comparable to the distribution of dated impact glass of different types. The distribution of melt associated with small craters (<100-m diameter; Figure 6a) is remarkably similar to that of radioisotopic ages of impact glasses spherules younger than 3 Ga (Figure 6b). If >100-m craters, which produce significantly more amount of melt than smaller ones, produce abundant glass spherules one would see randomly distributed age peaks like in Figure 4b. However, such peaks are absent in the glass spherules record of the past 3.0 Ga (Figure 6c). This suggests that the datable impact glass spherules in the returned lunar samples may have mainly been sourced from craters smaller than 100 m ( Figure 6a) and that >100m craters do not contribute abundant spherules. As indicated by the crater distribution in Figure 2b, the small craters that are the dominant source of glass spherules may be smaller than 20 m. This is in line with chemical compositions of the impact glass spherules (Zellner & Delano, 2015). It was found that impact glass spherules typically are representative of the local regolith in which they were collected, implying small and local impact events as the source Simon et al., 1982;Zellner et al., 2002). Some impact glasses, however, may not have a local origin but still represent relatively small craters. For example, Norman et al. (2019) suggested that relatively small impact events thousands of kilometers away from the Apollo 16 landing site might have contributed exotic glasses to the regolith at this landing site within the last 500 My. However, in our model, the most distant small crater that contributed material to the target region is only ∼5 km away (Crater 4; see Figure 2d), and most of the craters whose ejecta were transported to the central target region after their formations are located no farther than 1 km away. It is possible that the melt products of small craters are entrained in the ejecta of craters larger than tens of kilometers and then redistributed to the area far away from their origin. Such transport by a large impact is not considered in this model. In addition, the hot ejecta blanket of large craters might have an effect on the entrained older melt products (e.g., on the diffusive loss of argon) of small craters, and this would be visible in the age data.
A more distal origin is likely for impact melt shards found in the lunar regolith. In a study of >1,000 Apollo 16 impact glasses, Delano et al. (2007) showed that impact glass shards more often have chemical compositions not typical for the local regolith making it more likely that they were produced in larger more distant impact events. The age distribution of the impact glass shards from lunar samples is shown in gray in Figure 6d. There is no clear increase <0.5 Ga like there is for the impact glass spherules. Excluding some peaks, the shards of different ages display similar abundances through time. As can be seen from Figure 6c, the distribution of melt generated by craters >5 km is very similar to that of the dated glass shards. While Figure 6a indicates that <100-m craters may mainly contribute to the observed age distribution of impact glass spherules rather than impact shards, Figure 6c further suggests that >5-km craters may contribute a significant amount of glass shards in the regolith. This supports the notion from the compositional data that the datable impact glass shards in the return lunar samples may mainly be derived from larger craters.
We provide constraints for the likely main source of datable impact glass collected from the surface and relate it to the crater size. Each impact event can generate different forms of impact glass, but the relative abundance of spherules to shards could be different among differently sized craters. Impact glass spherules form as molten droplets entrained within the excavated material. Impact glass shards form when the existing impact glass spherules are destroyed by the shock waves of subsequent impacts. The capability of destroying spherules depends on the impact energy which is related to the impactor size, thus the crater size. Cratering efficiencies decrease with impactor size and the reach of the shock wave is greater for larger craters, and hence there would be larger shock zones that have magnitude exceeding the strength of impact glass spherules generating more shards. Therefore, the shards distribution is more likely related to large impacts and a significant amount of shards could be derived from a distance. Smaller impacts destroy fewer spherules and hence during the recent period when small impacts are predominant, the majority of their generated spherules could be shielded from destruction and better preserved in the near-surface. It is worth noting that we here do not assume that the K-Ar clock is fully reset during the breaking of spherules into shards. When a sphere is smashed into shards it is likely that the resulting shards are not big enough to be dated or that the K-Ar system cannot remain untouched. Thus, the datable shards are likely a small fraction of the original set of spherules that the shards are formed from. Small craters might just not produce enough original spherules for them to make up a sizable part of the shard population.

Implications for Sample Interpretations
While most impact glass ages record events <3 Ga (Zellner, 2019), our model predicts that the large volumes of melt products that are present in the surface regolith should have formed during the basin-forming period ( Figure 4b). As impact glasses only make up about ∼3-5% of the regolith (McKay et al., 1991), it is possible that a large part of these impact melt products can be found as melt rocks in the regolith. However, due to the small size of fragments in lunar regolith samples it is both difficult to identify impact melt rock fragments as well as to reliably date their formation. In most cases where impact melt rocks have been identified and dated in Apollo and Luna soil samples, the studies themselves concluded that the data did not constrain the nature of the impact event in which they were generated due to the difficulty in determining the provenance (Barra et al., 2006;Kirsten et al., 1973) or that the measured age might record a later thermal event (Fernandes et al., 2013;Podosek et al., 1973;Schaeffer & Husain, 1973). The case that we could find in which a study assessed an age to be reliable and interpreted it to be the formation age of the impact melt, is that of an impact melt fragment (3,895 ± 17 Ma) from Luna 20 by Swindle et al. (1991). A factor influencing the low probability of finding datable samples generated during the basin-forming period could be the small grain size of surviving fragments. These old melt products would mostly be pulverized by later cratering events, effectively grinding them into particle sizes too small to be dated using conventional dating techniques. Mean grain diameters in Apollo regolith samples range typically around 35-55 μm (McKay et al., 1991). Hartmann (2003) proposed high premare cratering rates and calculated extreme physical damage of lunar materials. He suggested that a melt generated 4.1 Ga ago would not be expected to survive until the present because >100 m of regolith has been pulverized to a mean particle size of around 60 μm. Due to the low chance of preserving old impact melt rocks large enough to date, it is not possible to assess whether there are indeed abundant old melt rock fragments in the Apollo and Luna regolith samples as our model predicts. Recently, with the development of dating techniques, old melt with small grain size has been dated (Crow et al., 2016;Vanderliek et al., 2018). According to our estimate, we expect more old melt contained in the returned samples could be found and dated in the near future.
Lunar meteorites are the next best source of samples for assessing the age distribution of impact melt in the upper regolith layer, as most of them come from the upper few meters of the regolith (e.g., Vogt et al., 1991;Warren, 1994). As our modeling results show the abundance of old impact melt increasing with regolith depth (Figure 5), one should expect an even stronger dominance of older impact melt in meteorites due to their deeper origin compared with the returned regolith samples. However, this predicted dominance appears to be absent when looking at the compilation of impact ages of meteorites by Fernandes et al. (2013) or Michael et al. (2018) (Figure 4a). As with the soil samples, dates that have been interpreted to represent impact events can represent resetting but not necessarily the melting process. We, therefore, further filtered those two compilations, keeping only the dates which the original studies interpreted to reflect the melting event (the gray histograms in Figure 4a; Cohen, 2008;Cohen et al., 2005;Daubar et al., 2014;Fernandes et al., 2000Fernandes et al., , 2009Gnos et al., 2004;Joy et al., 2011;Nyquist et al., 2006). The original and filtered data represent an upper and the lower limits of melt products that have been successfully dated, respectively. The filtered data show a higher relative abundance of melt components with ages of 2.5-3.8 Ga reflecting the increasing frequency of impacts around this period, which is consistent with the predicted abundance ( Figure 4b). However, even in the filtered data, no predominance of melt older than 3.8 Ga appears. Thus, we predict that lunar meteorites should contain a significant quantity of as yet unidentified old impact melt products. Similarly to the regolith samples, the reason for this current underrepresentation in the published data could be a lack of measurable argon (or other elements) and/or the loss of the original age information because of the small size of the surviving fragments. In addition, Figure 4b shows that the melt generated by large craters during the Moon's middle age (∼2-3 Ga; gray bars in Figure 4b) volumetrically eclipses the melt from the small craters (black bars in Figure 4b). But very few impact melts were suggested to have been derived from large impacts with these ages. However, for melt from small craters (the black bars), the decreasing abundance with age does reflect well the distribution of small craters (Kirchoff et al., 2013) and the distribution of ages of various kinds of samples (Zellner, 2017). It may result from the insufficient sampling of the melt from large impacts. The high-frequency small impacts result in most of the lunar surface containing their melt products; in contrast, during the late period, only a small part of the lunar near-surface has been gardened by larger impacts (0.7% Ga −1 , see Section 2). Therefore, although the greater amount of melt has been generated by large impacts, its smaller distribution area reduces its probability of being contained in the returned samples that were derived from the regional surface.

Implications for Sampling Strategy
Despite the lower impact flux and decrease in impactor size, the recent and smaller events dominate the gardening of the top lunar surface. We predict that there should still be abundant unheated ancient melt of the megaregolith in the top surface, although it has been intensively gardened by small impacts (Figure 4).
The returned lunar regolith samples were mainly collected from the near-surface. The origin of the lunar meteorites is typically deeper, but still they are probably mainly derived from the surface regolith (McKay et al., 1991). Samples that have not been ground to very small grain size can preserve age information (Section 4.2). Combining the predicted melt sitting depth ( Figure 5) and the age distribution of melt products that have been dated (Figures 4a, 6b, and 6d), a rough picture of the distribution of datable melt products at varying depth emerges (Figure 7): • <3 Ga: Our modeling results show that young melt is abundant in the near-surface, with some occurring at a greater depth of a few meters. In the returned samples, many impact melt products younger than 3 Ga have been dated (mainly impact melt glasses Figure 6) and some have also been found in meteorites ( Figure 4). The comparison of our results with dated impact melt products indicates that datable young impact melt glasses are representative of the expected age distribution of young melt products.
LIU ET AL.

10.1029/2020JE006708
13 of 17 Figure 7. A schematic showing the abundance of melt that can or cannot be radioisotopically dated with respect to depth. Note that the boundary between deeper regolith and megaregolith is not well-defined. The size of markers indicates the relative size of the datable melt products (not to scale).
• 3-3.8 Ga: Our simulation results present that melt products of 3-3.8 Ga are relatively more abundant than the younger melt in the near-surface although their great abundance is at greater depth ( Figure 5). In the returned regolith samples (i.e., at shallow sampling depth), some melt products generated between 3 and 3.8 Ga have been detected. This indicates that in the near-surface regolith a handful of the melt products that should be present from this time remain datable. In addition, these melt products are often found in meteorites (i.e., from somewhat the deeper regolith layer; Figure 4). The calculated abundance at varying depth ( Figure 5) suggests that >25% of melt products are deposited in the top ∼10 m; and >50% are at a depth of 10-100 m. This suggests that more of these melts might be found in the deeper regolith and megaregolith where they also would have a higher chance to remain datable, as those materials at those depths have been less affected by reworking (Figure 3). • >3.8 Ga: Our simulations suggest that the volume of the melt products remaining from the basin-forming period in the top surface should be much larger than that of the young melt. In addition, the calculated abundance at varying depth ( Figure 5) shows that about 25% of those old melt products are distributed in the top 100 m. However, dated samples from both the regolith samples and the meteorites are few for various reasons (see Section 4.2 above). That is, the likelihood of extracting their ages from near-surface samples appears (at least with conventional dating methods) low. Therefore, we predict that despite being present in great volume, the datable melt products of this age are limited at both the surface sampling depth (<10 cm) and in the plausible greater depth of meteorites' origin.
For future missions, to avoid the influence of small craters the sampling depth should be at least greater than the top meter ( Figure 3a). If the datable melt older than ∼3.9 Ga is the target, the samples should be collected from the depth deeper than the origin of most lunar meteorites (i.e., the deeper regolith). The typical regolith thickness of the highlands is about 10 m. The maria are not only covered with the soft regolith, but also have been filled with volcanic mare deposits of various thicknesses (McKay et al., 1991). However, using current sampling tools the megaregolith beneath the regolith layer is difficult to reach. The deepest core drilled during the Apollo 17 mission was ∼3 m. Due to such shallow sampling depth, the majority of the collected datable melt products is of small-crater origin and the older parts are likely to have been significantly gardened by small impact events. A smarter approach when aiming to collect old melt is to take advantage of the reexcavation zones of craters close to the landing site. If the surface is not covered by the mare material, those craters should be at least >100 m in diameter to assure the reexcavation of a significant amount of the deeper-seated old melt. More importantly, these targeted craters should have formed recently, to increase the chance of the collected material preserving a measurable formation age. During the Apollo missions, this approach has been shown to be successful. For example, rock material from the ejecta blanket of the young (50 Ma) North Ray crater contained impact melt products that record pre-3.9 Ga impacts preserved in stratigraphically deeper layers of the lunar regolith (Fernandes et al., 2013;Norman & Nemchin, 2014;Norman et al., 2016). If the largest recently formed crater near the landing sites is tens of meters in diameter, its excavation depth is about several meters, likely smaller than the local regolith thickness, and hence the majority of its reexcavated melt products older than 3.9 Ga are probably difficult to date. If the landing site is located on the maria, the crater should be even larger to excavate material beneath the mare deposit. Since the typical mare thickness is hundreds of meters (McKay et al., 1991), the crater should be larger than tens of kilometers, which, however, increases the chance that the radioisotopic systematic of the excavated material is affected during excavation. In addition, as the majority of impact melt has been sequestered as a lens beneath crater floors, reexcavated material from within a crater or basin will increase the chance to collect associated melt products.

Conclusions
For this study, a numerical model to investigate the mixing by small impacts is developed and added to the model of Liu et al. (2020), with the aim to produce a more comprehensive picture of the melt distribution in the lunar regolith. The simulation results provide insight into our knowledge of the recent lunar impact flux, the plausible main source of melt products (glass spherules and shards) in the near-surface, and the sampling strategies of the future missions that aim to collect the old datable impact melt.
The majority of the melt generated by craters <5 km is distributed in the top 1 m of the regolith. These craters mainly mix the material at and close to the surface, which creates abundant melt with age <0.5 Ga in the very top surface (∼10 cm); at greater depth the melt ages are more representative of the impactor flux.
The results indicate that the higher frequency of impact glass spheres <0.5 Ga in the returned lunar surface samples was likely caused by the extensive reworking by small impacts, and not an increase in the impact rate during that time. In addition, by comparing the distribution of melt generated by craters <5 km and by the larger craters with that of impact melt glass spherules and shards from the returned samples, the possible main origin of the datable impact melt glass in the near-surface can be constrained. To reconcile the observed age distribution of lunar impact glass samples with our model results requires that the melt products of <100-m craters are the main source of datable glass spherules and that >100-m craters do not contribute abundant datable spherules. The modeling results also show that at locations with a coverage of thick ejecta from large craters the earlier-formed melt of small craters should be buried or diluted, but the distribution of later-formed melt should be identical to the distribution at locations without such an ejecta blanket. The abundance of the affected melt that could be subsequently gardened to the top surface is related to the thickness and the age of the ejecta and the relative position to the source crater.
By combining the mixing effects of both large and smaller craters, the abundance of ancient melt in the near-surface is estimated. The results show that even in the top meter where the majority of the melt generated by small craters is deposited, the ancient melt from the megaregolith should still be predominant. However, this is not concordant with the current age record of regolith samples. We think that the grinding of near-surface materials by the intense regolith gardening makes it difficult to find particles of a suitable size for conventional radioisotopic dating, an effect not taken into account by the model. But with the development of the dating techniques, we expect more old melt in samples could be dated. By combining the predicted melt abundance and the distribution of melt products that have been dated, a picture of the distribution of datable melt at varying depths is inferred. For future sampling mission, whose goal is to collect datable old material, ejecta blankets of young and large craters (>100 m on highlands; >10 km on maria) that have excavated deep megaregolith appear to be the most fruitful targets.