The Moon's surface has experienced eons of bombardment by asteroid impacts. This bombardment is thought to have thoroughly shattered the lunar crust, the outermost layer of the Moon, to depths of 20 km or more. Despite this, the processes of fracturing and fragmentation during an asteroid impact have received little attention though. In this study, we simulated asteroid impacts and tracked how impacts shatter the lunar crust. We found that impacts break the lunar crust into roughly meter‐sized blocks. Impacts can break up the lunar crust down to approximately 20 km, and objects greater than 1 km in diameter are the most efficient at fragmenting the crust. This implies that the fragmented bedrock can be almost exclusively created by large‐scale high‐speed impacts and the lunar crust was thoroughly fractured early in its history. Similarly, impacts likely shattered the crust of the ancient Earth and Mars. Understanding how thoroughly fractured ancient planetary crusts are may help us determine if and when these environments became habitable.

The megaregolith of the Moon is the upper region of the crust, which has been extensively fractured by intense impact bombardment. Little is known about the formation and evolution of the lunar megaregolith. Here we implement the Grady‐Kipp model for dynamic fragmentation into the iSALE shock physics code. This implementation allows us to directly simulate tensile in situ impact fragmentation of the lunar crust. We find that fragment sizes are weakly dependent on impactor size and impact velocity. For impactors 1 km in diameter or smaller, a hemispherical zone centered on the point of impact contains meter‐scale fragments. For an impactor 1 km in diameter this zone extends to depths of 20 km. At larger impactor sizes, overburden pressure inhibits fragmentation and only a near‐surface zone is fragmented. For a 10‐km‐diameter impactor, this surface zone extends to a depth of ~20 km and lateral distances ~300 km from the point of impact. This suggests that impactors from 1 to 10 km in diameter can efficiently fragment the entire lunar crust to depths of ~20 km, implying that much of the modern day megaregolith can be created by single impacts rather than by multiple large impact events.

