Re-examining the main asteroid belt as the primary source of ancient lunar craters

It has been hypothesized that the impactors that created the majority of the observable craters on the ancient lunar highlands were derived from the main asteroid belt in such a way that preserved their size-frequency distribution. A more limited version of this hypothesis, dubbed the E-belt hypothesis, postulates that a destabilized contiguous inner extension of the main asteroid belt produced a bombardment limited to those craters younger than Nectaris basin. We investigate these hypotheses with a Monte Carlo code called the Cratered Terrain Evolution Model (CTEM). We find that matching the observed number of lunar highlands craters with Dc~100 km requires that the total number of impacting asteroids with Di>10 km be no fewer than 4x10-6 km-2. However, this required mass of impactors has<1% chance of producing only a single basin larger than the ~1200 km Imbrium basin; instead, these simulations are likely to produce more large basins than are observed on the Moon. This difficulty in reproducing the lunar highlands cratering record with a main asteroid belt SFD arises because the main belt is relatively abundant in the objects that produce these"megabasins"that are larger than Imbrium. We also find that the main asteroid belt SFD has<16% chance of producing Nectarian densities of Dc>64 km craters while not producing a crater larger than Imbrium, as required by the E-belt hypothesis. These results suggest that the lunar highlands were unlikely to have been bombarded by a population whose size-frequency distribution resembles that of the currently observed main asteroid belt. We suggest that the population of impactors that cratered the lunar highlands had a somewhat similar size-frequency distribution as the modern main asteroid belt, but had a smaller ratio of objects capable of producing megabasins compared to objects capable of producing ~100 km craters.

impact basins as a constraint on the ancient impactor population of the Moon. We find that matching the observed number of lunar highlands craters with D c ≃ 100 km requires that the total number of impacting asteroids with D i > 10 km be no fewer than 4 × 10 −6 km −2 . However, this required mass of impactors has < 1% chance of producing only a single basin larger than the ∼ 1200 km Imbrium basin; instead, these simulations are likely to produce more large basins than are observed on the Moon. This difficulty in reproducing the lunar highlands cratering record with a main asteroid belt SFD arises because the main belt is relatively abundant in the objects that produce these "megabasins" that are larger than Imbrium. We also find that the main asteroid belt SFD has < 16% chance of producing Nectarian densities of D c > 64 km craters while not producing a crater larger than Imbrium, as required by the E-belt hypothesis. These results suggest that the lunar highlands were unlikely to have been bombarded by a population whose size-frequency distribution resembles that of the currently observed main asteroid belt. We suggest that the population of impactors that cratered

Introduction
The heavily-cratered highlands of the Moon are uniquely suited as a terrain for studying the early bombardment history of the inner solar sys-tem. Little post-bombardment geological processing has modified highlands craters over the ∼ 3.5 Gy since the majority of them formed. This is in contrast to the heavily cratered terrains of Mars and Mercury which are much more heavily modified by geologic processes other than cratering. We also have at present a substantial sample collection from the Moon, including highlands rocks, thanks to both the lunar meteorite collection and the Apollo and Lunar programs (Vaniman et al., 1991). The combination of impactdominated geology and an abundant and diverse sample record means that the lunar cratering record has been the primary source for absolute-age crater chronology systems. It is also the best record currently available for testing models for the bombardment history of the inner solar system.
It has been hypothesized that the ancient lunar highlands craters were formed by impacts from a population that was identical to that of the main asteroid belt, based on similarities between the size-frequency distribution of asteroids and the size-frequency distribution of lunar highlands impactors derived using a crater size scaling relationship (Strom et al., 2005). This correlation between the lunar highlands impactors and the main asteroid belt was interpreted to suggest that some mass-independent process perturbed main belt asteroids onto terrestrial planet-crossing orbits. Because gravitational accelerations are mass-independent (Galilei, 1638), it was hypothesized that the sweeping of mean motion and secular resonances by the migration of giant planets (Fernandez and Ip, 1984;Malhotra, 1993;Hahn and Malhotra, 1999;Tsiganis et al., 2005) was responsible for generating the impactor populations that cratered the lunar highlands. Dating of lunar samples and evidence from geological superposition suggests at least two large lunar basins (where "basin" here refers to a crater with final rim-to-rim diameter D c > 300 km), Imbrium and Orientale, were formed 600-800 My after solar system formation (Wilhelms et al., 1987), and therefore the migration of the giant planets needed to be delayed by that amount of time in order to produce at least these basins from the asteroid belt.
A scenario for giant planet orbital evolution dubbed the "Nice model" was developed, in which the giant planets initially formed in a tightly packed, quasi-stable configuration with a massive icy planetesimal disk beyond Neptune that was disrupted when Jupiter and Saturn crossed a mean motion resonance Morbidelli et al., 2005;Gomes et al., 2005).
The Nice model provides a mechanism to allow a delay between planet formation and the onset of the Late Heavy Bombardment (LHB). It was estimated that the Nice model could have produced a bombardment of ∼ 10 22 g of asteroids onto the Moon, plus an equivalent amount of comets . This mass is enough to produce all of the of observed lunar basins on the Moon (c.f. Ryder, 2002). Assuming an asteroid mean density of ρ i = 2.5 g cm −3 and using a size-frequency distribution of asteroids obtained by the Wide-field Infrared Survey Explorer spacecraft (WISE) (Masiero et al., 2011), this corresponds to a number density of D i > 10 km impactors of approximately 9.7 × 10 −7 km −2 .
However, the original Nice model calculations for the magnitude of the lunar bombardment produced by disrupting the main asteroid belt was challenged when it was shown that if the migration rates of Jupiter and Saturn were too slow they would have caused the terrestrial planets to become disrupted (Brasser et al., 2009;Agnor and Lin, 2012). Also, a migration rate that was too slow would have altered the distribution of asteroid eccentricities (Minton and Malhotra, 2011) and inclinations (Morbidelli et al., 2010) in ways that are not supported by observation. A more rapid evolution of Jupiter and Saturn compatible with these inner solar system constraints occurs in a small subset of Nice model simulations dubbed the "Jumping Jupiter" scenario when the ice giants (Uranus and Neptune) have close encounters with Jupiter and Saturn (Brasser et al., 2009). However, a more rapid migration of the giant planets would not have produced the required amount of mass needed to crater the lunar highlands and produce all of the observed lunar basins. It was hypothesized that the main asteroid belt may have had a primordial inner extension, dubbed the E-belt, which is a region that would only have become unstable after giant planet migration occurred and the ν 6 secular resonance arrived at ∼ 2 AU . The destabilized E-belt, along with a small contribution from the main asteroid belt, would produce ∼ 1/3 of the lunar basins as a result of rapid giant planet migration under the Jumping Jupiter scenario, or the ∼ 14 youngest basins beginning with Nectaris (Fassett et al., 2012). The more limited view of the LHB suggested by the E-belt model postulates that the pre-Nectarian terrains may have been created by leftover planetesimals from the epoch of terrestrial planet formation , and possibly also include a contribution from cometary impactors, which should dominate total impactor mass , although it is not clear why a non-asteroidal signature is not apparent in the lunar cratering record .
Putative impactor remnants identified in lunar samples also favor asteroids over comets in both highly siderophile elements (Kring and Cohen, 2002) and oxygen isotopes (Joy et al., 2012). A more extensive review of this topic may be found in Fassett and Minton (2013).
Here we revisit the hypothesis that the lunar highlands were cratered by main belt asteroids. A powerful Monte Carlo code for studying the evolution of cratered terrains has recently been developed (Richardson, 2009).
The code is named, simply enough, the Cratered Terrain Evolution Model (CTEM). We used CTEM to model the cratering history of the lunar highlands by bombarding simulated lunar surfaces with a main belt asteroid impactor population. Our goal is to quantify the amount of cratering required to reproduce the observed lunar highlands crater record. We do this by performing a large number of CTEM simulations of impactors bombarding a lunar-like surface and then comparing the outcomes with published crater counts for the lunar highlands. For this study we use a catalog of craters produced by the Lunar Orbiter Laser Altimeter (LOLA) aboard the Lunar Reconnaissance Orbiter spacecraft . Our impactor population uses the most up-to-date size-frequency distribution (SFD) for the main asteroid belt as obtained by the Wide-field Infrared Survey Explorer (WISE) (Masiero et al., 2011). For impactor sizes below 5 km, where WISE data becomes unreliable due to biases, we use the main belt SFD from the Sub-Kilometer Asteroid Diameter Survey (SKADS) (Gladman et al., 2009), assuming 9% albedo. This paper is divided into two main sections. In the first section we will describe in detail our methods, including a description of the CTEM code as well as our efforts to calibrate it with a human crater counter. In the second section, we will describe the results of applying the CTEM model to the problem of the impact history of the lunar highlands.

