The Relative Effects of Surface and Subsurface Morphology on the Deflection Efficiency of Kinetic Impactors: Implications for the DART Mission

The Double Asteroid Redirection Test (DART) mission impacted Dimorphos, the moonlet of the binary asteroid 65803 Didymos, on 2022 September 26 and successfully tested a kinetic impactor as an asteroid deflection technique. The success of the deflection was partly due to the momentum of the excavated ejecta material, which provided an extra push to change Dimorphos’s orbital period. Preimpact images provided constraints on the surface but not the subsurface morphology of Dimorphos. DART observations indicated that Dimorphos contained a boulder-strewn surface, with an impact site located between a cluster of large surface boulders. In order to better understand the momentum enhancement factor (β) resulting from the impact, we performed impact simulations into two types of targets: idealized homogeneous targets with a single boulder of varying size and buried depth at the impact site and an assembly of boulders at the impact site with subsurface layers. We investigated the relative effects of surface morphology to subsurface morphology to put constraints on the modeling phase space for DART following impact. We found that surface features created a 30%–96% armoring effect on β, with large surface boulders measuring on the order of the spacecraft bus creating the largest effect. Subsurface effects were more subtle (3%–23%) and resulted in an antiarmoring effect on β, even when layers/boulders were close to the surface. We also compared our 2D axisymmetric models to a 3D rectilinear model to understand the effects of grid geometry and dimension on deflection efficiency computational results.


Introduction to the DART Kinetic Impactor Test
Understanding how to deflect an asteroid on a collision course with Earth is of great importance and a focus of research undertaken by the planetary defense community.Deflection of an asteroid by a kinetic impactor is one mitigation method that has a high degree of technological maturity (National Research Council 2010;Dearborn & Miller 2015).The Double Asteroid Redirection Test (DART) mission recently successfully demonstrated this mitigation method.On 2022 September 26, at 23:14 UTC, DART made history by executing the first fullscale technology demonstration of asteroid deflection via a kinetic impactor (Cheng et al. 2018;Daly et al. 2023).The DART spacecraft autonomously navigated to and impacted its intended target, the secondary body (Dimorphos) of a nonthreatening binary asteroid system (65803 Didymos).This impact resulted in massive rays of ejecta that contributed to a reduction in the binary orbital period of 33.0 +/−1.0 (3σ) minutes, measured via ground-based telescopes and radar (Thomas et al. 2023), indicating an orbital velocity change of 2.7 mm s −1 (Cheng et al. 2023).
The efficiency of a deflection by a kinetic impactor is characterized by the momentum enhancement factor, β (Ahrens & Harris 1994;Holsapple & Housen 2012;Rivkin et al. 2021;Statler et al. 2022).In its simplest form, β is the ratio of the change of momentum of the target following impact (ΔP T ) to the incoming impactor momentum (P i ): Equation (1) shows the ideal equation of β for a normal impact and is a function of the target's mass (m T ), the change in velocity of the target (Δv), and the momentum of the ejecta (P e ).It is helpful to rewrite β in terms of the ejecta momentum because, for hypervelocity impacts, the change in momentum of the asteroid can be amplified by the momentum of the crater ejecta that exceeds the escape velocity of the target asteroid (Housen & Holsapple 2011).Therefore, β greater than 1 indicates that excavated ejecta contributed in imparting additional momentum to the system (Nysmith & Denardo 1969).Despite efforts in numerical modeling, scaling rules, and laboratory-scale experiments, there exists significant quantitative uncertainty in predicting β for planetary-scale impacts.This is because β varies significantly depending on the target asteroid's material properties and composition (Holsapple & Housen 2012;Jutzi & Michel 2014;Stickle et al. 2015Stickle et al. , 2022;;Bruck Syal et al. 2016;Raducan et al. 2019;DeCoster et al. 2022;Luther et al. 2022).Based on previous work on the effects of surface morphology (e.g., surface boulders at the impact site; Barnouin et al. 2019a;Graninger et al. 2022) and subsurface morphology (e.g., layers; Quaide & Oberbeck 1968;Raducan et al. 2020;Ormö et al. 2022), it is likely that the target structure also plays a major role in determining β.
In the case of a threatening object on a collision course with Earth, one way to learn more about the threatening asteroid is to perform a flyby reconnaissance in advance of a mitigation mission.A flyby spacecraft would likely collect characterization data in the form of images to inform the mitigation mission design.These images would likely capture the surface morphology; however, characterizing the subsurface structure would be much more challenging.This is similar to the data taken by the DART spacecraft prior to impact with Dimorphos.The DRACO camera (Fletcher et al. 2018) took images up until <2 s before impact, revealing a rubble-strewn surface (Figure 1).Even with these images, much about the asteroid remained unknown, which complicates the interpretation of the results and future numerical simulations.
Previous missions to asteroids indicated that the surface structure and likely the substructure can vary significantly from one asteroid to another (Sullivan et al. 2002;Miyamoto et al. 2007; Barnouin et al. 2019a;Rozitis et al. 2020;Jutzi et al. 2022).Generally, asteroids are categorized into three categories: solid objects made of strong, cohesive material; rubble piles containing strong pieces within a weaker matrix held together by self-gravity, small particles, and cohesion; and fractured shards consisting of a relatively strong body with large-scale fractures resulting from previous impact events (Bottke et al. 2002).Although the internal structure of asteroids has not been directly observed, models of their formation suggest that each of these categories has unique internal characteristics, including different material properties that can significantly affect the outcome of an impact event (Scheeres et al. 2015b).For example, the strength and cohesion of asteroid material influence the formation of ejecta after impact, while the presence of fractures can modify the energy absorption and disruption within the asteroid (Holsapple 1993;Asphaug et al. 2002;Holsapple & Housen 2012;Bruck Syal et al. 2016;Luther et al. 2018;Barnouin et al. 2019b;Raducan et al. 2019).Therefore, comprehending the surface and subsurface structure and constituent materials of asteroids is vital for predicting deflection efficiency.
In the case of the DART impact, little was known about Dimorphos prior to impact, which led to preimpact simulation studies that considered unconstrained phase spaces (Bruck Syal et al. 2016;Raducan et al. 2019;Rainey et al. 2020;Stickle et al. 2022).The successful execution of the DART mission provided a wealth of observations (shown in light blue in Figure 2) from DRACO (Fletcher et al. 2018), the Italian Space Agency's Light Italian Cubesat for Imaging of Asteroids (Dotto et al. 2021), and terrestrial and space-based telescopes to constrain the preimpact unknown parameters.Following impact, it was evident that a significant part of the momentum change experienced by Dimorphos was due to momentum carried by ejecta (Cheng et al. 2023).
The momentum carried by the ejecta is a complicated, timedependent, nonlinear function of multiple parameters, including impact angle, surface geology, topography, and target material properties (Housen & Holsapple 2011;Hoerth et al. 2013;Raducan et al. 2019;Rivkin et al. 2021;Jutzi et al. 2022;Statler et al. 2022;Stickle et al. 2022).Therefore, high-fidelity impact simulations are required to understand the relationship between the ejecta momentum and the kinetic impact scenario.Figure 2 shows how general observables obtained from a kinetic impactor planetary defense mission like DART constrain simulation input parameters.Typically, little is known about the target asteroid before impact; however, the impactor properties are well constrained.After impact, close-up images of the asteroid along with standoff observations from Earth and space telescopes can provide information about the asteroid's shape, surface morphology, and material properties (bulk porosity, approximate density, etc.).For example, DRACO images showed that Dimorphos is an oblate spheroid with a boulderstrewn surface resembling a rubble pile (Figure 1; Daly et al. 2023).Additionally, DRACO provided critical information about the impact conditions occurring between two large surface boulders (Daly et al. 2023).Therefore, some material model parameters (target material parameters like strength, internal friction, and surface morphology) can be constrained from observations (Figure 2), and some (e.g., porosity, density, and subsurface morphology) cannot.In Dimorphos's case, future observations from the European Space Agency's (ESA) Hera mission will further constrain forthcoming impact models, including subsurface radar sounding measurements (i.e., lowfrequency radio frequency) that will provide details about the internal structure of Dimorphos.Even more recent asteroid exploration has revealed a layered structure of Ryugu (Arakawa et al. 2020;Kadono et al. 2020) and Bennu (Bierhaus et al. 2023), which further motivates the investigation of substructure effects on kinetic impacts.
The effects of subsurface morphology on impact cratering have been experimentally and computationally investigated for layered targets (Oberbeck & Quaide 1967;Quaide & Oberbeck 1968;Senft & Stewart 2007, 2008;Stickle & Schultz 2012, 2013;Prieur et al. 2018;Raducan et al. 2020).They suggest that layer thickness and properties (e.g., strength, porosity, etc.) affect both the resulting crater morphometry (e.g., Quaide & Oberbeck 1968) and the ejecta dynamics, including the velocities, angles, and volumes of ejecta (Senft & Stewart 2007;Raducan et al. 2020).Weaker, porous layers have been shown to reduce β and crater size, which is indicative of armoring (Raducan et al. 2020).In contrast, the presence of a stronger substrate layer close to the surface has resulted in an amplification (antiarmoring effect) of the amount of ejected mass (and β).However, this amplification is sensitive to the material properties of the layers, where a substrate with small porosity reduced or even negated the antiarmoring effect.Raducan et al. found that the optimum layer thickness that produced the maximum ejected momentum depended on the strength and density ratios of the layers (Raducan et al. 2020).
Boulders, rather than layers, can also result in armoring or antiarmoring behavior, depending on the boulder arrangement within the target.Simulations by Graninger et al. (2022) showed that boulders at the surface reduced β (armoring effect), while the presence of many boulders just below the surface (with no boulders at the surface) resulted in an enhanced β (antiarmoring effect), though this effect was relatively minor (∼10% change).
Further, surface morphology and local boulder field distributions can also affect the ejecta dynamics, ejecta mass, deflection efficiency, and total delivered Δv magnitude (Scheeres et al. 2015a).Ormö et al. (2022) concluded that boulders reduced the momentum enhancement compared to a pure sand target (Ormö et al. 2022).Additionally, Barnouin et al. (2019b) experimentally investigated kinetic impacts into coarsely fragmented material with physical heterogeneities at the surface and found nearly a factor of 2 variability on the impact efficiency due to natural variation in impact conditions due to variation in impact coupling (Barnouin et al. 2019b).
Altogether, it is not clear from previous work what the relative effects of heterogeneities in the surface and subsurface morphology are on the deflection efficiency observables relating to a kinetic impact.DART provides an interesting example case because the surface morphology of the impact site was well constrained from DART mission data, but the subsurface morphology was not.Here, the rubble-strewn surface of the DART target Dimorphos provides inspiration for the setup of a rubble-pile asteroid target to better understand how subsurface morphology would affect deflection by a kinetic impactor.The problem is complex; therefore, an idealized study of the effects of an individual boulder at the impact point provides an anchor to better understand the physical (anti)armoring process before extrapolating to more complex asteroid-like targets and impact by spacecraft.