1 Introduction Impact cratering is a ubiquitous process that dominates the evolution of ancient planetary crusts (Melosh, 1989). During a hypervelocity impact a shockwave propagates into the projectile and target and is followed by a rarefaction wave that propagates from free surfaces. This succession of shock and rarefaction sets up the flow that excavates the crater (Melosh, 1984). The rarefaction can produce strong tensile stresses and high strain rates, triggering dynamic fragmentation (Melosh et al., 1992). The volume of this fragmented material is much larger than the volume excavated by the impact (Melosh et al., 1992). The role of in situ fragmentation during crater formation, however, has received little attention. As the Moon's ancient highlands are the best preserved and best studied example of an ancient planetary crust, our work focuses on in situ impact fragmentation of the lunar crust, making it the first study of its kind. The thoroughly fragmented upper 10 km or more of the lunar crust is known as the megaregolith (Hartmann, 1973; Heiken et al., 1991), which can be subdivided into three layers. The uppermost stratum is the regolith layer, which is meters thick and is composed of fine‐grained heavily reworked material (Bart et al., 2011; Papike et al., 1982). Underneath the regolith layer in the lunar highlands lies the “large‐scale ejecta layer,” which is composed of material ejected during the formation of large basins and is approximately a few kilometers thick (Heiken et al., 1991). Below this is the thickest layer, composed of material fractured or fragmented by impacts but not ejected (Heiken et al., 1991). Here we refer to this layer as the “in situ fragmentation layer.” While the near‐surface regolith has been directly sampled by the Apollo missions, much of the deeper structure has been determined by early seismic studies from the Apollo missions (Heiken et al., 1991; Latham et al., 1973; Toksöz et al., 1972). It should be acknowledged here that the megaregolith, as defined by Hartmann (1973), is the fractured crust below the regolith layer, that it is not formed in the same way the regolith is formed, and is distinct from the regolith, in that it is not just reworked material (i.e., the ejecta of basins). This has been shown by Gravity Recovery and Interior Laboratory (GRAIL), as it found that the fragmented crust is too thick for it to be composed mostly of ejecta and that the Moon is not saturated in basins (Neumann et al., 2015). As megaregolith is produced by impacts and a planetary crust becomes more fractured and porous, its thermal conductivity is reduced and permeability increased, and magma ascent becomes less favorable (Jaeger et al., 2007; Melosh, 2011; Warren & Rasmussen, 1987). Thus, the timing of formation and evolution of the megaregolith is important for thermal, magmatic, and, where applicable, the hydrological evolution of planetary bodies (Clifford, 1993; Haack et al., 1990; Shearer et al., 2006; Zhang et al., 2013). Previous studies of megaregolith formation have focused on large‐scale basin ejecta (i.e., the large‐scale ejecta layer), which makes up the upper few kilometers of the lunar megaregolith (Hörz, 1977; Rolf et al., 2017; Short & Foreman, 1972). The state of the bedrock below craters (i.e., the in situ fragmented layer), however, has been largely unstudied until now (Hartmann, 2003). The 2011–2012 GRAIL mission provided an unprecedented look into the interior structure of the Moon (Zuber et al., 2013). GRAIL revealed that at least the upper few kilometers of the lunar crust has an average crustal porosity of 12% (Wieczorek et al., 2013). Although this porosity decreases with depth, typical porosity is ~8% and ~3–4% at depths of 10 and 20 km, respectively (Besserer et al., 2014). The GRAIL results suggested that by focusing on the upper few kilometers of large‐scale basin ejecta we are only considering the surface of the megaregolith. Appreciable porosity present to the base of the lunar crust (Besserer et al., 2014; Wieczorek et al., 2013) is likely the result of in situ fragmentation, but this process has received little attention. We aim to show that the process of dynamic fragmentation is key to the development of the megaregolith. During dynamic fragmentation, tensile stresses cause multitudes of cracks to grow and link up, producing independent fragments. The number of cracks that are active depends on material properties and accumulated strain (Grady & Kipp, 1987). Ultimately, the typical fragment sizes produced during dynamic fragmentation depend on strain rates and material properties. This is unlike quasi‐static fragmentation, which operates on the growth of the weakest available crack in a material (Melosh et al., 1992). To describe the process of dynamic fragmentation, Grady and Kipp (1987) developed a model that tracks the cumulative effect of crack growth due to tensile stresses. The implementation of the Grady‐Kipp fragmentation model into the iSALE shock physics code grants us the ability to quantify the contribution of tensile forces in fragmentation, and, in particular, it allows us to explore in situ fragmentation and its role in the development of the lunar megaregolith for the first time. In this work, we focus on the creation of in situ fragmented lunar crust to better understand the structure and formation of the megaregolith. We outline the parameters and tools we use to explore impact fragmentation below, with details of the Grady‐Kipp fragmentation model and its validation given in Appendix A. We apply these methods to a range of situations, the results of which are given in section 3. In section 4, we discuss the results and how they relate to the overall process of impact fragmentation, and how it relates to the creation of lunar megaregolith. We make our concluding statements in section 5.