Methods
For this project we employ a Monte Carlo cratering code called the Cratered Terrain Evolution Model (CTEM) (Richardson, 2009). Numerous improvements have been made to CTEM since the Richardson (2009) study, which we will describe in detail here. First, several performance improvements have been made to the code, including the addition of parallelization, allowing us to explore a much wider parameter space for a given study than was available previously. The code was also re-written in Fortran 2003, and allows a much greater degree of extensibility and modularity than the earlier Fortran 77 version. Lastly, we have significantly revised the crater-counting algorithm within the code and performed a comprehensive calibration study with a human crater counter (co-author Fassett), allowing us to better apply the computer model results obtained from CTEM to the interpretation of remotely observed planetary surfaces.
CTEM models a given planetary surface as a grid of points (heretofore referred to as pixels), where the number of pixels in the grid is set by the user. Each pixel contains a variety of data, one of which is the current elevation at that location, allowing the construction of a Digital Elevation Model (DEM) of the surface. For each impact that is modeled, the DEM of the surface is altered to represent the topographic changes brought about in the formation of a crater. Because CTEM models the topography of craters, rather than idealizing them as simple shapes such as rings or circles, a number of important crater degradation processes are naturally produced by the code, including unstable slope collapse, impact ejecta, and diffusive erosion.
Determining whether a crater is countable or not is done by measuring how much the topography of a given crater deviates from its original shape.
In order to quantify crater "countability" within the code, we performed a crater counting calibration study with CTEM generated DEMs by supplying them to a human crater counter (co-author Fassett), who analyzed CTEM-generated DEMs in a way similar to how craters were counted using LOLA-generated lunar DEMs. Therefore the crater counting algorithm of the CTEM code has been calibrated using the same process that our comparison data set was generated from.
Here we summarize all of the physical processes that CTEM models, and further details regarding the code's capabilities are described in Richardson (2009). Once an impactor size, velocity, and impact angle are determined, the following sequence of steps is performed to model the resulting crater and its effects on the topography of the simulated terrain: 1. Perform a crater scaling calculation to determine the dimensions of the resulting crater.
2. Place the crater at a random location on a surface and modify the terrain to reflect the final crater shape. We use a repeating boundary condition to mitigate against edge effects.
3. Emplace ejecta onto the surface and degrade the pre-existing terrain affected by ejecta deposition.
4. Collapse any slopes that are above the angle of repose.
Once at the end of a run, as well as periodically throughout the run, CTEM will examine the simulated surface and tally up the number of observable craters. In Subsection 2.1 we will we will describe the crater scaling relationship used by CTEM as well as its method for emplacing ejecta onto the surface. Because including estimates of the number of large lunar impact basins is a critical component of our present work, in Subsection 2.2 we describe how we incorporated recently published scaling relationships for lunar basins based on computer hydrocode modeling into CTEM. Next, in Subsection 2.3 we will describe a number of ways that craters may degraded over time on a CTEM-generated surface. Then, in Subsection 2.4 we will describe the study that we performed in order to calibrate the crater recognition algorithm within CTEM with a human crater counter. Finally, in Subsection 2.5 we will describe the impactor population model that we adopted for this study.