Methods/Numerical Models
We used the iSALE (Amsden et al. 1980;Wünnemann et al. 2006) shock physics code in two dimensions and CTH (McGlaun et al. 1990;Crawford 1999) in two and three dimensions to model DART-like impacts on asteroid surfaces to quantify the effects of surface and subsurface morphology on deflection efficiency.Simulations using iSALE-2D examined an idealized study of the effects of a single boulder surrounded by weak regolith material.In this study, we varied the size and burial depth of the boulder at the impact site.Next, we performed simulations using iSALE-2D and CTH (2D and 3D) to understand the effects of an assembly of boulders at the impact site in a rubble-pile target.We also modified the subsurface morphology by including core layers at varying depths from the rubble-pile surface.The Supporting Information (S.1) provides full details of the material models and the analysis procedure for the calculated β, the ejecta mass and velocity distributions, and the crater size and is summarized below.

2D iSALE
The iSALE code is based on the SALE hydrocode solution algorithm (Amsden et al. 1980), which was modified to include fragmentation and elastic-plastic constitutive models, multiple materials, and various equations of state (EOSs; Melosh et al. 1992;Ivanov et al. 1997).Later on, a modified strength model (Collins et al. 2004) and a porosity compaction model Figure 2. Schematic detailing the parameter phase space for modeling a kinetic impact for planetary defense in reference to the DART mission.The known parameters are in pink, and the unknowns are in orange.For DART, observations made before and after impact (light blue) provided constraints on the preimpact unknowns, which allowed for further constraints on impact model parameters (green).Future observations from the ESA's Hera mission (dark blue) will further constrain forthcoming impact models.(Wünnemann et al. 2006;Collins et al. 2011) were included.iSALE has been benchmarked against other shock physics codes (Pierazzo et al. 2008;Stickle et al. 2020;Luther et al. 2022).Validation tests against laboratory hypervelocity impact experiments have been performed in various studies to confirm the accuracy of the numerical results in terms of crater formation (Güldemeister et al. 2015;Ormö et al. 2015;Wünnemann et al. 2016;Winkler et al. 2018) and material ejection (Wünnemann et al. 2016;Luther et al. 2018;Raducan et al. 2019).The iSALE code has also been applied to study shock metamorphism in rubble-pile asteroids, for which heterogeneities have been explicitly resolved (Güldemeister et al. 2022).In addition, the code has been used to study the effect of target heterogeneities in various materials including the investigation of mesoscale (i.e., resolving heterogeneities) and macroscale (i.e., bulk material behavior) models conducted with iSALE-2D and 3D that were successfully compared to observations from laboratory experiments and meteorite samples (Güldemeister et al. 2013;Kowitz et al. 2013;Davison et al. 2016;Kowitz et al. 2016;Moreau et al. 2017).