2 Methods To simulate the development of deep megaregolith through impacts, we use the multimaterial, multirheology shock physics hydrocode iSALE‐2D (Collins et al., 2004; Ivanov et al., 1997; Melosh et al., 1992; Wünnemann et al., 2006), which is based on the SALE hydrocode (Amsden et al., 1980), with an implementation of the Grady‐Kipp fragmentation model (Grady & Kipp, 1980). For a description of our implementation and testing of the Grady‐Kipp fragmentation model see Appendix A. In the Grady‐Kipp model, k and m are material parameters that describe the number of active flaws per unit volume, N = kεm, as a function of strain ε. We set k = 1036 and m = 9.5, which provides a good fit to experimental data (Takagi et al., 1984) for basaltic targets (Appendix A). The use of a 2‐D hydrocode instead of a 3‐D hydrocode is more than sufficient for the simulations outlined below, as planetary curvature is not important, and through late stage equivalence (the concept that far from the point of impact, the dynamics of an impact with the same coupling parameter are similar; the coupling parameter depends on impactor size and velocity; e.g., a small fast impactor can have the same coupling parameter as a large slow impactor; Dienes & Walsh, 1970; Holsapple, 1993), only at very oblique impact angles would the difference between 2‐D and 3‐D become apparent (Pierazzo & Melosh, 2000). Here we simulate vertical impacts of spherical basaltic impactors on flat, nonporous, basaltic targets. We use the ANEOS (Analytic equations of state) equation of state (EOS) for basalt (Pierazzo et al., 2005), and the rock‐like strength and damage model of Collins et al. (2004) for both the impactor and the target. The basalt EOS has a similar reference density to the grain density of lunar anorthosite and should be a reasonable choice to simulate the behavior of intact lunar crust (Kiefer et al., 2012; Pierazzo et al., 2005; Taylor & Wieczorek, 2014). Material strength parameters come from Potter et al. (2012) fits to gabbro rock strength (Table 1). It is important to mention here that our current implementation of the Grady‐Kipp model into iSALE is only stable for one material type for each simulation. Table 1. iSALE Material Input Parameters Parameter Value Reference EOS ANEOS basalt 1 Poisson's ratio 0.25 2 Melting temperaturea 1513 K 2 Thermal softening coefficienta 1.2 2 Simon A parameterb 1,840 MPa 2 Simon C parameterb 7.27 2 Frictional coefficient (damaged)a 0.71 2 Frictional coefficient (undamaged)a 1.1 2 Cohesive strength (damaged)a 0.01 MPa 2 Cohesive strength (undamaged)a 31.9 MPa 2 Strength limit (damaged)a 2.49 GPa 2 Strength limit (undamaged)a 2.49 GPa 2 k 1036 Appendix A m 9.5 Appendix A Each of the vertical impact runs was performed with a constant gravitational acceleration set to 1.62 m/s2 to match that of the Moon and an impactor velocity of 15 km/s (a reasonable estimate for lunar collisions; Le Feuvre & Wieczorek, 2011; Yue et al., 2013), unless otherwise specified. The runs were kept at a constant resolution of 10 cells per projectile radius (cppr; Elbeshausen et al., 2009; Pierazzo et al., 2008; Wünnemann et al., 2006, 2008), as the computational expense of higher‐resolution simulations would be prohibitive. This is primarily due to the large high‐resolution zones of our simulations, where, unlike many simulations, our high‐resolution zone extends well past the crater rim and covers the entirety of the fractured regions. To examine the effects of impactor size on the development of fragments and megaregolith, the impactors' diameters were varied from 10 m to 10 km in logarithmic steps (10 m, 100 m, 1 km, and 10 km). The details of these simulation setups (i.e., grid size) is given in Table 2. We additionally examined the effects of overburden pressure, different impact velocities, and varying k. It should also be noted that we did not include a dilatancy model in any simulations presented here. No acoustic fluidization model was included in these simulations because we end our simulations before the modification stage of crater formation. Table 2. Simulation Parameters Impactor diameter (km) cppr Grid size in X (from symmetry axis) and Y (down from surface) directions X (cells) X (km) Y (cells) Y (km) 0.01 10 1,800 0.9 2,600 1.2 0.1 10 1,500 7.5 2,200 10.5 1 10 1,200 60 1,200 55 10 10 950 475 500 200