Crater and ejecta scaling model
The impact crater volume scaling relationships that we use to relate the size of a primary impactor to the size of its resulting crater on a particular target surface, given a set of projectile and target material properties and impact parameters, are reviewed in Holsapple (1993). Previously, most applications of such relationships dealt strictly in either the gravityor strength-dominated cratering regimes; for example, Vickery (1986). However, cratering on small scales or on small target bodies often falls into neither regime: gravity and target strength are both important to the size of the final crater (Ivanov, 2001;Richardson et al., 2005). We have therefore adopted the general solution to the transient crater volume scaling relationship given by Eq. 19 in Holsapple (1993), which includes both gravity and strength terms.
The application of a general solution to the crater volume scaling relationship permits us to also utilize the general solution to the ejecta velocity scaling relationships developed in Richardson et al. (2007). These revised ejecta scaling relationships permit us to directly compute surface ejection velocity as a function of projectile and target material properties, impact parameters, and distance from the impact site, all the way out to the transient crater rim. This analytical solution is applied to a discretized model (by a factor of 10 3 to 10 4 ) of the excavation flow-field's hydrodynamic streamlines (Maxwell and Seifert, 1974;Maxwell, 1977), in order to compute excavated mass as a function of distance from the impact site and ejection velocity. Finally, these discrete ejected mass segments are followed in post-ejection flight through a set of ballistics equations in order to compute ejecta blanket thickness as a function of distance from the final crater rim, as described in Richardson (2009). This scaling-relationship based excavation-flow properties model is quite general in nature, and can be applied to craters ranging in size from the laboratory (Richardson, 2011) to impacts on small solar system bodies (Richardson et al., 2007) to large, planetary scales (Richardson, 2009). This allows us to model ejecta deposition over a wide range of gravity and projectile/target material properties and impact parameters.
The principal equations utilized in this impact ejecta model are as follows (from the Appendix to Richardson (2011)). The central analytical expression of the "excavation flow properties model" describes the effective particle ejection velocity v ef as a function of radial distance r from the impact site as follows (from Secs. 2.2 & 2.3 of Richardson et al. (2007) where ρ t is the target density, g is the gravitational acceleration, andȲ is the effective target strength with respect to crater excavation. The uncorrected ejection speed v e is given by the gravity-scaled ejecta speed equation derived by Housen et al. (1983): where R g is the gravity-scaled transient crater radius; that is, the radius that would result if the target was strengthless. µ is a commonly used exponential material constant (Holsapple, 1993), and the constant C vpg is given in Richardson et al. (2007) as: This expression contains constant C T g , which has an experimentally determined value of 0.8-0.9 (Schmidt and Housen, 1987;Holsapple and Housen, 2007) as discussed in Sec. 2.2 of Richardson et al. (2007).
The constant C vps in the strength term of Eq. 1 is given by: where R s is the transient crater radius due to the effects of both gravity and strength. The transition strength Y t (between gravity and strength dominated cratering) is expressed as: where v i is the vertical component of the impactors speed, a i is the impactor's mean radius, and ρ i is the impactor density.
In order to estimate the size of the transient crater produced by a particular impact, we make use of the general solution to the crater-size scaling relationship developed in Holsapple (1993), which contains the effects of both gravity and strength: where K 1 , µ, andȲ are experimentally derived properties of the target material. We adopt K 1 = 0.22 and µ = 0.55 for our study, as these are expected values for hard rock target material based on experiments (Holsapple, 1993).
In practice, the effective target strengthȲ is roughly a measure of "cohesion," and usually lies somewhere between the laboratory-measured tensile and shear strengths of the material (Holsapple, 1993;Holsapple and Housen, 2007), for low-to medium-porosity targets. In high-porosity targets, the effective target strengthȲ can be thought of in broader terms: as the energy per unit volume (rather than the force per unit area) required to both crush pore spaces and break the material apart for excavation. Because in our study we are only comparing our model against D c > 20 km craters, our scaling relationship is insensitive to effective target strength,Ȳ .
The transient crater volume V can be related to the more easily measured transient crater (rim-to-rim) diameter D t or radius R t by: where we assume that the transient crater depth H t is roughly 1/3 its diameter D t : in experiments this is somewhat variable, between 1/4 and 1/3 (Schmidt and Housen, 1987;Melosh, 1989 as described in Chapman and McKinnon (1986) and Melosh (1989). For complex craters the conversion follows a power-law, as described in Croft (1981), Chapman and McKinnon (1986), Melosh (1989), and Schenk et al. (2004). This power-law begins at a simple crater diameter known as the simple-to-complex transition point D tr , and where this transition point varies inversely with the surface gravity of the body involved (D tr ∝ 1/g).
Fitting the data shown in Fig. 12.18 of Schenk et al. (2004), we developed and expression for the simple-to-complex transition for silicate rock: where the gravity g is given in terms of m s −2 , and the transition crater diameter D tr is given in term of meters.
In addition to the above method for determining the diameter at which simple craters transition to complex craters, we also need a method for computing the size of the final complex crater given the transient crater diameter from Eq. 7. For a silicate rock body, the final diameter of a complex crater is given by: where the simple crater diameter D s is computing using Eq. 8, and the above expression is used only when D s > D tr ; for simple craters D c = D s . The silicate body power-law exponent is taken from Croft (1985).
Finally, we set a maximum depth for the floors of complex craters to a depth, h max given by Pike (1977): This creates complex craters with a characteristic flat-floored shape. The current version of CTEM does not attempt to model central peaks or peak rings.

Scaling relationships for large lunar impact basins
Large lunar basins present a challenge for crater scaling relationships, as the formation of a lunar basin, which we define as any crater with a final diameter D c > 300 km, involves energies far beyond those approachable in laboratory or field studies. Even the largest known impact craters on Earth are somewhat below the size of a "small" lunar basin. Nevertheless, a number of advancements have been made in recent years in understanding the processes involved in lunar basin formation using numerical hydrocodes. Potter et al. (2012b) developed a scaling relationship for lunar basins based on hydrocode simulations that relates the size of the transient crater to a crustal feature called the "crustal annulus," which is observed in the crust beneath lunar basins (Wieczorek and Phillips, 1999;Hikida and Wieczorek, 2007;Melosh et al., 2013). The size of a basin's crustal annulus is often easier to determine than its rim-to-rim size. Potter et al. (2012b) also showed showed that the transient crater diameter was sensitive to the thermal profile of the lunar upper mantle, and developed two scaling relationships to account for changes to the lunar thermal gradient over time as well as the major hemispherical dichotomy of the Moon. They define two thermal profiles, TP1 and TP2, where TP1 represents a weaker, warmer upper mantle, and TP2 represents a colder, stronger upper mantle. Because the near-side crust is thinner and possibly warmer due to higher abundances of radioactive nuclides within the Procellarum KREEP terrain (Jolliff et al., 2000), this thermal dependence means that the size of basins that form on the near side tend to be larger than those on the far side for a given impactor size.
This effect has been seen in data returned by the Gravity Recovery and Inte- where diameters are given in units of km. The cold/strong TP2 scaling relationship is given as: Although CTEM has the capability of modeling crustal thickness, our crater counting algorithm is calibrated for surface topography, and therefore we need an estimate of the final crater rim. Potter et al. (2012b) do not explicitly relate the size of the crustal annulus to the size of the final rim, however they do provide measurements for three well-studied craters, Orientale, Imbrium, and Serenitatis. Imbrium was modeled using TP1, and Orientale was modeled using the TP2. The Orientale final rim diameter of 930 km used in our crater catalog is ∼ 56% larger than the crustal annulus.
Imbrium and Serenitatis rims are ∼ 30% larger and ∼ 10% larger than their respective crustal annuli. Therefore to obtain the final crater rim diameter, D c , from the crustal annulus diameter, D ca : Because of the statistical nature of our study, we choose to apply either the TP1 or TP2 scaling relationship for a given crater with a 50% probability to reflect the hemispherical dichotomy of the Moon. We apply our basin scaling relationship for impactors with D i > 35 km. To test our basin scaling relationship, we generated a test crater the size of Orientale. Because h e = 2.9 km (4/R cr ) −2.8 In order to reproduce a crater the size of Orientale with D c = 930 km and an ejecta blanket thickness at the rim of 2.9 km, using cold/strong mantle thermal profile TP2, we require a D i = 75 km impactor with a velocity of 15 km s −1 , and scaling parameters of K 1 = 0.22 and µ = 0.55. This is consistent with the velocity and size of the Orientale impactor of 50-80 km found by Potter et al. (2013), using the iSALE hydrocode and constraints from crustal thicknesses derived from GRAIL spacecraft data. Using a scaling relationship from Holsapple (1993), Ryder (2002) derived a 4 × 10 20 g impactor traveling at 20 km s −1 for Orientale, which corresponds to a 64 km impactor for ρ i = 2.5 g cm −3 . Our adopted scaling relationship for basins therefore errs on the upper end of previously estimated ranges of Orientale-impactors, which makes it a conservative model for our problem. We will discuss the ejecta blanket profile of this test crater in more detail in Section 2.3.3.
As large basin impacts on the scale of ∼ 2500 km South Pole-Aitken basin are important for our study, we tested our scaling relationship for these large megabasins. Potter et al. (2012a) found using hydrodcode modeling that the best match to SPA was found for D i = 170 km traveling at 10 km s −1 and impacting into mantle with a high thermal gradient (more similar to TP1 than TP2). For this size impactor and velocity, using the TP1 crater scaling relationship we obtain a final diameter of of D c = 2033 km, and for TP2 we obtain D c = 1355 km. This impactor size is also consistent with other hydrocode simulations of the formation of SPA (Collins and Melosh, 2004).
Therefore, CTEM reproduces basins at the scale of SPA, and in fact we may be somewhat conservative in the size of basins produced by impactors similar to the one that Potter et al. (2012a) modeled.

Crater degradation modeling with CTEM
In order for CTEM to properly model a surface in or near the state of cratering equilibrium, we must first calibrate the crater detection algorithm to ensure that the craters that are counted by CTEM are equivalent to the craters that a human would count on a real surface with an equivalent cratering history. Here we first describe how CTEM handles different mechanisms of crater erasure, and then describe the work we have done to calibrate CTEM's crater identification algorithm with a human crater counter.
The difficulty in implementing accurate crater obliteration modeling may be restated as the difficulty in identifying "countable" model craters in a way that mimic what a human crater counter would also consider countable on real terrain with an equivalent bombardment history. This difficulty is one that has been encountered by all researchers that have attempted to con-struct Monte Carlo cratering codes. A limitation in past attempts has been the confluence of limited computing power and a problem that requires a great deal of computer resources to study, due to the huge range in spatial scales and the power law nature of impacting populations. We take advantage of advances in computing power in order to model the cratering process as a sequence of topographic changes produced on a simulated surface represented by a finite grid of elevation points, or digital elevation model (DEM).
This allows us to model the processes of crater formation and obliteration naturalistically.
For the global and regional scale lunar cratering simulations performed for this study, old craters are obliterated by four primary mechanisms, each of which is modeled within the present version of CTEM: 1) Cookie-cutting, where new craters destroy craters directly by forming on the same terrain, 2) Sandblasting, where numerous small craters erode an old crater enough that the old crater is no longer recognized, 3) Seismic shaking, where the surface vibrations produced by each new impact destabilize and degrade slopes and craters, and 4) Ejecta burial/ballistic sedimentation, where ejecta deposits fill in and erode old craters. In the following subsections we will describe how CTEM handles the three most important crater obliteration mechanisms for our study (cookie-cutting, sandblasting, and ejecta burial/ballistic sedimentation). We will also describe how we have tested scenarios wherein one of the three mechanisms dominates over the others. On large bodies, such as the Earth's moon, impact-induced seismic shaking likely plays only a minor role in crater obliteration, particularly on the regional and global scales considered in this particular study (Houston et al., 1973). Therefore we do not employ the seismic shaking model here.