Idealized Study of the Effect of an Individual Boulder at the Impact Site
Our simulations of an idealized impact onto a target comprising a single boulder embedded in a porous weak regolith followed previous simulations of impacts in a regolith halfspace (Raducan et al. 2019).The impactor was modeled as a hollow, thin-shelled aluminum cylinder (Raducan et al. 2022a) with a diameter and height of 80 cm, a shell thickness of 10 cm, and a mass of 628 kg.The vertical impact velocity was 6.5 km s −1 .The thermal behavior of the impactor material was approximated by a Tillotson EOS for aluminum (Tillotson 1962) and the mechanical response by a Johnson-Cook model (Johnson & Cook 1983).
In these simulations, we used a Tillotson EOS for basalt, coupled with the epsilon-alpha porosity compaction model to describe the volumetric response of the boulder and regolith.The solid basalt material had a grain density of 2.650 g cm −3 , and the bulk porosities for the boulder and regolith were 7.4% and 20.0%, respectively.A simple pressure-dependent, damage-independent yield curve with a zero-pressure cohesion of 10 Pa was used to describe the strength of the regolith (Lundborg 1968;Raducan et al. 2019).Boulder strength was modeled using a more complex damage-and pressure-dependent strength model (Collins et al. 2004) with an intact cohesion of 1 MPa (and a damaged cohesion of 10 Pa).Table A1 lists the complete suite of material model parameters used in the idealized study.The intact boulder was substantially stronger, stiffer, and denser than the surrounding regolith.The simulation employed the regridding technique described by Raducan et al. (2019).The simulation started with a cells per projectile radius (cppr) of 40 (the shell was resolved by 10 cells) and then was regridded to coarsen the mesh by a factor of 2 four times during the simulations so that the final resolution for most of the crater growth was cppr = 2.5 (16 cm cell -1 ).

Study of the Effect of Boulder Arrangements and the Subsurface Core
In addition to an idealized study, we also performed a study of the effects of surface boulder field arrangements and varying subsurface morphologies.In these simulations, which sought to approximate the DART experiment, we defined a 579 kg impactor.Our impactor was modeled by three aluminum spheres, with a more massive central sphere of 503 kg that represented the spacecraft bus and two smaller spheres of 38 kg each representing the solar panels (Owen et al. 2022).The centers of the small spheres were 2.3 m distant from the center of the main sphere.The main sphere was resolved by 15 cppr and a cell size of ∼3.3 cm.In order for our projectile to represent the spacecraft bus with a diameter of about 1 m without explicitly resolving any detailed structures, the projectile density was set to 1000 kg m −3 .The impact velocity was set to 6.1446 km s −1 .Note that we chose a different projectile geometry in this section with respect to the idealized study (three spheres in series versus a single hollow cylinder) in order to better represent the complicated DART spacecraft in a larger computational domain.This is because Owen et al. (2022) showed that in both iSALE and CTH, three spheres (in series) resulted in a closer approximation for β to the full DART spacecraft than a solid cylinder (which tended to lead to overpredictions; Owen et al. 2022).Although Owen et al. (2022) did not investigate hollow projectiles, Raducan et al. (2022a) did in iSALE and found that solid cylinders produced similar βs to thick-and thin-shelled cylinders, implying that three spheres in series is the best simplifying representation of the complicated DART spacecraft in a 2D axisymmetric grid (Raducan et al. 2022a;Owen et al. 2022).
The target asteroid consisted of many boulders, which we placed in a relatively weak matrix.We used a basalt Tillotson EOS (Benz & Asphaug 1999;Luther et al. 2022) to define the boulders and a quartzite ANEOS (Melosh 2007) to define the matrix material.The strength of the boulders was simulated using a complex strength model with transition from an yield envelope of the intact state to the one of a damaged state, according to the accumulated damage (Collins et al. 2004).The matrix strength was described by a Drucker-Prager model (Drucker & Prager 1952) with 10 Pa cohesion and a coefficient of friction of 0.6.Table A2 summarizes the material model parameters.The initial conditions of a subset of the iSALE simulations are illustrated in Figure 3. Additionally, the Supporting Information (S.1.2.2) provides further details about the gridding and the schematics of the initial conditions of the simulations (Figure A1), as well as a discussion of how the rubble-pile targets were constructed.

2D CTH Simulations: The Effect of Boulder Arrangements
We performed 2D axisymmetrical (2DC) numerical simulations using the CTH hydrocode to evaluate the effects of surface boulder arrangements and subsurface structure on DART-like impacts.CTH is a multidimensional, multimaterial, large-deformation, strong shock-wave physics code developed by Sandia National Laboratories (McGlaun et al. 1990;Crawford 1999).CTH is a two-step Eulerian code that uses a continuum representation of materials and has an adaptive mesh refinement (AMR) capability, which allows selective refinement of areas of interest to enable increased computational efficiency.We refined the mesh along the shock front and in areas where there was a material interface.The development history and description of the models and novel features of CTH are described elsewhere (McGlaun et al. 1990;Trucano & McGlaun et al. 1990).
CTH has been extensively validated for hypervelocity impacts, as well as benchmarked against both laboratory-scale impact experiments and other numerical shock physics codes (Grady & Winfree 2001;Pierazzo et al. 2008;Caldwell et al. 2019;Stickle et al. 2020).
In general, these previous studies showed that the prediction of momentum enhancement is highly resolution-dependent, requiring higher resolutions (cppr 10) to converge compared to the crater size.We performed our 2D simulations at cppr = 10, calculated with respect to the radius of the central sphere of our projectile.See the Supporting Information, section S.1.3,for details about the grid resolution convergence study.Further, we note that the choices of the strength model and material porosity have larger effects on the resulting crater size and β compared to the variation among codes (Stickle et al. 2020).Therefore, we used similar material EOSs and strength models between our simulations in CTH and iSALE to maintain consistency.We performed 2DC simulations of rubble-pile targets using CTH and added a core layer at varying depths beneath the rubble-pile surface.The rubble piles consisted of competent boulders placed in a granular-like continuum matrix material.The EOS for basalt was used for the boulders, which is commonly used for asteroid impact simulations (Stickle et al. 2022).The 2D rubble-pile targets were constructed using ParticlePack, a program created by Lawrence Livermore National Laboratory that generated geometry models of particles packed into an enclosure for use in hydrocode simulations (Friedman 2012).Rubble-pile realizations were created by randomly placing circular boulders (ranging in diameter from 1.0 to 10.0 m) within a circular asteroid (163 m diameter).The boulder sizes followed a power-law distribution derived from scaling the rock size distribution on Itokawa (Tancredi et al. 2015).Note that this boulder size distribution was different from what was reported for Dimorphos (Weibull distribution; Pajola 2023).The boulder fill fraction was ∼62%, and the remaining asteroid volume was filled with the matrix material (60% porosity), resulting in a total asteroid density of 2.18 g cm −3 .We note that this density was lower than the reported nominal Dimorphos density (2.40 g cm −3 ); however, it was still within the 1σ confidence interval reported for the Didymos system (Davison et al. 2016).We chose to run our models at this lower density because a parallel study showed that a highly porous matrix material leads to more boulder ejection (DeCoster et al. 2023), which was a major observable following the DART impact (Jewitt et al. 2023).The choice of material properties was motivated by DART mission observations, and they are listed in Table A3.
We used multiple clustered spherical projectiles to represent the DART spacecraft, which is a better approximation to the more complex spacecraft model compared to a single spherical projectile while balancing computational efficiency (Owen et al. 2022).Note that Owen et al. (2022) recommended representing the complicated DART spacecraft as three spheres in parallel; however, we represented the projectile as three spheres arranged in series due to the axisymmetric nature of our 2D grids (three spheres in parallel would result in a toroid-shaped projectile).We note that DART did not hit solar panel, then bus, then solar panel; however, Owen et al. showed that representing the complex DART projectile as three spheres in series more closely captured the deflection efficiency response of the complex DART spacecraft than a single spherical projectile, so we adopted this projectile form for our study.The central sphere had r = 35.35cm (88% of the total mass), and the two auxiliary spheres each had r = 15.30cm (each contained 6% of the total mass), which was the approximate ratio of the main DART spacecraft bus to the Roll Out Solar Arrays (ROSAs).The auxiliary spheres were placed 2.33 m away from the central sphere (measured center to center), which equated to the width of a single ROSA.We individually varied the surface morphology (i.e., boulder size at the impact point), the subsurface morphology (i.e., addition of a strong/competent or weak granular basalt core at different thicknesses below the rubble-pile surface), and the cohesive and tensile strengths of the matrix material to understand their relative effects on β, ejecta mass-velocity distribution, and crater size.We chose basalt as the representative EOS for the asteroid boulders and core because it is generally well understood and commonly used in other impact simulation studies (Raducan et al. 2019;Stickle et al. 2020Stickle et al. , 2022)).The boulder sizes at the surface ranged from diameters of 1.0 m (small surface boulders) to 6.0 m (large surface boulders).Boulders smaller than 1.0 m in diameter were too computationally expensive to simulate, and discrete boulders below 1.0 m were assumed to make up the matrix material, which was represented as a continuum here.The larger surface boulders were consistent with the size of the surface boulders observed at the DART impact site on Dimorphos (Davison et al. 2016).The subsurface core layer ranged from 0.88 to 15.00 m below the rubble-pile surface.We report the core placement as a function of the rubble-pile top-layer thicknesses (h) over the central sphere projectile radius (a) so that 2.47 h/a 42.25, allowing for the investigation of a wide range of subsurface effects.Figure 4 illustrates examples of the different target morphologies that were simulated in CTH, which map to 13 individual simulations detailed in Table 1.Table A3 lists the parameters used in the CTH material models.