4 Discussion Our current implementation of Grady‐Kipp fragmentation model, in addition to currently working with just one material type and being relatively computationally expensive due to the very large high‐resolution zones, only tracks accumulation of fracture area when material is in tension. Thus, fragment size for material failing in compression or shear failure are inaccurate. This results in a significant volume of damaged target material not being tracked as it fragments due to shear stresses. Therefore, in all simulations discussed above, there are zones of large fragment sizes adjacent to the transient crater, which should be viewed as an absolute upper bound of fragment sizes (for more information, see Appendix A). It is possible that impacts with differing impactor sizes and velocities, but the same coupling parameter, will result in different fragment sizes in the regions close to the point of impact. Determining if such a dependence exists may be important for understanding the dynamics of crater collapse and its apparent sensitivity to impact velocity (Silber et al., 2017). An accurate calculation of fragment sizes for this region will require estimates of how fracture area accumulates during shear failure. Some preliminary work has been done (Melosh et al., 2017), assuming a fixed shear strain from the shock rather than a constant tensile strain rate, which led Melosh et al. (2017) to obtain a vastly lower k value than we have presented here. It is not clear this assumption is valid when material exceeds the Hugoniot elastic limit, as plastic failure could be ductile in some cases or fractures may be healed before material is put into tension by the rarefaction. A potential work‐around is to let shear damage and tensile damage accumulate separately, which could provide a better estimate of fragment sizes in the shear zone. Regardless, higher‐resolution simulations with the implementation of shear fragmentation calculations will provide more insight into the shear zones, mixed shear and tension zones, and the interaction of both shear and tensile fragmentation. We have not included a dilatancy model in the majority of our simulations, where dilatancy is the volume change experienced by material undergoing shear deformation. A dilatancy model is necessary to account for the creation of porosity by impacts, as was shown in previous work on the formation of lunar craters (Collins, 2014; Milbury et al., 2015). Inclusion of the dilatancy model of Collins (2014) in the 1‐km impactor simulation (with all parameters the same as the original 1‐km simulation but with the dilatancy model included) has resulted in no appreciable difference in fragment sizes or extent of the HTFZ, with a porosity ≤0.1% at 7 km deep down to the end of the HTFZ. This simulation likely represents a minimum estimate for porosity produced by such an impact since it neglects the creation of porosity during tensile failure if shear is limited. The inclusion of preimpact porosity, however, would affect the results of our simulations by causing the fragmentation field to be smaller. Impacts into the porous upper crust of the Moon have been implicated in the finding of positive free air gravity anomalies in midsized lunar craters (Milbury et al., 2015). However, we leave the implementation of porosity for a future study. Our results have implications for the creation and the current understanding of the lunar megaregolith. Classically, the megaregolith is viewed as the upper ~10 km of the lunar crust (Heiken et al., 1991); however, GRAIL has revealed that the lunar crust must be more fractured and porous than previously thought (Besserer et al., 2014; Wieczorek et al., 2013). Our two largest simulations, 1‐ and 10‐km impactors, match what we interpret to be the extent of fracturing predicted for the lunar megaregolith. Impactors larger than 10 km would show a similar pattern to what was observed with the 10‐km run, with the majority of the fragmentation occurring near the surface far away from the contact site. The largest two impactors (1 km and 10 km) fragment the upper 20 km of crust, and over time we expect that the crust of the lunar highlands will be completely fragmented down to this depth and will reach a fragmentation saturation threshold. The smaller impactors will further fragment the upper layers of the crust until the impact induced strains are insufficient to break apart already fragmented rock at depth. It also appears that the two largest impactor sizes studied here, the 1‐ and 10‐km impactors, are able to account for a substantial amount of in situ fragmentation, which implies that multiple impacts may not be necessary for the creation of the deep lying megaregolith and that in fact a single large impactor can locally create much of the modern megaregolith. However, we still need porosity estimates to more definitively probe the structure of the upper lunar crust. Porosity estimates, when combined with a Monte Carlo code that tracks how terrains evolve due to impacts, such as the Crater Terrain Evolution Model (CTEM) of Minton et al. (2015), can provide insight into how the lunar crust specifically evolved with time.

5 Conclusion We have implemented Grady‐Kipp fragmentation in iSALE‐2D allowing us to simulate the process of in situ impact fragmentation. Our study shows that impact cratering on the Moon is an efficient process for fragmenting the lunar crust to depths of ~20 km. Overburden pressure impedes fragmentation at even greater depths. The fragment sizes within the HTFZ are generally on the order of 0.1 to 1 m, with increasing sizes further away from the point of impact reaching sizes on the order of 10 m. For all impactor sizes studied here that produced HTFZs, the average fragment sizes within them are on the order of 1 m. However, it is important to note that the average fragment sizes within the HTFZs do appear to change with increasing impactor size, demonstrating a slight dependence upon impact energy. Even moderately sized craters, ~25 km in diameter, can form megaregolith down to depths exceeding 20 km, which explains the ubiquity of the megaregolith in the lunar highlands. Basin‐forming impacts are capable of fragmenting the lunar crust to a depth of ~20 km, far from the point of impact, and may be a dominant agent of megaregolith formation.

Acknowledgments We thank Gareth Collins for the helpful feedback and advice regarding the implementation of our code. We gratefully acknowledge the developers of iSALE‐2D, the simulation code used in our research, including G. Collins, K. Wünnemann, T. Davison, B. Ivanov, and D. Elbeshausen. We also acknowledge Tom Davison, the developer of the pysaleplot tool, which was used to create many of the figures given in the work. S. E. W. acknowledges funding from the Rhode Island Space Grant Consortium. Thank you to the reviewers, including Auriol Rae, for their helpful criticism and insight. All data associated with this study are listed in tables and shown in figures. Scientists interested in using or developing iSALE should see https://www.isale‐code.de/redmine/projects/isale/wiki/Terms_of_use for a description of application requirements. Simulation inputs and outputs are available on Harvard Dataverse (https://doi.org/10.7910/DVN/ENCZGA).