Crater erasure mechanisms: Cookie Cutting
Cookie cutting describes the most direct way that new craters destroy old craters, which is by simple geometric overlap (Woronow, 1977a(Woronow, ,b, 1978.
When an impactor strikes a cratered surface, any pre-existing craters within the final rim of the new crater are obliterated during the crater formation process. CTEM implements cookie cutting in a simple way, but in order to describe it we must first digress and explain a feature of the code that facilitates crater identification during the counting stage.
The surface that CTEM models is represented by a grid of pixels with a repeating boundary. Every pixel on a CTEM surface grid contains a variety of information in addition to elevation. When CTEM generates a crater centered at a particular position on the grid, the pixels interior to a circle of diameter 1.05D c , where D c is the crater's final rim diameter, are tagged with the crater's unique identifier. The value of D c is used as the unique identifier for each crater, as it is unlikely that any two craters will have precisely the same diameter (to double precision). In order for CTEM to successfully identify topographic depressions as possible craters, the code must allow for any given pixel to contain the tags of multiple craters, otherwise a large crater could become completely erased by multiple small impacts despite the original crater being identifiable as as topographic depression. We implement this capability in a layering system, where a single pixel contains multiple crater tags. The existence of a crater's tag on the grid is a necessary, but not sufficient, condition for it to be counted.
If pre-existing craters within the new crater's final rim are deemed to be removed due to cookie cutting, then the tag that identifies the old crater is removed from the pixel. The cookie cutting mechanism generally only applies if the new crater is comparable in size or larger than the pre-existing craters. So when a pixel is tagged with a new crater, the tags for all craters smaller than the new crater within the region where the two craters overlap are removed. We show an example of craters affected by cookie cutting in Small craters that form inside of larger ones don't necessarily erase the larger ones. A large crater can become completely covered in small craters, and still be recognized as a topographic depression and counted, since the smaller craters did not remove the larger crater's tags. Pre-existing craters that are larger than newer overlapping craters can only be removed by substantial topographic alteration.

Crater erasure mechanisms: Sandblasting
The term "sand-blasting" has been used to describe the process by which many small impacts can erode the surfaces of asteroids (Chapman, 1976), but was initially adopted for the analogous process by which numerous small craters erode larger craters to the point where the larger craters are not recognizable as a crater by a human crater counter (Soderblom, 1970;Hartmann, 1984). Human crater counters regularly identify topographic depressions on planetary surfaces as degraded craters. Therefore, in order for CTEM to count as craters those features which would be counted by a human, it must be able to successfully identify sandblasted craters.
When multiple small craters impact the surface occupied by a larger old crater, the new craters transport ejected material and, over time, the topography of the large old crater relaxes back toward the mean elevation of the terrain. This topographic diffusion is seen to occur on the lunar mare (Fassett, 2013). This process can be seen clearly in a CTEM simulation as shown in the sequence of images in Fig

Crater erasure mechanisms: Ejecta burial
The ejecta from a fresh crater can obliterate old craters beyond its rim by burying them in ejecta. Intuitively, it seems obvious that if the ejecta thickness is greater than the depth of a crater, that the crater will be buried, and thus will be invisible to a crater counter. This concept has been employed to estimate the thickness of ejecta deposits from lunar craters ( The reason for the failure of the ejecta to erase craters can be seen in Fig. 4a, which plots the pre-impact and post-impact surface profiles along the centerline of the grid from the simulation shown in Fig. 3a. CTEM has effectively "painted" the surface with ejecta, preserving the topography of the craters beneath. This is clearly not a physically-plausible ejecta deposition model. Real ejecta deposition is highly energetic, and the deposition process scours and transport material some distance away from the original ejecta impact point.
To model the "blanketing" effect of the ejecta blanket, we use topographic diffusion to smooth the terrain prior to emplacing ejecta. We use the diffusion equation of the form where z is the elevation, x and y are the spatial dimensions on the surface, and K d is a diffusion coefficient, which can vary spatially. The diffusion model is discretized on our grid using a 2nd-order central-difference scheme. grid (not shown here). Our empirically determined diffusion coefficient is where b e is the local ejecta thickness and l px is the length of a single pixel. We note that each crater's ejecta profile is computed within CTEM using our excavation and ballistic transport model, and so they do not scale in a simple way with the crater size. This gives us confidence that our model of crater erasure through proximal ejecta burial is robust.

Crater counting calibration study
Periodically CTEM performs a tally step to count all visible craters on the surface. CTEM's crater counting algorithm employs the use of identifying tags for craters that have formed on the surface, as described Section 2.3.1.
This tag system has an advantage over algorithms that attempt to count craters on real surfaces: CTEM knows where all the craters formed. This simplifies crater counting, in that CTEM does not need to employ any kind of feature-recognition algorithm in order to find potential craters. However, CTEM still needs to determine whether or not those craters that have escaped cookie-cutting are still features that a human would identify as a crater.
To accomplish this, we performed a calibration study with a human crater counter in order to determine what measurable qualities of CTEM-generated craterforms predicted their countability.
We based our calibration study on the classic experimental study described by Gault (1970). In that study, craters were produced in 2.5 m wide "sandbox" at the NASA Ames Research Center using projectiles and explosives to simulate heavily cratered lunar terrains. Gault formed the cratered surfaces using six projectiles that produced craters between 5 mm and 17 cm, such that a crater was approximately twice as large as the next smallest. The number of craters of a given size was constrained to be k times as large as the next largest size class. The terrains produced by the experiment were given names depending on the value of k: "Terra Alta" for k = 6, "Mare Exemplum" for k = 10, and "Mare Nostrum" for k = 16. For the Terra Alta experiment, a seventh crater size of 35 cm (produced using explosives rather than a projectile) was added. These crater size distributions may be approximately described by a cumulative power law given by N >D ∝ D −p , where p = log (k + 1) / log 2.
For our calibration study, we produced 5 simulated surfaces designed to reproduced variations of the sandbox experiments of Gault (1970). Our simulations were constructed using either the "Terra Alta" (k = 6) or "Mare Nostrum" surfaces (k = 16), defined by Gault (1970). These simulated surfaces were designed to match as close as possible the experimental setup of Gault. For instance, we modeled the input size distribution as a step function, rather than a smooth power law as is normally done, to model the craters as discrete size classes. Three of the simulations had a grid size of 2000 2 , one had a grid size of 500 2 , and one had a grid size of 100 2 .
The lowest resolution runs were performed so that we could obtain complete counts of very small craters on a given surface without overburdening our human counter. We produced 10 4 m −2 craters on the Terra Alta surfaces and 10 5 m −2 on the Mare Nostrum surfaces. We also included one high resolution Terra Alta simulation where we produced twice as many craters in order to obtain a highly saturated surface.
The collection of simulated surface DEMs were given to co-author Fassett, who was asked to count craters on each of the surfaces and report their diameters and center locations on the grid. He identified a total of 744 craters. After receiving the list of identified craters for each simulated surface, we matched them with the list of craters produced by CTEM on that surface. The identified craters contained some error in size and position determination, and we had to adopt tolerance values for these two errors. Robbins et al. (2014) found that the typical position and diameter errors of professional crater counters was 5-10%. We adopted a somewhat more generous position tolerance of 30% relative to crater diameter and a diameter tolerance of 1/ √ 2D true < D reported < √ 2D true , which is the width of a standard R-plot bin (Crater Analysis Techniques Working Group, 1979). This diameter tolerance is also larger than the factor of 2 difference between the diameters of each crater size class, and therefore there is no ambiguity that a particular identified crater belongs to a particular diameter class.
Four craters were identified that did not correspond to any known crater produced by CTEM. We dub these types of detections "false positives," and because only a small number of these were reported, they are not expected to influence the overall results of our calibration and were discarded. The calibration test results are reported in Table 1 lists the features of each of the simulation types, including the size class of the smallest reported craters, the number of craters produced that were larger than the smallest reported crater size, and the total number of craters that were counted.
The goal of our analysis was to determine which measurable properties of a CTEM-generated crater correlated with countability. We identified two measures that both strongly correlate with the countability of a crater. The first we call the depth ratio measure, R d , which we define as: where h rb is the average rim-to-bowl height of the crater at the time of measurement, and h rb,0 is the original average rim-to-bowl height. These measurements are made relative to a reference plane that is determined by the average slope of the terrain just prior to crater emplacement. This allows us to accurately measure rim-to-bowl heights on craters that form on slopes (such as small craters that form on the inner walls of larger craters).
The cumulative fraction of craters below a given R d is plotted in Fig. 5a. We have separated out the subset of 111 test craters as red lines in this figure, showing that there was no substantial difference between the test craters and the complete set. Fig. 5a shows that there is a strong correlation between R d and the countability of a crater. This correlation can further be quantified by placing craters in bins of R d value, and then calculating the fraction of craters in each bin that were counted. We plot the fraction of craters counted in bins of R d , with a bin width of 0.1, in Fig. 5b. We then perform a least squares fit to the data to produce a model of the probability of being counted as a function of R d , which we call the depth ratio probability function p d . We found that a second order polynomial was a better fit than a linear function.
The best fit functional form of the depth ratio probability function is given as equation (18).
The value of p d for any given crater is a probability that a crater with its corresponding R d value will be counted by a human or not.
We also identified a second measure that also correlates strongly with the countability of the craters. The second measure we call the shape deviation measure, R σ , which we define as: Here, σ is the standard deviation of the average difference between the elevation of the crater pixels at the time of measurement and its original elevation after formation. Again, measurements are made relative to a reference plane that is at the average terrain slope prior to crater emplacement. For any given crater that occupies N pixels, σ is defined as: where h i is the elevation of pixel i and h i,0 is the elevation that pixel i had just after the crater formed, and The cumulative fraction of craters below a given R σ is plotted in Fig. 6a. Fig. 6a shows that there is a strong correlation between R σ and the countability of a crater. We can perform a similar analysis that we did for R d by placing craters in bins of R σ value, and then calculating the fraction of craters in each bin that were counted. We plot the fraction of craters counted in bins of R σ , with a bin width of 0.25, in Fig. 6b. We then perform a least squares fit to the data to produce a model of the probability of being counted as a function of R d , which we call the depth ratio probability function p d . We found that a linear function fit the data well. The best fit functional form of the shape deviation probability function is given as equation (22).
The two measures that we identified, R d and R σ , are correlated. Therefore we cannot treat them as independent probabilities. However, we have found that taking the minimum value of either p d or p σ , which we dub the "p-score," predicts the probability that any given crater is countable better than either measure alone. A histogram of the fraction of craters that were counted in bins of p-score (for all 560 craters in our calibration set) is plotted in Fig. 7.
The dashed line is the y(x) = x line, and would be the ideal case where the p-score exactly predicts the probability of any crater being counted. Our calibration parameters match the probability of counting quite well, closely adhering to the y(x) = x line.
We implement our calibration results into CTEM using the p-score. For every crater, we evaluate its p-score using by first calculating R d and R σ using equations (17) and (19), respectively. Then we calculate p d and p σ using equations (18) and (22), respectively. The minimum of the two becomes the p-score. Ideally, the p-score would be used as a probability that the crater is counted or not. Implementing the p-score as a pure probability function yielded complications to the code, because low-p craters could have their pscores artificially inflated if a new crater happened to form with a similar size and position as an old, degraded crater. This results in an artificial secular increase in the number of craters on a heavily cratered terrain over time. In order to prevent this, we use a p-score of 0.5 as a threshold for countability.
If the p-score is less than 0.5, then the crater is not counted, and the record of the crater's existence is obliterated from the grid so that it cannot be later confused for another similar-sized crater that might later happen to form close to the old crater's location.
We verified our calibrated counting using the "classic" example of a terrain in cratering equilibrium: Sinus Medii. Sinus Medii is a small mare deposit at the sub-Earth point of the Moon that was emplaced between 3.63-3.79 Gy ago (Hiesinger et al., 2010). This mare is often used as a case study in crater equilibrium (or saturation equilibrium), due to the "break" in the power law slope of the crater SFD at ∼ 100 m crater diameter (Gault, 1970;Marchi et al., 2012). We performed a test in CTEM designed to reproduce the cratering history of Sinus Medii. An impactor population with of the form N i ∝ D −3.25 i was used to generated the craters, and a "dry soil" model for the regolith material properties was assumed for the crater scaling relationship was used where µ = 0.41, K 1 = 0.24, andȲ = 0.18 MPa (Holsapple, 1993). terrain at this scale. We also do not see a secular increase in the small-crater abundance once it reaches the equilibrium line. The mismatch at the largest crater size is a resolution effect due to the limited size of our grid.

Impactor population
In this work we wish to test whether the highlands impacting population originated in the main asteroid belt, as suggested by Strom et al. (2005).
This constrains the impactor size frequency distribution, impactor density, and impact velocity distribution. We use the observed main asteroid belt size Our impact velocity distribution shown in Fig. 10 is derived from an Nbody simulation of the dynamical diffusion of main belt asteroids into the NEA region (Minton and Malhotra, 2010;Yue et al., 2013). The RMS velocity of this distribution is 18.3 km s −1 . We adopted 2.5 g cm −3 as the density of our impactors in line with typical densities of S-type asteroids (Britt et al., 2002), assuming that most lunar impactors derive from the S-dominated inner main asteroid belt (Bottke et al., 2006).
For each impact, CTEM draws an impactor size and velocity from the input distributions and determines a random location on the surface to place the impactor. The code also selects a random impact angle based upon the canonical formula whose derivation is summarized in Pierazzo and Melosh (2000): The normalized, integrated (cumulative) form of this equation is simply: where Eq. 24 is inverted to generated a random impact angle θ as a function of a randomly generated probability (P between 0 and 1) for each simulated impact in the model.

The effects of our velocity and angle distributions on the scaling relation-
ship is shown in Fig. 11. Here we plot the average number of craters produced in 500 CTEM simulations as a density map. The log of the average number of craters within a bin with a size of 0.1 × 0.1 log diameter (km) is plotted as the contour value. This also includes the effects of our basin scaling, which we use for D i > 35 km. Because the scaling parameters depend differently on size (mass) and velocity, there is a size dependence in how the various distributions shape the distribution of impactor sizes for a given crater size.

Lunar cratering simulations
Our goal is to test the hypothesis that impactors with the main belt asteroid SFD can produce the observed lunar highlands crater SFD. Our comparison data set is the catalog of all observed lunar craters with D > 20 km obtained using the Lunar Orbiter Laser Altimeter (LOLA) aboard the Lunar Reconnaissance Orbiter spacecraft . We studied this problem in two steps. First we performed a series of regional simulations designed to determine the minimum level of cratering (i.e. impactor flux exposure) needed to reproduced the observed abundance of craters found on the lunar highlands. Next we performed a series of global simulations designed to determine whether or not the best fit level of cratering obtained in the regional simulations is consistent with the observed number of large lunar basins. The regional constraint is a lower limit on the abundance of mid-sized craters, while the global constraint is an upper limit on the number of large basins. Because of the nature of the Monte Carlo technique, the "best fit" for any given constraint is 50%. That is, half the runs for a given level of cratering satisfy the constraint, while half do not. The total amount of cratering on a given run is parameterized by the quantity of expected impactors larger than 10 km in diameter,N pf (D i > 10 km), per unit area. The actual number of impactors, N i (D i > 10 km), will vary from run to run due to the nature of the Monte Carlo method, but over many runs the average number will equal the expected number. We step through values ofN pf (D i > 10 km) until we match the observed abundance of craters on the lunar highlands from our crater catalog. We make use of the relative crater density, or R-value, which is defined by (Crater Analysis Techniques Working Group, 1979) as: where N c is the number of craters within bins with boundaries (b 1 , b 2 ). We use standard R-plot bins where the bin boundaries are given as 1 km×(2) n/2 , and n is an integer (Crater Analysis Techniques Working Group, 1979). The geometric mean diameter,D c , of the bin is defined as: where D c,j are the individual craters in the bin. Alternatively, is used when the individual crater diameters are not available.
We will also make use of R-plots, which are plots of the log of R as a function of the log of D c , as a way of comparing size distributions of observed craters and modeled craters. For the regional simulation step, we use only the subset of craters in the LOLA catalog that occur on the lunar highlands, excluding Orientale and SPA basins. This region covers an area of 2×10 7 km 2 of the Moon, and is shown in Fig. 12 In particular, craters in the size range 90.5-128 km make for a very useful diagnostic for determining how many impacts are required to match the regional lunar highlands cratering record, because the relative crater density, or R-value, of lunar highlands craters in this size range is at a peak.
Once we have determined the probability that a given value ofN pf (D i > 10 km) will reproduce the observed crater density, we next performed a similar set of runs for a global lunar surface. In these runs, we define a basin constraint based on the observed abundances of lunar basins in our LOLA crater catalog, where we adopt the definition that a basin is any crater with D c > 300 km. The basins are observable, even if they are nearly completely filled with mare, such is the case with Imbrium. Recently measurements of the gravity signature of basins on the Moon obtained from the GRAIL spacecraft have revealed that the observed number of basins is a largely complete inventory of all basins that have formed on the Moon since the crust solidified . In particular, the largest impact structure on the Moon is the ∼ 2500 km South Pole-Aitken basin (SPA), which is a very prominent and deep topographic feature, and acts as a strong constraint on the total impact history of the Moon since the crust solidified. SPA appears to pre-date all other features on the Moon (Wilhelms et al., 1987), so it is somewhat ambiguous whether its formation is linked to the rest of the lunar basin formation process. Nevertheless, it is unlikely that any other basin as large or larger than SPA postdates its formation, and its nearest rival is the ∼ 1160 km Imbrium basin. Because SPA apparently pre-dates the entirety of the lunar highlands cratering record, we would be well justified in ignoring it in our study. Instead, we will adopt the conservative assumption that the SPA impactor originated from the same population that the rest of the highlands did, but by chance impacted near the beginning of the observable lunar cratering record. For our simulations, we adopt the constraint that we must produce no more 1 basin with D c > 1200 km (the size of Imbrium).
For the E-belt simulations our basin constraint is that we must produce no basins with D c > 1200 km.
For each set of conditions we performed 100-1000 CTEM simulations of the lunar surface. We tally the countable craters in each simulation, and bin them. However, we plot our results in a way that is somewhat unusual for crater counting. Usually crater counts are reported using error bars that are scaled by ±N 1/2 in size to reflect the assumption that the variability in the number of craters for a given amount of flux is Poissondistributed (Crater Analysis Techniques Working Group, 1979). However, we make no a priori assumption as to how crater variability is characterized.
While impact events are well described with Poisson statistics, because the area of a cratered surface is finite and craters may obliterate one another, craters are not strictly Poisson-distributed. Large craters are more effective obliteration agents than small craters, and as craters approach a size comparable to the counting area, Poisson statistics are a poor model for the intrinsic variability of craters. By performing a large number of runs, we can directly estimate the variability of crater counts directly without having to assume that they are Poisson-distributed.
To report our CTEM-derived results for multiple runs with the identical parameters (except for the random number generator seed), we make use of "box and whisker plots." This style of plot contains a box within an error bar. The box is drawn to span the 25% of data points above the median and the 25% of data points below the median. The value of the median is drawn as a horizontal line within the box. The error bars enclose 99% of the data. Because we report our model statistics, we have no need to estimate statistical variations in the observed data set. Therefore our observational data set (the catalog of LOLA crater of D c > 20 km) is plotted simply as points with no error bars. An example of data plotted in this fashion is given in Fig. 13.
Our regional CTEM simulations were performed on a 2000 × 2000 pixel grid at a resolution of 2.24 km px −1 . Our catalog of highlands craters includes craters on an area of 2 × 10 7 km 2 of the lunar surface that excludes SPA and Orientale (see Fig. 12). This catalog contains no craters larger than 628 km in diameter. In these simulations we deliberately restricted our cratering model to produce only craters less than this maximum size. The craters in this region were not greatly effected by mare emplacement, and therefore more closely reflect the crater density of the pre-mare lunar surface than the global catalog. We bombarded a simulated lunar surface with an amount of craters equal toN pf (D i > 10 km) = 3.5 × 10 −6 km −2 -6 × 10 −6 km −2 .
We found that the bin that spans D c = 90.5-128 km to be a useful diagnostic of the success of cratering of a given terrain, as it has the highest value of R in the regional data set. The total number of observed craters in LOLA lunar highlands regional data for this bin is 198, and we set as our constraint that our simulated surface must also have this many craters in this size range to be considered a success. Fig. 13 shows results of four sets of CTEM simulations in R-plot form. ForN pf (D i > 10 km) < 4.0 × 10 −6 km −2 , none of the runs produced the observed number of craters in the 90.5-128 km bin. ForN pf (D i > 10 km) = 5.0 × 10 −6 km −2 , 55.6 ± 2.4% of runs satisfied the regional mid-sized crater constraint, and therefore is close to the best fit cratering abundance for the main belt size distribution. We plot these statistical results as red points in Fig. 14.

Global lunar surface runs with the basin constraint
Our global CTEM simulations were performed on a 2000×2000 pixel grid at a resolution of 3 km px −1 . We bombarded a simulated global lunar surface with an amount of craters equal toN pf (D i > 10 km) = 5 × 10 −7 km −2 -6 × 10 −6 km −2 . We set a constraint that there should be no more than 1 basin with D c > 1200 km and ended the run if a crater this large was generated. Fig. 14 shows the fraction of simulations that fit our basin constraints as a function ofN pf (D i > 10 km). As we saw earlier only runs withN pf (D i > 10 km) > 4 × 10 −6 km −2 satisfied the regional constraint, with 5 × 10 −6 km −2 being close to a best fit value. However, for N pf (D i > 10 km) = 4 × 10 −6 km −2 , only 0.8 ± 0.28% of runs satisfied the basin constraint. Because we end the runs at the moment that a crater larger than 1200 km in diameter is generated, we can determine the mean value ofN pf (D i > 10 km) at the time that the first constraint-violating basin formed. We find that the global constraint is violated on average at N pf (D i > 10 km) = 9.2 × 10 −7 km −2 . This is approximately a factor of 5 less cratering than the best fit value ofN pf determined using the regional constraint.
We then compared the average global cratering record for those runs that satisfied the basin constraints and compared the crater size distribution of smaller craters. Our two constraints were only very rarely satisfied forN pf (D i > 10 km) = 5 × 10 −6 km −2 , where 8 out of 1000 global simulations satisfied the basin constraint. We plot the resulting crater R-plot distributions for the runs withN pf (D i > 10 km) = 1 × 10 −6 km −2 and N pf (D i > 10 km) = 5 × 10 −6 km −2 in Fig. 15. AlthoughN pf (D i > 10 km) = 1×10 −6 km −2 is very near the best fit value of cratering that satisfies the basin constraint, none of the runs had enough craters to satisfy the regional constraint. As we showed in our regional simulations in Section 3.1, we needed N pf (D i > 10 km) > 4 × 10 −6 km −2 in order to match the observed abundance of craters in the 90.5 km-128 km size range. The global constraint, by contrast, is better matched atN pf (D i > 10 km) = 9.2 × 10 −7 km −2 .

Runs that test the E-belt hypothesis
The region inward of ∼ 2 AU is currently unstable due to the presence of the ν 6 secular resonance (Williams and Faulkner, 1981). If the giant planets formed in a more compact configuration, the ν 6 would have been farther from the Sun than its present location (Minton and Malhotra, 2011;Agnor and Lin, 2012), and the region between the inner main asteroid belt and Mars could have been filled with asteroids. Under the E-belt hypothesis, described in , the arrival of the ν 6 resonance after giant planet migration would have been responsible for destabilizing this innermost portion of the asteroid belt (the E-belt), and the impactors associated with the LHB would have primarily come from this region. The E-belt, plus a small contribution from the main belt, could only supply enough large impactors to produce the sequence of basins beginning with Nectaris. This population was assumed to have had a similar size-frequency distribution as the current inner main asteroid belt, however the impact velocity of this population was somewhat higher, with a median velocity of 21 km s −1 instead of the 18.3 km s −1 that we used in the previous sections.
Based on LOLA crater counts, crater density on Nectaris basin is N(64) = 17 ± 5 and N(20) = 135 ± 14 (Fassett et al., 2012). Here N(D) refers to the N >Dc per 10 6 km 2 surface area. Fassett et al. (2012) also report 14 basins younger than Nectaris, plus Nectaris itself. The largest crater younger then Nectaris is the D c ≃ 1200 km Imbrium basin. We performed a series of regional runs similar to those described in Section 3.1, as well as a series of global runs similar to those described in Section 3.2. The results are shown in Fig. 16. The fraction runs for a given value ofN pf (D i > 10 km) that produced N(20) and N(64) densities within the range determined by Fassett et al. (2012) for Nectaris are plotted as green triangles and red squares respectively.
For clarification, because here we are using crater densities as our observational constraint, rather than total number of craters as in Sections 3.1 and 3.2, we do not use 50% match as our criteria for "best fit." Here we simply use the reported uncertainty ranges of the crater densities to determine whether model runs at a given number of impactors are a good fit or not.
The N(64) densities suggest values ofN pf (D i > 10 km) = 1.5-2.25 × 10 −6 km −2 in order to reach the observed crater densities on Nectaris. The N(20) densities suggest somewhat higher values ofN pf (D i > 10 km) = 2.25-2.75 × 10 −6 km −2 , however this could simply imply that the dataset is either overabundant in ∼ 20 km craters, or the slope our SFD for craters with D c < 100 km is too shallow. To avoid any ambiguity in the D c ≃ 20 km crater counts we adopt the N(64) densities as our constraint. For our E-belt basin constraint, we require that runs produce no basins with D c > 1200 km.

Discussion & Conclusions
We defined two constraints that must be satisfied simultaneously in order for the main asteroid belt size-frequency distribution to produced the observed lunar highlands crater size-frequency distribution. Our two constraints are a regional one (we must produce at least the total number of observed D c ≃ 100 km craters on the highlands), and a global one (we must produce no more than the observed number of basins and none larger than the largest observed basin). Our results indicate a very low probability in matching both of these constraints with an impactor population resembling the modern main asteroid belt. The modern main asteroid belt SFD is therefore a poor model for producing the observed lunar highlands crater population. This is due to the relative abundance within the main asteroid belt of objects that would produce lunar basins larger than Imbrium (and some larger than South Pole-Aitken) if they collided with the Moon, as compared with main belt objects that would produce D c ≃ 100 km craters.
From Fig. 14 we show that the regional constraint is not met in any of our runs atN pf (D i > 10 km) < 4 × 10 −6 km −2 . Our best fit value to the regional constraint isN pf (D i > 10 km) = 5 × 10 −6 km −2 . However, at that level of cratering, very few of our runs satisfy the global constraint. They nearly all produce at least one basin larger than Imbrium. The average value ofN pf at the moment that a basin larger than 1200 km is produced is N pf (D i > 10 km) = 9.2 × 10 −7 km −2 .
A simple numerical exercise can be used to demonstrate why this is so by comparing the relative abundance of impactors that produce craters of these sizes within the main asteroid belt. In Fig. 17 we plot the crater diameter vs. impactor diameter for craters produced in a CTEM global lunar surface simulation atN pf = 5 × 10 −6 km −2 for two different size ranges of crater. We can see from this figure that to make a crater in the diameter range 90.5-128 km requires an impactor of ∼ 3-15 km, with most produced by impactors of ∼ 5-6 km in diameter. To make a basin larger than the ∼ 1200 km Imbrium basin requires an impactor of ∼ 30-300 km. On average an impactor with D i < 70 km created at least one megabasin in these runs.
In our LOLA crater catalog for the lunar highlands, we have 296 craters larger than 90.5 km in diameter, and globally there is only one lunar basin larger than Imbrium basin, which is the SPA basin. If we extrapolate the number of craters in the 2.0 × 10 7 km 2 area of highlands to reflect what the 3.8 × 10 7 km 2 global pre-mare lunar surface may have experienced we have ∼ 560 craters with D c > 90.5 km. Not every crater produced is observed due to the variety of crater erasure mechanisms described in Sec. 2.3. In our regional runs atN pf (D i > 10 km) = 5 × 10 −6 km −2 , we calculate that 89 ± 2.4% of all craters produced with D c > 90.5 are observable. This gives us a ratio of N >90.5 km /N >1200 km ≃ 630. From the main asteroid belt model size-frequency distribution plotted in Fig. 9, we can estimate that the corresponding ratio of N >5.5 km /N >70 km impactors is ∼ 100. Therefore it is highly unlikely that the lunar surface could have been bombarded by enough main belt asteroids to produce the observed number of D > 90.5 km craters without producing many basins larger than Imbrium. We quantified this probability in Fig. 14, which we find to be 0.8 ± 0.28%. We also showed in Section 3.3 that the likelihood that the Nectarian crater densities could be produced while not producing any basins larger than Imbrium was 16%. We should note that we obtained this results assuming that the Ebelt SFD was similar to the average main belt, rather than the inner main asteroid belt. However, the inner main asteroid belt may even more topheavy than the main belt as a whole. From Masiero et al. (2011), for the inner main belt N >5.5 km /N >70 km = 74. Using the inner main belt as the impactor population would therefore reduce the probability that the regional and global constraints could be met simultaneously.
One possible solution to this problem could be that our scaling relationships between impactor size and final crater size are incorrect. If, for instance, the impactors that generated the D c ≃ 100 km craters were smaller (more numerous) than the D i ≃ 5 km that we have estimated, or if the impactors that generated the D c > 1200 km megabasins were larger (less numerous) than the D i 70 km that we have estimated, then it is possible that the main belt SFD could produce the required ratio N >90.5 km /N >1200 km = 630. Using our SFD, we estimate that in order to achieve this ratio, the impactors that generated the D c ≃ 100 km craters would have to be D i = 2.2 km, while keeping the basins scaling relationship unchanged. Alternatively, the megabasin impactors would have to be D i = 155 km, while keeping the 100 km scaling relationship unchanged. These values can be parameterized as a gravity-scaled impactor size, π 2 = 1.61gD i /v 2 i , and dimensionless transient crater diameter π D = D t (ρ t /m i ) 1/3 (Schmidt and Housen, 1987;Melosh, 1989;Wünnemann and Ivanov, 2003;Wünnemann et al., 2006;Elbeshausen et al., 2009). We plot these alternate scaling results along with the experimentally and numerically derived scaling relationship of common materials in Fig. 18.
Both of these alternative scaling relationships fall well outside of experimentally and numerically determined scaling relationships. Therefore it is unlikely that the scaling relationships for these craters are uncertain enough to explain our results.
One result of our study is the revelation that the small body population that cratered the lunar highlands was more depleted in D > 70 km objects than the main asteroid belt (or, equivalently, more abundant in D > 5 km objects). The relative abundances of D > 70 km asteroids has been used as evidence that the initial size of planetesimals in the protoplanetary disk was this large . A similar feature in the size distribution is suspected for the Kuiper belt (Fraser, 2009;Shankman et al., 2013). In the case of the Kuiper belt, the ratio of these large objects to smaller objects is still quite uncertain, though work is ongoing to address this uncertainty Minton et al., 2012).
There is certainly a great deal of similarity between the shape of the main belt size distribution and that of the population of impactors that produced the lunar highlands. For instance, using our scaling relationship parameters of µ = 0.55 and K 1 = 0.22, which are within the range expected for impacts rocky objects onto rocky targets (Holsapple, 1993), and our velocity distribution derived from main belt asteroids, the resulting craters both exhibit a peak on an R-plot near 100 km and a trough near 400-500 km (see Fig. 13, for instance). Caution must be used in over interpreting this similarity. The critical specific energy required to disrupt a small body, Q * D , changes from strength-dominated for small objects to gravity dominated for large ones (Benz and Asphaug, 1999). Bodies become weaker per unit mass as they become larger in the strength-dominated regime, but stronger as they become larger in the gravity-dominated regime. This change in the size dependence of strength sets up standing waves in a collisionally evolved population, and the amplitude, wavelength, and phase of the waves depend on the shape of the strength law and the mutual collision velocities of the bodies in the population (O'Brien and Greenberg, 2003). Therefore, the similarities seen between the lunar highlands impactors and the main asteroid belt (when megabasins are ignored) may simply reflect that these two populations had similar values of material strength (i.e. they were made of rock) and mutual impact velocities (i.e. they collided with each other while orbiting in the inner solar system). Nevertheless, we conclude that it is unlikely that the lunar highlands were bombarded by a population with a size-frequency distribution identical to that of the modern main asteroid belt.  crater, but only partially cut the smaller offset crater. CTEM no longer counts the 1 km centered crater, but still counts the partially cut offset crater. The grid was 1000 2 px, and the resolution was 6 m px −1 . Lunar gravity was assumed as g = 1.62 m s −2 and target density was ρ t = 2.7 g cm −3 . The grid was 1000 2 px, and the resolution was 6 m px −1 . Lunar gravity was assumed as g = 1.62 m s −2 , target density was ρ t = 2.7 g cm −3 , and projectile density was ρ p = 2.5 g cm −3 .   (Hiesinger et al., 2010). This mare is often used as a case study in crater equilibrium (or saturation equilibrium), due to the "break" in the power law slope of the crater SFD at ∼ 100 m crater diameter (Gault, 1970;Marchi et al., 2012). i was used to generated the craters, and a "dry soil" model for the regolith material properties was assumed for the crater scaling relationship was used where µ = 0.41, K 1 = 0.24, andȲ = 0.18 MPa (Holsapple, 1993).  Figure 11: Example of the crater-size scaling relationship (from impactor to crater diameter) used in our models (Holsapple, 1993). The scaling relationship is for the Lunar surface with a surface gravity of 1.62 m s −2 , including the effects of our impact velocity distribution shown in Fig 10, and impact angle distribution given by Eq. 24. We used K 1 = 0.22 and µ = 0.55 as our scaling parameters (see Eq. 6). The contours represent the average numbers of objects within bins of 0.1 × 0.1 log diameter (km) per simulation using 500 simulations. Figure 12: The color region is the terrain used in our regional simulations. This is the lunar highlands excluding SPA basin. 4.5 × 10 −6 km −2 , 5.0 × 10 −6 km −2 , and 5.5 × 10 −6 km −2 . The CTEM data is plotted as a box and whisker plot. The box is drawn to span the 25% of data points above the median and the 25% of data points below the median. The value of the median is drawn as a horizontal line within the box. The error bars enclose 99% of the data. Global: Basin constraint Regional: N 90.5-128 km constraint Figure 14: Fraction of CTEM runs that two of our constraints as a function of number of impactors,N pf (D i > 10 km). The circles are the global runs that satisfy the basin constraint: No more than 1 basin with D > 1200 km. The squares are the regional runs that satisfy the N 90.5−128 km constraint: The bin spanning 90.5-128 km must contain at least 198 craters, which is the number observed in the LOLA lunar highlands data set.

Tables
The impactors were drawn from the main asteroid belt size distribution. 10 km) = 5 × 10 −6 km −2 (right). In these runs, the simulations were ended when a basin larger than 1200 km diameter was generated, and therefore the average value of N pf (D i > 10 km) = 9.2 × 10 −7 km −2 . The CTEM data is plotted as a box and whisker plot. The box is drawn to span the 25% of data points above the median and the 25% of data points below the median. The value of the median is drawn as a horizontal line within the box. The error bars enclose 99% of the data. The left panel shows the global crater counts for the set of runs near the best fit value ofN pf for the global basin constraint.
The right panel shows the global crater counts for the set of runs near the best fit value ofN pf for the regional 90.5 < D c < 128 km constraint. Global: Basin constraint Regional: N(64) constraint Regional: N(20) constraint Figure 16: Fraction of CTEM runs that two of our E-belt constraints as a function of N pf (D i > 10 km). The black circles are the global runs that satisfy the basin constraints: No basins with D c > 1200 km. The red squares are the regional runs that satisfy the constraint that N (64) = 17 ± 5, and the green triangles are the runs that satisfy the constraint that N (20) = 135 ± 14, based on crater densities on Nectaris (Fassett et al., 2012). The impactors were drawn from the main asteroid belt size distribution. The right panel shows the distribution of impactors of a given size (in bins of 10 km) that produce craters with D c > 1200 km. The plot shows the effect of all parameters that effect the scaling law in a given simulation, including impactor size, velocity, and impact angle, and in the case of impactors with D i > 35 km, our two crustal thermal profiles. ing relationships (Melosh, 1989;Wünnemann and Ivanov, 2003;Wünnemann et al., 2006;Elbeshausen et al., 2009).