3D CTH Simulations
The 2D simulations are useful for examining the influence of parameters and allow for larger parameter-space sweeps and examination of specific trends as parameters change (Caldwell et al. 2020, p. 16;Pierazzo et al. 2008;Heberling et al. 2017;Raducan et al. 2020Raducan et al. , 2022a;;Rainey et al. 2020;DeCoster et al. 2022;Luther et al. 2022).However, a 2DC representation does not faithfully replicate all stress states, especially for complex targets, and can lead to artifacts such as accentuation of damage along the centerline (Trucano & McGlaun et al. 1990).Thus, it is important to use 3D simulations to evaluate effects when possible.To this end, we performed a 3D rectilinear (3DR) simulation of a DART-like impact into a rubble-pile target constructed with a similar geometry to the 2D rubble piles and  Note.Large boulders at the surface had radius = 3 m, and small boulders at the surface had radius = 0.5 m.See Table A3 for material model parameters.
using the same boulder-packing algorithm in order to compare and contrast the 2DC and 3DR results.The 2D and 3D rubble-pile targets were both constructed using ParticlePack (Friedman 2012).For the 3D version of the DART-like impact model, spherical "boulders" were randomly packed in a spherical "asteroid" with a diameter of 163 m.Dimorphos's true shape is an oblate spheroid, where the extents along the principal x-, y-, and z-axes are 177, 174, and 116 m, respectively (Davison et al. 2016).We used a spherical diameter of 163 m, which overpredicts the total volume derived for Dimorphos by ∼26% but more closely approximates the extent along the principal x-and y-axes, which is important for precisely capturing the effects of boulder packing in the interior of the rubble-pile structure near the impact location (Raducan et al. 2022b).The 3D asteroid was packed to a fill fraction of 59% with 64,363 boulders that ranged in diameter from 2.0 to 10 m, with a distribution of sizes following the power-law distribution derived from Tancredi et al. (2015).We increased the minimum boulder diameter size compared to the 2D simulation due to computational resource limitations.Figure 5 compares the 3D rubble-pile target to the 2D equivalent.Note that the boulder spatial distribution was variable from 2D to 3D even though the boulder fill fraction and distribution were held roughly constant.

Idealized Study: The Effects of an Individual Boulder at the Impact Site
We used iSALE-2D to perform an idealized study of a hollow, thin-shelled cylindrical projectile directly impacting a single boulder, where both the boulder radius (r) and the elevation of the top of the boulder relative to the surface (z) were varied.Figure 6 plots β−1 (at 2 s postimpact) as a function of the normalized boulder position for boulders of varying size (the radius ranged from 0.2 to 1.5 m).The case where no boulder was present and the projectile impacted only sandy regolith represented the boundary between two regimes in momentum enhancement factor effects: antiarmoring and  armoring.Antiarmoring occurred when the target morphology caused an increase in the deflection efficiency contributed by the ejecta (β−1) from the nominal no-boulder case, and armoring resulted in a reduction.Figure 6 shows that the placement of a boulder at or below the surface had an antiarmoring effect; the closer the buried boulder was to the surface, the larger the effect.As the boulder protruded above the surface, an armoring effect occurred that scaled with the protrusion of the exposed boulder.The magnitude of these affects was dependent on the boulder size, where larger boulders (r = 1.0-1.5 m) had a greater effect on (anti)armoring than smaller boulders (r = 0.2-0.4m).
Figures 7(a)-(c) compare the provenance of the ejecta as a function of the vertical ejection velocity for three scenarios: a half-buried boulder, a fully buried boulder, and no boulder at the impact site.Compared with the no-boulder case, the fully buried boulder increased the volume of ejecta at all speeds, resulting in an antiarmoring effect that enhanced momentum transfer (Figure 7(d)).Comparison of Figures 7 (b) and (c) shows that the ejection speed has been enhanced throughout the excavated zone, so there is a greater volume of ejecta launched at a given speed (as can also be seen in Figure 7(d)).Comparing where the ejection speed contours intersect the target surface shows that the contours of ejection speed each extend slightly further from the impact point and encompass a slightly greater volume of ejecta.The high-impedance boulder enhanced the coupling of impactor kinetic energy with the target.In the half-buried boulder scenario, however, we observed two competing effects.The volume of low-speed ejecta was dramatically reduced by the dissipation and deflection of impact energy in the boulder protruding above the surface.On the other hand, the higher-impedance boulder material at the contact zone enhanced the speed of the fastestmoving ejecta compared to the case where only regolith material was present.Thus, in this case, we observed a greater mass (and momentum) of high-speed ejecta than in the noboulder scenario but a much-reduced mass (and momentum) of low-speed ejecta (Figure 7(d)).
This idealized study of a projectile directly impacting a single boulder of varying size and burial depth indicated that surface effects (large boulders at the surface) produced an armoring effect in β−1 and that subsurface effects (buried boulders) produced an antiarmoring effect.The magnitude of each effect was dependent on the boulder size.In general, surface effects (i.e., a large boulder protruding above the surface) were more pronounced than subsurface effects (buried boulders).Specifically for the largest boulder (r = 1.5 m), in the buried regime studied, we saw that the highest antiarmoring effect (23% increase compared to the no-boulder case) occurred when it was buried right at the surface (z/r = 0).However, the largest armoring effect (50% decrease compared to the noboulder case) occurred when the boulder was only half-submerged on the surface (z/r = 1).This indicated that surface morphology affected β more than subsurface morphology.