Appendix A: Details of Grady‐Kipp Fragmentation and Implementation 1980 1992 (A1) (A2) Following Grady and Kipp () and Melosh et al. (), we track the evolution of tensile damage using the following equation:where c g is the crack growth speed taken to be 0.4 times the longitudinal sound speed (Benz & Asphaug, 1994 m and k are Weibull distribution parameters, such that at a given strain ε the number of active flaws is N = kεm . Following Benz and Asphaug ( 1994 σ max , to calculate the strain in ε (A3) D is material damage and E is the Young's modulus of undamaged material, E = 9Kμ/(3K + μ), where K is the bulk modulus and μ is the shear modulus. Damage is only allowed to accumulate according to equation (Melosh et al., 1992 1994 (A4) V is the volume of material in our simulation and p(i) is a random number between 1 and N . N = V/ min (v cell ) where min(v cell ) is the volume of the smallest cell in a simulation. Because our simulations are axisymmetric, cell volume increases with distance from the symmetry axis. When the volume of a cell, v cell , exceeds min(v cell ) , we calculate ε m for that cell m = v cell / min (v cell ) times. We then take the minimum value of ε m calculated to represent the weakest flaw in the cell's volume. As Benz and Asphaug ( 1994 V, is active at , but locally the minimum failure strain does not depend on the size of the target or mesh. Figure Hereis the crack growth speed taken to be 0.4 times the longitudinal sound speed (Benz & Asphaug,),andare Weibull distribution parameters, such that at a given strainthe number of active flaws is. Following Benz and Asphaug () we use the maximum (tensile) principal component of stress,, to calculate the strain inwhereis material damage andis the Young's modulus of undamaged material,whereis the bulk modulus andis the shear modulus. Damage is only allowed to accumulate according to equation A1 if some minimum failure strain is exceeded. The minimum failure strain ensures that at least one flaw is active for large enough volumes, the minimum failure strain is given by(Melosh et al.,). Following Benz and Asphaug (), we set the minimum failure strain of each element of volume in our mesh towhereis the volume of material in our simulation andis a random number between 1 andwhereis the volume of the smallest cell in a simulation. Because our simulations are axisymmetric, cell volume increases with distance from the symmetry axis. When the volume of a cell,, exceeds, we calculatefor that celltimes. We then take the minimum value ofcalculated to represent the weakest flaw in the cell's volume. As Benz and Asphaug () describe, the use of equation A4 is for a local minimum failure strain which ensures the weakest flaw for a large enough volume,, is active at, but locally the minimum failure strain does not depend on the size of the target or mesh. Figure A1 shows a plot of the minimum failure strain for a simulation replicating a laboratory impact fragmentation experiment. As expected cells nearer to the symmetry axis tend to have higher minimum failure strains due to the nature of cylindrical symmetry of the mesh (cells closer to the symmetry axis have smaller volumes than those further away). Therefore, it is important to view the trend for minimum failure strain in Figure A1 as representative for all other simulations carried out in this work. Figure A1 Open in figure viewer PowerPoint 1984 A plot of minimum failure strain for the modeling of the (Takagi et al.,) experiments. The color bar on the right represents fragment size, and material is colored accordingly. (A5) 2004 2004 2004 2004 2004 If the minimum failure strain in the cell is exceeded, damage is evolved according to equation A1 in two half‐time steps. In between the two damage half steps is a full step calculation for fracture area:The fracture area and minimum failure strain are tracked as a cell‐centered variable and advected in the same way as material damage (Collins et al.,). Stresses are then updated according to the damage increment in the same way as described by Collins et al. (). One difference from Collins et al. () is that if material is in tension (i.e., when a principal stress is tensile), we assume failure is in tension and do not consider shear failure. This avoids the possibility that the shear failure routine can double count failure leading to accumulation of shear damage even if material is in uniaxial tension. Once material is completely damaged, we no longer evolve damage or fracture area. Put simply, the tensile failure model of Collins et al. () has been changed to now use the Grady‐Kipp model and the shear failure model of Collins et al. () has remained unchanged. (A6) t f is the time at which the damage is greater than 0.99 (~1). The peak of the fragment size distribution curve is then calculated aswhereis the time at which the damage is greater than 0.99 (~1). To test our implementation of the Grady‐Kipp model and determine appropriate k and m material parameters, we performed simulations that recreate an experiment of Takagi et al. (1984). The same approach was used by Melosh et al. (1992), except that Melosh et al. (1992) used a Lagrangian mesh, whereas we are using an Eulerian mesh. We model ~2‐cm‐diameter basalt impactors striking ~6.5‐cm‐diameter spherical basalt targets at ~600 m/s, at resolutions of 5, 10, 20, and 40 cppr. To decide which k value to use, we vary the values from k = 1034 to k = 1036 to k = 1038, with an m value kept constant at 9.5. Jutzi et al. (2009) found that k and m are not independent and for a given target size if ln(k)/m is kept constant, their results were largely unaffected and it is not important to explore the effects of m and k separately. We apply the same EOS, strength model, and damage model as those used in the simulations described in the main text (see Methods, section 2 in the main text; Table A1). Table A1. Material Parameters (Takagi et al., , Modeling) Parameter Value Reference EOS ANEOS basalt 1 Poisson's ratioa 0.25 Melting temperaturea 1360 K 2 Thermal softening coefficienta 0.7 3 Simon A parameterb 4.5 GPa 2 Simon C parameterb 3.0 2 Frictional coefficient (damaged)a 0.6 3 Frictional coefficient (undamaged)a 1.4 3 Cohesive strength (damaged)a 0.01 MPa 2 Cohesive strength (undamaged)a 20.0 MPa 3 Strength limit (damaged)a 2.5 GPa 2 Strength limit (undamaged)a 2.5 GPa 2 k 1034, 1036, 1038 m 9.5 1984 k and resolution affect the final distribution of fragment sizes. The fragment size results of the high‐resolution runs are shown in Figure m = 9.5 and k = 1036 are in good agreement with experimental results, justifying their use in the simulations performed in this work. The fragment sizes plotted in Figure (A7) The results of modeling the Takagi et al. () experiment is shown in Figure A2 . The local definition of minimum failure strain (equation A4 ) allows us to resolve some fragments with damage less than 1 and an intact core. Figure A2 shows how varyingand resolution affect the final distribution of fragment sizes. The fragment size results of the high‐resolution runs are shown in Figure A3 . Figure A3 shows that9.5 and10are in good agreement with experimental results, justifying their use in the simulations performed in this work. The fragment sizes plotted in Figure A3 use the following equation: Figure A2 Open in figure viewer PowerPoint k and resolution (cppr), for the modeling of the Takagi et al. ( 1984 k = 1034, 1036, and 1038 runs occupy the top, middle, and bottom rows, respectively. The cppr steps vary from 5 (first column), to 10 (second column), to 20 (third column), and to 40 (fourth column). cppr = cells per projectile radius. A comparison of fragment size and how it changes with differing values ofand resolution (cppr), for the modeling of the Takagi et al. () experiments. The color bar on the right represents fragment size, and material is colored accordingly. Material is left gray if it is 99% or less damaged. The10, 10, and 10runs occupy the top, middle, and bottom rows, respectively. The cppr steps vary from 5 (first column), to 10 (second column), to 20 (third column), and to 40 (fourth column). cppr = cells per projectile radius. Figure A3 Open in figure viewer PowerPoint M fragment /M target ) of fragments. The blue diamonds represent the results of the Takagi et al. ( 1984 1992 1984 m values kept at 9.5 and k = 1038, 1036, and 1034, respectively. The values used in the teal, gold, and purple lines were gathered from our 40 cells per projectile radius runs. Plot of cumulative number of fragments versus the normalized mass () of fragments. The blue diamonds represent the results of the Takagi et al. () experiment; the orange triangles represent the results of Melosh et al. (). The teal, gold, and purple lines represent the results of our modeling of the Takagi et al. () experiment with allvalues kept at 9.5 and= 10, 10, and 10, respectively. The values used in the teal, gold, and purple lines were gathered from our 40 cells per projectile radius runs. L, in the ith cell, where , and where V cell is the volume of the cell. Thus, the net fragment size distribution is the sum of , given by (A8) n is the number of fractured cells. The size of resolved fragments is not included in the calculation. Producing and testing an algorithm to estimate the sizes of the resolved fragments may be the subject of future work and would likely account for the poorer fit at larger fragment sizes in Figure which is the number of cumulative fragments larger than the fragment diameter,in theth cell, where, andwhereis the volume of the cell. Thus, the net fragment size distribution is the sum of, given bywhereis the number of fractured cells. The size of resolved fragments is not included in the calculation. Producing and testing an algorithm to estimate the sizes of the resolved fragments may be the subject of future work and would likely account for the poorer fit at larger fragment sizes in Figure A3 (the axisymmetric nature of iSALE‐2D means that all resolved fragments are actually annuli and it is not clear what the size of such a fragment should be). The general shape of the fragmented region and calculated fragment sizes are relatively insensitive to model resolution (Figure A2 ). One notable change is that more small fragments are seen at higher resolution. The plotted fragment sizes are actually the peak (equation A6 ) in the fragment size distribution (equation A7 ). The smaller fragments appearing at higher resolution are actually part of the distribution of fragments contained in the larger cells of the lower‐resolution simulations. However, it is clear in Figure A4 that the difference between the 5 and 40 cppr is quite stark, but between the 20‐ and 40‐cppr runs the difference is less pronounced; once well resolved, resolution effects are minimal. Figure A4 Open in figure viewer PowerPoint M fragment /M target ) of fragments. The blue diamonds represent the results of the Takagi et al. ( 1984 1992 1984 m = 9.5, k = 1036, and cpprs of 5, 10, 20, and 40, respectively. cppr = cells per projectile radius. Plot of cumulative number of fragments versus the normalized mass () of fragments. The blue diamonds represent the results of the Takagi et al. () experiment; the orange triangles represent the results of Melosh et al. (). The gold, purple, green, and teal lines represent the results of our modeling of the Takagi et al. () experiment with= 9.5,= 10, and cpprs of 5, 10, 20, and 40, respectively. cppr = cells per projectile radius. 3For the sake of completeness, we felt it prudent to show the results of two additional 1‐km‐diameter impactor runs that highlight how Grady‐Kipp fragmentation affects the distribution of damage as compared to iSALE‐2D Dellen using the damage model of Collins et al. (2004), as well as how resolution can affect the distribution of fragment size on a larger scale than in the simulations of the Takagi et al. (1984) experiments. These two runs use the same parameters that were outlined in section 2 of the main text (see Table 1). In Figure A5, the unaltered iSALE‐2D damage field is plotted on the left, while the damage field of the run using iSALE‐2D with Grady‐Kipp implementation is on the right. The extent of damaged material from the point of contact is greater with Grady‐Kipp, and the material is less uniformly damaged, while the iSALE‐2D Dellen run is almost all either 100% damaged or not at all. Aside from these variances, however, there is not an appreciable difference between the two models. A comparison of 10 and 20 cppr for the 1‐km impactor simulation is given in Figure A6. This plot highlights that the typical fragment size in the 20‐cppr run is smaller than that in the 10‐cppr run. Other than this change in typical fragment size, the two runs are relatively similar, with nearly equivalent shape and damage distribution. The runs simulating lunar impacts are less sensitive to changes in resolution than the simulations of the Takagi et al. (1984) experiments (e.g., Figure A6 compared to Figures A2 and A4). This observation when combined with the long run times associated with higher‐resolution simulations justifies our use of 10 cppr for our simulations of lunar impacts. Figure A5 Open in figure viewer PowerPoint Mirrored plot of damage 10 s after a 1‐km‐diameter impactor strikes a planar target at 15 km/s. On the left is the result using iSALE‐2D with no implementation of the Grady‐Kipp fragmentation model and with a tensile strength of 10 MPa, and on the right is the result using our implementation of Grady‐Kipp fragmentation. The results are similar but with a slightly less damaged hemispherical tensionally fragmented zone that extends further from the point of contact. Figure A6 Open in figure viewer PowerPoint Mirrored plot of fragment size 10 s after a 1‐km‐diameter impactor strikes a planar target at 15 km/s, with a resolution of 10 cppr on the left, and 20 cppr on the right. All colored material is material that was at least 99% damaged. With greater resolution, the typical fragment size within the HTFZ is lower, while the overall shape of the HTFZ is approximately equivalent. cppr = cells per projectile radius; HTFZ = hemispherical tensionally fragmented zone.