The Effects of Boulder Arrangements at the Impact Site
DART mission observations showed that Dimorphos exhibited a rubble-pile surface, where the DART spacecraft impacted between an assembly of three large surface boulders.The previous section investigated an idealized case of a projectile directly impacting a single boulder of varying size and buried depth in a regolith matrix.This study demonstrated the armoring and antiarmoring effects that boulder size and placement can have on β−1.However, because DART impacted between multiple large surface boulders, we moved beyond the idealized study and investigated the effects of an arrangement of boulders at the impact site on deflection efficiency.In this section, we investigate the relative effects of surface and subsurface morphology on deflection efficiency by changing the size of surface boulders and including a subsurface core layer at varying depths below the surface.In line with the results from our idealized study, we found that the effects of surface morphology with respect to an assembly of boulders outcompeted the effects of subsurface morphology (the presence of a core or buried boulders) when material properties (i.e., strength parameters, porosity, boulder fill fraction) were held constant.Further, large surface boulders resulted in an armoring effect, and subsurface layers resulted in an antiarmoring effect on β−1.

The Effects of Surface Morphology
We studied the surface morphology effects on the deflection efficiency by varying both the size of the surface boulders (d = 1.0 and 6.0 m) at the impact site and the strength (cohesion) of the regolith matrix material in the rubble-pile targets (Y o = 1 and 10 Pa). Figure 8 shows the effects of surface morphology on deflection efficiency.We found that the targets with the weakest regolith (homogeneous regolith and rubble-pile targets) resulted in the largest βs compared to both homogeneous and rubble-pile targets with stronger regolith.This trend shows the important influence of material strength properties, in which a weaker regolith resulted in more ejected material and higher β−1 (Figure 8(a)).We saw that increasing the boulder size at the impact site resulted in a drastic armoring effect in β−1.The magnitude of this effect depended on the strength of the regolith.Large surface boulders reduced β−1 by ∼30% when the regolith had a stronger cohesion (Y o = 10 Pa) compared to ∼96% when it was weaker (Y o = 1 Pa).We found that a weaker regolith material enhanced the armoring effect of large surface boulders on β−1. Figure 8(b) indicates that the case comprised of larger surface boulders and stronger regolith resulted in less total ejecta mass compared to the case with weaker regolith; however, the stronger regolith resulted in a larger population of faster-moving ejecta ensuing more momentum enhancement.We plot the cumulative normalized ejecta mass as a function of the normalized ejecta velocity in Figure 8(b) to understand the effects of surface boulder size on the ejecta dynamics.The plot shows that large surface boulders significantly reduced both the total ejecta mass and the ejecta velocity.In particular, the fastestmoving ejecta from 224 to 597 m s −1 (ejecta velocity normalized by impact velocity = 0.036-0.097)present for the small surface boulder targets were eliminated when large surface boulders were present.Consistent with previous findings, we saw that the weaker targets resulted in relatively more ejecta mass than their stronger counterparts (Bruck Syal et al. 2016;Raducan et al. 2019).However, more ejecta mass did not necessarily map to larger β−1s.For example, for the rubble-pile targets with large surface boulders, the weaker rubble pile resulted in more than ∼2× the total ejecta mass of the stronger rubble pile.However, its corresponding β−1 was 83% smaller compared to the stronger rubble pile with large surface boulders.The difference in β−1 stemmed from the increased population of fast-moving ejecta from ∼42 to 135 m s −1 (normalized ejecta velocity 0.006-0.022)for the rubble pile with stronger regolith.Further, the rubble-pile targets with small surface boulders did not display this regolith strength-dependent trend in ejecta mass and velocity.Figure 9 shows that the total ejecta mass-velocity profiles for rubble-pile targets with small surface boulders and varying regolith strength were similar.Additionally, the individual contributions from the boulder and regolith material to the ejecta were similar, with the boulder material dominating the total ejecta mass.This indicated that the more homogeneous surface morphology provided by the smaller surface boulders (the strength of the boulders was constant in all cases) washed out the effects of the regolith material strength on the ejecta velocity.This trend is in stark contrast to the cumulative ejecta mass-velocity distribution when large surface boulders were present.
Figure 10 plots the resulting depth and width of the transient craters at 130 ms postimpact for the homogeneous and rubblepile targets.Figure B2 and Table B2 provide illustrations and measurements of the resulting transient craters.We found that the homogeneous targets resulted in wider craters than the rubble piles, in which the weaker sandy targets resulted in the largest craters.Crater size did not directly map to β−1, however.For example, the strong homogeneous basalt target resulted in the largest β−1 but one of the smaller crater sizes.Even further, the rubble pile with large surface boulders and weak regolith had by far the smallest β−1 of the targets and resulted in one of the deepest craters.Further, the transient craters resulting from all of the rubble-pile targets clustered together.This implied that surface morphology had a small effect on the crater size.We saw that the large surface boulders reduced the width and increased the depth of the transient crater (compared to the other rubble piles) in the case with weaker (Y o = 1 Pa) regolith.When the regolith was stronger, this effect became marginal.

The Effects of Subsurface Layers
We investigated the effects of subsurface morphology by placing homogeneous cores into the rubble-pile target and varied its extent below the rubble-pile surface from 0.88 to 15.00 m.The core layers took the form of either strong competent basalt (Y o = 1 MPa), which we refer to as a "strong core," or loosely consolidated dry sand (Y o = 1 Pa), which we refer to as a "weak core." Figure 11 plots the temporal evolution of β−1 for simulations performed in CTH and iSALE with variable subsurface morphology.The rubble-pile targets (with and without cores) all exhibited β−1 higher than the homogeneous targets investigated (which were representative of the regolith composition in the rubble piles).This indicated that the addition of subsurface boulders within the rubble pile in addition to core layers imposed an antiarmoring effect (increase in β−1).Next, Figures 11 and 12 show that the addition of a core at varying depths under the rubble-pile surface slightly enhanced (antiarmoring; ∼3%-15%) β−1 compared to the purely rubble-pile structure.Generally, the closer the core layer came to the surface, the greater the enhancement.This trend was more apparent in CTH than in iSALE, which could be a result of differences in simulation times.Generally, these results align well with Raducan et al. (2020).In both studies, thick regolith layers (h/a > 14) had little effect on β; however, thinner regolith layers (h/a < 14) exhibited more complex behavior, as layering influenced both the velocity and the mass of the ejecta.Both studies showed that β was sensitive to the material properties of the layers (i.e., the regolith and the core/substrate).Raducan et al. (2020) modeled rubble-pile targets where the porosity of the regolith and substrate varied from 35% to 50% and 0% to 10%, respectively.In our study, the regolith and core (substrate) porosities were held constant at 60% and 0%, respectively.Raducan et al. found that when the regolith had a porosity of 50% and the substrate had no porosity (similar to our target), layers with h/a < 2 resulted in a large amplification (∼60%) in β−1 compared to impacts into a 50% porous homogeneous half-space.We measured a smaller amplification (∼10%-20%) in β−1 in the regime where h/a < 3.5 (core layer buried at 0.88 and 1.23 m below the surface) compared to the homogeneous weak sand target (Y o = 1 Pa).The differences in amplification likely come from the differences in the regolith porosity (50% versus 60%), where Raducan et al. showed that porosity (which includes the specific crush curve parameters used in the model) can greatly reduce or even negate amplification.
Figure 13 plots the cumulative ejecta mass as a function of the ejecta velocity at 130 ms and 4 s postimpact for the simulations in CTH and iSALE, respectively, for simulations where the subsurface was varied.In both CTH and iSALE, the ejecta mass-velocity profiles for all cases with a subsurface core were very similar, where the effects of both the core size and the core strength (strong basalt versus weak sand) made little difference in the ejecta mass and velocity distributions.We saw previously that the addition of a core layer had very little (∼3%) effect on β; however, we did observe that cores close to the surface resulted in an antiarmoring (enhancement) effect in β up to 15%. Figure 13 shows that in CTH, the rubble-pile with no core resulted in the highest cumulative ejecta mass and a larger population of the slower-moving ejecta from 10 to 172 m s −1 (normalized ejecta velocity = 0.002-0.028)compared to the rubble-piles with cores.Further, the presence of a core resulted in an increased population of the fastest-moving ejecta from 300 m s −1 (normalized ejecta velocity = 0.049) to 386-546 m s −1 (normalized ejecta velocity = 0.063-0.089).This result was similar to what was reported by Raducan et al. (2020), in which a competent core close to the surface created an amplification in β due to increased ejecta velocity (Raducan et al. 2020).In total, we saw that the presence of a heterogeneous core layer close to the surface resulted in an antiarmoring effect on β by increasing the velocity of the excavated ejecta.We saw this trend for both the strong competent basalt core and the weak granular sandy core.
Figure 14 plots the transient crater depth and width (diameter) for both CTH and iSALE, with variable subsurface morphology.Figures B1 and B3 illustrate the transient craters and the material damage for iSALE (3 s postimpact) and CTH (130 ms postimpact).We saw that the damage profiles for the sandy homogeneous targets were very similar despite the weaker homogeneous target producing a larger crater.Damage occurred in boulders within a radius of ∼10 and ∼20 m from the impact point in CTH and iSALE, respectively.The difference in the extent of the damage was likely due to the different length of time the simulations were run.The addition of a layer subtly changed the damage profile in the target.Specifically, the addition of subsurface layers decreased the extent of the damage depth, where the closer the layer was to the surface, the more localized the damage was to the surface (Figures B1 and    crater depth (compared to the homogeneous targets), and the addition of a near-surface competent core reduced the crater depth more than affecting the crater width.We also saw that a near-surface competent core layer changed the morphology of the crater from a bowl shape to a flat-floor crater (Figure B3).In line with the CTH results, craters from iSALE showed a decrease in depth for a shallow subsurface layer at 1 m depth and a similar crater shape (flat-bottom crater).In this case, the crater width significantly increased, while the width remained relatively constant for all other layer depths modeled.The case with a 5 m deep layer showed a slightly increased crater depth, where it extended down to the top of the subsurface layer.In total, we found that subsurface morphology played a much larger role than surface morphology in influencing the transient crater morphometry.

Comparison of 2D and 3D Simulations
We performed 2DC and 3DR simulations of kinetic impacts into similar rubble-pile targets to understand the effects of grid dimensionality on our kinetic impactor study of rubble-pile targets (Figure 5).In the 2DC grids, we were simulating a halfspace with a boundary condition to reflect the simulation across the centerline.This meant that in 2DC, the boulders were not actually spheres, but toroids, which created hoop stresses (normal stress in the tangential direction) that were not present in the 3DR simulations, where the boulders were spherical.Therefore, boulder ejection and the stress/damage profiles in the boulders were different in 2D versus 3D.Further, the symmetry boundary condition (at the centerline) caused accentuated stress at the axis.In all, these differences in stress states can change based on boundary conditions and projections, which can cause the   B1 and B3 for crater values.Note that the vastly different timescales between CTH (130 ms) and iSALE (4 s) means the CTH simulations had not evolved to a "final" crater; therefore, comparisons between the data should be avoided.material to behave differently.Therefore, to understand these differences, we compared the 2DC rubble-pile (with small surface boulders) to a 3DR rubble-pile (with small surface boulders).Figure 15 plots the temporal β−1 and the cumulative normalized ejecta mass versus the normalized ejecta velocity at 130 ms postimpact for the 2DC and 3DR rubble-pile targets (see Figure 5 for a schematic of the initial conditions).The boulder size and spatial distribution were not exactly the same between the 2D and 3D cases, so we did expect some variability in β and the ejecta mass-velocity distributions due to surface morphology effects; however, the initial conditions were similar enough to perform a comparison to understand the effects of grid dimensionality.Figure 15 shows that the 2DC case overpredicted β−1 by ∼36% compared to the 3DR case; however, it appears that the 3D case may approximate the 2D case better at longer times.The ejecta mass-velocity distribution indicated that the 2D and 3D simulations resulted in approximately the same amount of excavated ejecta material (∼1.6 × 10 8 g); however, the ejecta mass-velocity distribution profiles were quite different.Most notably, the 2DC case resulted in a much larger population of faster-moving ejecta compared to the 3DR.Extra stresses in the 2D grid that do not appear in 3D, like hoop stress within the boulders, would be expected to consume energy and promote lower ejection velocities.Instead, we see that the ejecta velocity (ranging from a normalized velocity of 0.002 to 0.05) in 2D was larger.This could be due to the buried boulders at the impact site that are present in 2D and not in 3D (see Figure 5).We found in the idealized study that a boulder present near the surface results in enhanced ejection speeds throughout the excavated zone, leading to a greater volume of ejecta launched at a given speed (see Figure 7(d)).It may also be possible that hoop stresses in undamaged boulders hindered the horizontal movement of material so that more material moved upward (vertically), which would affect the ejecta angle and momentum.
Figure 16 shows the damage profiles in 2D and 3D at 130 ms postimpact.The extent of the damage was greater in the 2DC case compared to 3DR, partially due to the higher density of boulders at the impact site and partially due to the extra hoop stress present in the grid and centerline effects.In total, the 2DC simulations resulted in an overprediction of ∼36% in β−1 (compared to 3D).In both cases, the extent of the damage was larger in the 2DC case compared to the 3DR case.These differences demonstrated the importance of 3D runs for replicating experiments (and the DART impact).However, we found that a 2DC representation of a rubble pile was still valid for understanding the relative trends between surface and subsurface morphology presented in this work.

Summary
The DART mission successfully demonstrated the effectiveness of a kinetic impactor for asteroid deflection.A major component of the success of the technique came from the contribution from the ejecta momentum that was released after impact, which contributed significantly to the Δv and orbital period change of Dimorphos following impact.The ejecta momentum is a complicated function of many variables, including but not limited to the target material properties and the surface and subsurface morphology.Following the encounter, a wealth of observational data gathered by DRACO before DART's impact was used to describe the surface morphology; however, the subsurface of Dimorphos remained unconstrained.Therefore, we computationally investigated the relative role of surface and subsurface morphology on the deflection efficiency of Dimorphos-like targets.DRACO observations revealed that the DART spacecraft impacted between multiple large surface boulders.This motivated two separate studies to understand the surface and subsurface effects.The first was an idealized study, where we directly impacted a single boulder of varying size and buried depth to understand the effect of a boulder at the impact site.The second study investigated the surface and subsurface morphology effects of an impact into an assembly of boulders.In the second study, we investigated the relative effects of large (boulder radius was on the order of the spacecraft bus radius) and small (boulder radius was smaller than the spacecraft bus radius) surface boulders, in addition to a core layer of differing strengths and compositions at varying depths beneath the rubble-pile surface.For both studies, we found that large surface boulders led to an armoring effect, while both buried boulders and the presence of a subsurface core led to an antiarmoring effect on β.The surface morphology had a larger effect (30%-96%) on β than the subsurface morphology (3%-23%); however, the subsurface morphology had a larger effect on the transient crater morphometry.The magnitude of the surface morphology effect on β demonstrated a dependence on the strength of the regolith material.For the stronger-regolith cases, large surface boulders reduced β by ∼30%.When the regolith was an order of magnitude weaker, large surface boulders reduced β by ∼96%.We found that large surface boulders held together by weaker regolith significantly reduced both the total ejecta mass and the ejecta velocity compared to large surface boulders held together by stronger regolith or small surface boulders.In order to assess the accuracy of our 2D study, we compared our results to a 3D rubble-pile simulation.We found that 2DC simulations resulted in a ∼36% overprediction of β−1 compared to the 3D simulation, which we ascribe to additional stresses present in the grid and variations in the subsurface structure that affected the ejecta velocity distribution but not the total ejecta mass.We concluded that 2DC simulations were useful for understanding the general trends between target morphologies; however, we acknowledged that 3D simulations were necessary for accurately representing the full DART impact and replicating its observables.Overall, there is a need to verify the validity of 3D simulations of geologic rubble-pile targets against experiments.Currently, there is a paucity of 3D code validation studies for hypervelocity impacts into rubble-pile targets (Jutzi et al. 2022;Ormö et al. 2022).However, we note that a 3D validation study is ongoing for iSALE (Ormö et al. 2023).Future work may investigate code validation against experimental data, like recent laboratory impact experiments conducted on targets consisting of various sizes of boulders (Barnouin et al. 2019;Kadono et al. 2020;Okawa et al. 2022;Yasui et al. 2022).

A.1.2. CTH
We performed a grid resolution convergence study in order to balance the coarseness of the grid resolution with our limited computational resources.We simulated the full rubble-pile case with small surface boulders (and no subsurface core layer), where the boulder cohesion was set to 1 MPa and the regolith cohesion was set to 1 Pa.We ran this simulation at four different AMR resolutions: cppr = 10, 15, 20, and 30.The results from our convergence study are plotted in Figure A2.We see that higher resolutions result in lower β−1 values, where a grid resolution of cppr = 10 overpredicted β−1 by 26% compared to cppr = 30.Figure A3 plots the normalized cumulative ejecta mass as a function of the normalized ejecta velocity.We see that the coarser grid resolution (cppr = 10) resulted in an overprediction of the mass of ejecta material, particularly the slower-moving ejecta from 6.14 to 59.3 m s −1 (normalized ejecta velocity = 0.001-0.009).The lower-resolution runs also produced more fastermoving ejecta in the grid, but these ejecta accounted for a small amount of the total ejected mass in the grid.
For the 2DC runs, we used a symmetry boundary condition on the left (inner edge), an outflow/vacuum boundary condition on the right (outer) and top edges, and a semi-infinite boundary condition on the bottom edge of the mass.In 3D, we employed an outflow boundary condition on all edges except the bottom, where a semi-infinite medium sound speed-based absorbing boundary condition was used.The selection of the    A3 for the projectile and the target constituents (boulders and regolith).All simulations were performed using CTH's AMR at cppr = 10 (3.5 cm cell -1 ).We used the brittle damage with localized thermal softening (BDL) to model the strong basalt boulders and the geological yield    surface (GEO) model to represent the granular regolith matrix material.
A.1.3.Analysis Procedure For the CTH simulations, we used the analysis procedure described in detail by DeCoster et al. (2022) to calculate β, the ejecta mass-velocity distributions, and the transient crater size (DeCoster et al. 2022, p. 20).For iSALE, we determined the ejecta characteristics based on the approach described by Luther et al. (2018).In brief, β was calculated by defining an "ejecta plane" above the target surface to track ejected material.As material passed through this plane, it was marked as "ejecta."To avoid erroneously marking material that might have remained attached to the crater rim as ejecta, this plane was defined above the surface.The ejecta plane was set to 5.0 and 2.6 m above the preimpact target surface for CTH and iSALE, respectively.The ejecta plane was code-specific due to subtle differences in the response of the target material as it formed the crater rim postimpact.β was then calculated by summing up the total ejecta momentum in the computational domain.To reduce advection errors and optimize computational efficiency, the computational mesh was necessarily limited in size.Here, ejected material that left the computational grid was accounted for at each time step using virtual witness plates placed along the edges of the mesh.The witness plates tracked the mass and velocity of each ejecta particle as it left the mesh, so that the total ejecta momentum was where m i and v i are the respective masses and velocities of the ith piece of ejecta material, and m i,e and v i,e were the masses and velocities of the ith piece of ejecta that have left the grid, respectively.The velocity used was the velocity component perpendicular to the impact plane (y-velocity for 2D and zvelocity for 3D cases).In iSALE, we used Lagrangian tracer particles, which represented the mass of the cell, their initial location, and their location throughout the simulation time.We stored the mass and ejection characteristics for each tracer that fulfilled the ejection criterion (Luther et al. 2018).Finally, β was calculated following where P impactor was the initial impactor momentum.
The transient crater was calculated by measuring the depth and width of the crater at the last time step of the simulation.The transient crater occurred during the contact/compression and excavation stages of impact crater formation, during which the energy from the impact forced the target material down and compressed it, and then material was ejected from the rapidly expanding crater.It is "transient" because it formed during the early stages of impact cratering and is expected to change at later stages.The depth was measured from the impact plane, and the diameter was measured as the distance in the xy-plane for 2D and both the xz-and yz-planes for 3D simulation between cells along the crater boundary profile that intersected the impact plane.

Figure 1 .
Figure 1.DART DRACO image of Dimorphos 11 s before impact; the yellow dot shows the approximate impact site, and the arrow indicates the direction of the Dimorphos +Z (north) axis.Credit: NASA/Johns Hopkins APL.

Figure 3 .
Figure 3. Initial conditions showing the different target surface and subsurface morphologies considered in the iSALE-2D simulations for (a) the full rubble-pile target, (b) the rubble-pile target with a substrate core layer 5 m below the surface, and (c) the rubble-pile target with a substrate core layer 1 m below the surface.The left side of each image plots the damage, and the right side plots the material density.

Figure 4 .
Figure 4. Initial conditions showing the different target surface and subsurface morphologies considered in the 2D CTH simulations.(a) The full rubble-pile target, (b) zoom-in of the surface from the full rubble pile with large surface boulders (diameter r = 6.0 m), (c) zoom-in of the rubble-pile surface with small surface boulders (diameter r = 1.0 m), (d) small surface boulders with a substrate core ranging from 0.88 to 15 m below the surface, and (e) a homogeneous target.

Figure 5 .
Figure 5.Comparison of the 3D spherical rubble-pile and 2D rubble-pile targets.The boulder fill fraction for the 2D and 3D targets was 62% and 59%, respectively.Panel (a) shows the full 3D spherical rubble-pile; (b) shows the 3D xz cut-through with the impact point at x = 0, z = 81.5 m; (c) shows the 2D xy cut-through with the impact point at x = 0, y = 81.5 m; and (d)-(f) show the impact site and projectile geometry for the (d) 3D yz-plane, (e) 3D xz-plane, and (f) 2D simulations.

Figure 6 .
Figure 6.Relative total momentum of ejecta launched in the first 2 s postimpact as a function of boulder size and normalized position (elevation of the top of the boulder relative to the surface, z, divided by boulder radius, r).The antiarmoring regime is above the black horizontal line indicative of β−1 with no boulder present, and the armoring regime falls below this line.

Figure 7 .
Figure 7.Comparison of provenance of ejecta as a function of vertical ejection velocity (V z ) normalized by the impact velocity (U) for (a) a half-buried boulder, (b) a fully submerged boulder, and (c) no boulder.(d) Cumulative ejected momentum as a function of normalized ejection velocity for a boulder with a 3 m diameter at various burial depths.The blue lines indicate where the boulders extended above the surface, and the green lines indicate where the boulders were buried below the surface.

Figure 8 .
Figure 8. CTH results showing (a) the temporal profile of β−1 and (b) the cumulative ejecta mass (M ejecta ) normalized by the impactor mass (m) vs. the ejecta velocity (V ejecta ) normalized by the impactor velocity (u) at 130 ms postimpact.The plots show the results for homogeneous targets and four rubble-pile targets where the boulder cohesion (Y o,b ) was held constant, the cohesion of the regolith material (Y o,reg ) was varied from 1 to 10 Pa, and the boulder size at the surface was varied from d = 1.0 to 6.0 m.

Figure 9 .
Figure 9.The cumulative ejecta mass (M ejecta ) normalized by the impactor mass (m) vs. the ejecta velocity (V ejecta ) normalized by the impactor velocity (u) separated by the boulder (dashed) and regolith (dotted) contributions to the ejecta at 130 ms.

Figure
Figure 10.Plot of the transient crater at 130 ms postimpact for the CTH simulations.Table B2 contains tabulated values, and Figure B2 shows illustrations of the transient craters.
Figure 10.Plot of the transient crater at 130 ms postimpact for the CTH simulations.Table B2 contains tabulated values, and Figure B2 shows illustrations of the transient craters.

Figure 11 .
Figure 11.Temporal evolution of β−1 for simulations with varying subsurface morphology performed in CTH (top) and iSALE (bottom).The shallowest iSALE simulation (blue line) had a core layer placed 1 m below the target surface rather than 0.88 m below the surface (which was the shallowest layer modeled in CTH).Note the vastly different timescales between CTH (150 ms) and iSALE (4 s).The difference in timescale means the CTH simulations had not evolved to a "final" crater.

Figure 12 .
Figure 12. β−1 vs. the normalized regolith layer thickness, where h was the thickness of the regolith layer on top of the core, and a was the radius of the center sphere of the projectile (35.3 cm).CTH results were plotted at 130 ms (filled symbols), and iSALE results were plotted as open symbols 4 s postimpact.Note that the vastly different timescales between CTH (130 ms) and iSALE (4 s) means the CTH simulations had not evolved to a "final" crater; therefore, comparisons between the data should be avoided.

Figure 13 .
Figure 13.The cumulative ejecta mass (M ejecta ) normalized by the impactor mass (m) vs. the ejecta velocity (V ejecta ) normalized by the impactor velocity (u) at 130 ms postimpact for CTH (left) and 4 s postimpact for iSALE (right).

Figure 14 .
Figure 14.Crater width and depth metrics.CTH data are plotted as filled symbols at 130 ms postimpact, and iSALE data are plotted as open symbols at 4 s postimpact.See TablesB1 and B3for crater values.Note that the vastly different timescales between CTH (130 ms) and iSALE (4 s) means the CTH simulations had not evolved to a "final" crater; therefore, comparisons between the data should be avoided.

Figure 15 .
Figure 15.Comparison of (top) β−1 and (bottom) the cumulative normalized ejecta mass vs. the normalized ejecta velocity for the 2DC rubble-pile with small surface boulders to the 3DR rubble-pile where the cohesion strengths of the boulders and regolith were set to Y o,b = 1 MPa and Y o,reg = 1 Pa, respectively.

Figure A1 .
Figure A1.Initial conditions showing the right half of the rubble-pile setup in iSALE.The gray scale denotes the materials, with dark gray for the projectile, medium gray for the matrix material, and light gray for the boulders.

Figure A3 .
Figure A3.Cumulative ejecta mass normalized by the projectile mass as a function of the ejecta velocity normalized by the impact velocity for the 2D rubble pile with small surface boulders and regolith cohesion (Y o ) set to 1 Pa as a function of grid resolution.

Figure B2 .
Figure B2.Transient crater for 2D CTH simulations 130 ms postimpact: (a) plots the damage for the weak sandy homogeneous basalt target (Y o = 10 Pa), (b) plots the material for the large surface boulder rubble pile (Y o = 10 Pa), (c) plots the material for the small surface boulder rubble pile (Y o = 10 Pa), (d) plots the damage for the weak granular sandy target (Y o = 1 Pa), (e) plots the material for the large surface boulder rubble pile (Y o = 1 Pa), and (f) plots the material for the small surface boulder rubble pile (Y o = 1 Pa).We plot the damage instead of the material for (a) and (d) so that the transient crater can be clearly visualized, which was not possible when plotting the material in the grid (too much low-density material obscures the crater cavity).

Figure B1 .
Figure B1.iSALE crater snapshots 3 s postimpact for (a) the rubble pile, (b) the 1 m subsurface layer, and (c) the 5 m subsurface layer model.The left side of each image shows the accumulated damage, and the right side shows the density.Craters are still evolving at this time.The cavity along the symmetry axis in panel (b) is a numerical artifact.

Figure B3 .
Figure B3.Images of the damage profile and transient craters in 2D CTH 130 ms postimpact for (a) homogeneous weak sand Y o = 10 Pa, (b) homogeneous weak sand Y o = 1 Pa, (c) rubble pile with regolith Y o = 10 Pa, (d) rubble pile with regolith Y o = 1 Pa, (e) rubble pile with strong core 0.88 m from the surface, (f) rubble pile with strong core 1.23 m from the surface, (g) rubble pile with strong core 2.12 m from the surface, (h) rubble pile with core 5 m from the surface, (i) rubble pile with strong core 10 m from the surface, (j) rubble pile with strong core 15 m from the surface, and (k) rubble pile with weak core (Y o = 1 Pa) 0.88 m from the surface.Note that there was no damage model applied to the weak core material, which is why it appears white in color.20

Table A1 Model
Parameters for Target Materials in Simulations of Idealized Impact on Single Boulders

Table A2 Model
Parameters for Target Material

Table A3 Material
Model Parameters Used for the 2D CTH Simulations