# Constraining surface properties of asteroid (162173) Ryugu from numerical simulations of Hayabusa2 mission impact experiment

Martin Jutzi, Sabina D Raducan, Yun Zhang, Patrick Michel, Masahiko Arakawa

<figure><img src="image_1.png" /></figure>

## To cite this version:

Martin Jutzi, Sabina D Raducan, Yun Zhang, Patrick Michel, Masahiko Arakawa. Constraining surface properties of asteroid (162173) Ryugu from numerical simulations of Hayabusa2 mission impact experiment. Nature Communications, 2022, 13, pp.7134. 10.1038/s41467-022-34540-x. hal-04264150

HAL Id: hal-04264150

https://hal.science/hal-04264150v1

Submitted on 30 Oct 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_3.png" /></figure>

# Constraining surface properties of asteroid (162173) Ryugu from numerical simulations of Hayabusa2 mission impact experiment

Received: 28 March 2022

Accepted: 27 October 2022

Published online: 30 November 2022

<figure><img src="image_4.png" /></figure>

Check for updates

Martin Jutzi Ⓔ<fnref n="1" />✉, Sabina D. Raducan Ⓔ<fnref n="1" />, Yun Zhang Ⓔ²,³, Patrick Michel Ⓔ³ & Masahiko Arakawa Ⓔ⁴

The Hayabusa2 mission impact experiment on asteroid Ryugu created an unexpectedly large crater. The associated regime of low-gravity, low-strength cratering remained largely unexplored so far, because these impact conditions cannot be re-created in laboratory experiments on Earth. Here we show that the target cohesion may be very low and the impact probably occurred in the transitional cratering regime, between strength and gravity. For such conditions, our numerical simulations are able to reproduce the outcome of the impact on Ryugu, including the effects of boulders originally located near the impact point. Consistent with most recent analysis of Ryugu and Bennu, cratering scaling-laws derived from our results suggest that surfaces of small asteroids must be very young. However, our results also show that the cratering efficiency can be strongly affected by the presence of a very small amount of cohesion. Consequently, the varying ages of different geological surface units on Ryugu may be due to the influence of cohesion.

Both the projectile's properties (e.g., mass, velocity), and the target's properties (e.g., strength, gravity, porosity, structure) determine the size and morphology of impact craters on a planetary surface. Crater sizes and morphologies on the surface of asteroids can provide a direct diagnosis of their material properties and their near-surface structure. Studies of the cratering process have important implications for understanding the chronology of the geological and geophysical evolution of asteroids, as well as for the design of planetary defence missions aimed at deflecting them by a kinetic impact.

The crater size is known to be controlled by gravity, strength, or a combination of the two, but it can also be controlled by porosity in the case of compaction-controlled craters<fnref n="1" />. Small asteroids like Ryugu have a very low gravity, but are also very weak, making it ambiguous which of these parameters controls the crater growth. While some mechanical properties of the surface of Ryugu can be inferred from the Hayabusa2 data, one large uncertainty is what is the surface strength. For the same impact energy, stronger surfaces are expected to have smaller craters than weaker surfaces.

The small carry-on impactor (SCI) experiment on the surface of Ryugu had known impact conditions (i.e., projectile mass, impact velocity) and represents an ideal test case to validate impact cratering models in the actual conditions of an impact onto an asteroid. Furthermore, it provides important insights to the properties, origin and evolution of the small body populations. The analysis of the cratering process induced by the recent SCI experiment on asteroid Ryugu concluded that the impact might have taken place in the gravity-dominated regime² and therefore the surface cohesion on Ryugu must be very small (<-1 Pa). These conclusions were based on the extrapolation of laboratory-derived scaling relationships applied to the surprisingly large SCI crater size (with a diameter of about 14.5 m³), as well as on the behaviour of the ejecta curtain. However, because the Hayabusa2 impact experiment was the first of its kind, it is unclear whether traditional impact models and crater scaling-laws are applicable to the low-strength low-gravity cratering environment. How much the cratering process and its outcome, such as the crater size and morphology, were influenced by target surface inhomogeneities (i.e.,

<footnote n="1">¹Space Research and Planetary Sciences, Physikalisches Institut, University of Bern, Bern, Switzerland. ²Department of Aerospace Engineering, University of Maryland, College Park, USA. ³Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Nice, France. ⁴Graduate School of Science, Kobe University, Kobe 657-8501, Japan. ✉e-mail: martin.jutzi@unibe.ch</footnote>
the boulders located close to the impact point) is another important question.

Using shock-physics codes, it has not been feasible so far to model the entire crater formation in the gravity-regime on small (100–1000 m) asteroids, because of the vastly different time-scales of the shock-wave propagation and the crater formation.

Here we use the Bern Smoothed-Particle-Hydrodynamics (SPH) impact code$^{3–5}$ to study the sensitivity of the SCI impact outcome to target cohesion and target inhomogeneities (i.e., boulders). The late-stage evolution of the crater is simulated using a recently developed fast integration scheme. Using this approach, the outcome of the SCI impact on Ryugu can be well reproduced. Our simulations show that the cratering efficiency in the investigated low-strength low-gravity cratering environment is very large, but can be strongly affected by the presence of a small amount of cohesion.

## Results

The shock physics code Bern SPH includes material models suitable to simulate the behaviour of geological materials, various equations of state and a porosity compaction model. We simulate the late-stage evolution by applying the numerical approach developed in$^6$ as it allows for faster calculation times and the simulation of the entire process (see Methods section, subsection Modelling approach for the late stage evolution). This is achieved via application of a transition to a low-speed medium at some advanced time of the simulation. In fact, after the initial shock has passed, the late-stage evolution is governed by low-velocity granular flows, which can be accurately modelled using a low bulk sound speed material (i.e., a material with a small $c_s$, allowing for a larger time step $dt$). This approach has been extensively tested and validated against laboratory experiments into homogeneous and heterogeneous targets$^7$. Furthermore, the comparison with an alternative, SPH-SSDEM modelling scheme, in which the late stage crater evolution is modelled using the soft-sphere discrete-element method (SSDEM) to directly simulate the granular material dynamics, shows a very good agreement (see Methods section, subsection Modelling approach for the late stage evolution).

The target material response to shear deformation is described by a simple pressure-dependent strength model$^{8,9}$. A measure of strength is given by the material's ability to withstand different types of stress states. A granular material, for example, can still have a significant amount of shear strength originating from van der Waals forces and the inability of the interlocking particles to move apart and slide over one another$^{10}$. Here we study the effects of the shear strength at zero pressure, which is often referred to as cohesion. Since the surface strength is the unknown quantity in the study, we first varied the initial cohesion ($Y_0 = 0–0.5$ Pa) and coefficient of friction ($f = 0.4–1.0$) (see Methods section, subsection Shock physics code model). For our nominal simulations, we keep the target initial density constant at 1300 kg/m$^3$ (initial porosity $\phi_0 = 50\%$), which is slightly higher than the bulk density of Ryugu (1190 kg/m$^3$) measured by Hayabusa2$^{11}$. We include 50%, microporosity (as described in refs. 4, 12), which is consistent with recent analysis of Hayabusa2 samples returned from Ryugu$^{13}$.

Our simulations show that the resulting final crater is very sensitive to our choice of cohesion and coefficient of internal friction (Fig. 1). The solution is not unique and the same crater size can be produced by impacting targets with a combination of mechanical properties (represented by the friction coefficient and the cohesion). In order to match the observed size of the SCI crater, a cohesionless ($Y_0 = 0$ Pa) target with the nominal porosity of 50% requires a relatively high friction coefficient ($f > 0.8$). On the other hand, the required crater size is also produced with a cohesion of -0.2 Pa and a small friction coefficient ($f = 0.4$).

It is not possible to constrain the cohesion at the required level, i.e., to distinguish between 0.0 Pa and 0.2 Pa cohesion, based on the

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 1** | Crater radius as a function of material properties. a Crater radius as a function of cohesion $Y_0$ for different friction coefficients ($f = 0.4$ and $f = 0.55$) and porosities (30, 50, 60 %). b Crater radius as a function of the coefficient of internal friction for 50% porous, cohesionless targets ($Y_0 = 0$). The horizontal line shows the measured SCI crater size at the pre-impact level (7.25 ± 0.8 m), with its associated error$^2$. Source data are provided as a Source Data file.</figcaption></figure>

available observations. However, reasonable assumptions can be made regarding the friction coefficient. The angle of repose of cohesionless materials has been measured to be 22° for glass beads, e.g. refs. 14, 15 and 30° for quartz sand, e.g. refs. 15, 16. For lunar regolith simulant with a small but significant cohesion in the range of 0.5–2 kPa, the angle of repose is significantly higher, in the range of 35–45° e.g. refs. 15, 17. The corresponding friction coefficients are about $f = 0.4$ (glass beads), $f = 0.55$ (quartz sand) and $f = 0.7–0.9$ (regolith simulant; Methods section, subsection Shock physics code model). We consider cohesionless quartz sand as the best representative of Ryugu's subsurface material because of the evidence for a very low cohesion (<1 Pa), and we define $f = 0.55$ as nominal friction coefficient. The friction coefficient for glass beads is regarded as a lower limit for geological materials and our results for targets with $f = 0.4$ and 50% porosity suggest that the cohesion present on Ryugu at the SCI impact site cannot be higher than -0.2 Pa (Fig. 1a).

However, in addition to friction and cohesion, the target porosity is also an important parameter. e.g., ref. 1. The relation between porosity and cohesion depends on the scale and type of porosity. For instance, a target made of (porous) boulders or small particles can have a relatively high porosity, but an overall small bulk cohesion (contact forces between the particles). Such targets would have a bulk crushing strength much larger than the bulk cohesion.

For a surface with a very low bulk cohesion ($\lesssim$1 Pa) and significant crushing strength ($>>$1 Pa), the SCI impact is not expected to be in the so-called compaction cratering regime (in which the crater forms mostly by crushing and permanent compaction of the target$^1$). Our results also confirm this, even for highly porous targets (i.e., the amount of compacted material in the final crater is negligible). This is because of the relatively small scale of the SCI impact. In craters formed at small scales, outward expansion caused by the late-time material flow dominates over compaction$^1$. The crater would form in the compaction cratering regime only if the compressive strength was of the same order as the cohesion ($\lesssim$1 Pa), while weak soils have
compressive strength > 10⁴ Pa¹⁸. However, porosity is still relevant for the crater size, as it affects the coupling to the target and it also determines the bulk density, which is an important parameter for the cratering efficiency in both the strength and gravity dominated cratering regimes, e.g., ref. 19. In order to investigate the effect of porosity and to obtain an upper limit of cohesion on Ryugu's surface, we perform impact simulations into targets with f = 0.4, varying cohesion (Y₀ = 0-0.5 Pa) and varying initial porosity (φ₀ = 30-60%). Including these porosity variations, our results suggest an overall upper limit for the cohesion on Ryugu's surface of ~0.5 Pa (Fig. 1a).

For a subset of our simulations, boulders with varying sizes are distributed within the target (Fig. 2), in a similar way as that observed at the SCI impact zone. The boulders have the same density as the matrix

<figure><img src="image_7.png" /><figcaption>**Fig. 2** | End-to-end SPH simulation of the SCI impact. a-c Snapshots of the simulation are shown at different times. At t = 1200 s the crater evolution is finished. d Observed SCI crater on Ryugu². The local height is measured relative to point F. The crater dimensions are also indicated by the dotted orange (crater wall) and dashed black (crater rim) lines in c. The overall outcome in terms of crater size and boulder displacement is well reproduced. In the simulation, the largest block (4) is moved by 10-20 cm and the medium boulder (3) by 2 m, which is roughly consistent with the observed displacement of the large boulders at the SCI impact site (36 cm for the large Okamoto boulder; 2.58 m for the medium sized lijima boulder;²⁰). In the simulation, the small boulders (1,2) are ejected out of the crater, as observed in the case of the SCI impact for some of the boulders (e.g. boulders mb01 and mb02; ref. 20).</figcaption></figure>

material but a cohesion of 100 MPa. The boulders are assigned a high strength because we are interested in the movement of intact boulders rather than their fragmentation (fragmentation has not been observed in the SCI impact). The homogeneous matrix material is considered to be cohesionless (Y₀ = 0 Pa) and with a coefficient of friction f = 0.55. The sizes and positions of these boulders (see Methods, subsection Boulder position and impact geometry) are representative of the large, medium and small boulders present close to the SCI impact point²⁰ on Ryugu.

The overall outcome of the SCI impact, including the displacement of the boulders (Fig. 2), is well reproduced in our simulations. The presence of the large block leads to an asymmetric crater similar to the observed one. The comparison to a simulation without boulders but otherwise the same target properties (Fig. 3) shows that the presence of the boulders within the target affects the resulting crater size only by a minimal amount (<5%). In the simulations, the transient crater is reached at about t ≈ 200-300 s (Fig. 2), which is consistent with the observations².

Estimates of an asteroid's surface age rely on the relationship between the projectile's size and the crater's size, which can be derived from the crater scaling-laws. The general form of the scaling relationship for the crater radius, R, depends on velocity, U; radius, a; density, δ; mass, m of the impactor and the target properties: density, ρ; strength, Y, and gravitational acceleration, g. It is given by, e.g., refs. 19, 21:

$$R\left(\frac{\rho}{m}\right)^{1/3} = K_{R1}\left\{\frac{ga}{U^2}\left(\frac{\rho}{\delta}\right)^{\frac{(6\mu-2)\mu}{3\mu}}\right. \\ \left. + \left[K_{R2}\left(\frac{Y}{\rho U^2}\right)\left(\frac{\rho}{\delta}\right)^{\frac{(6\mu-2)}{3\mu}}\right]^{\frac{2+\mu}{2}}\right\}^{\frac{\mu}{2\mu}} \quad (1)$$

where the constants μ, ν, K<sub>R1</sub> and K<sub>R2</sub> are determined empirically e.g.²¹.

Using the nominal coefficient of friction (f = 0.55), we investigate the dependence of the impact outcome to the impactor size for different values of cohesion. The simulation results are then used to calibrate Eq. (1), which yields the fitting constants K<sub>1</sub> = 0.60, K<sub>2</sub> = 0.25, μ = 0.44 and Y = Y₀ = 0.05 Pa.

Our analysis shows that even in this low-strength low-gravity regime – which is not accessible by laboratory experiments for the size and timescales relevant to asteroid cratering – the general form of the impact cratering scaling (Eq. (1)) allows for a very good fit of the crater sizes as a function of the impactor size and target strength (i.e., cohesion). Our results also suggest that a cohesion of 0.05 Pa is required to match the SCI crater size for nominal material parameters. This means that the SCI impact may have taken place at the transition

<figure><img src="image_8.png" /><figcaption>**Fig. 3** | Comparison of SCI impact simulations with varying target surface properties. a Targets with boulders and b targets without boulders, with otherwise the same impact conditions. In both cases, the homogeneous target material is</figcaption></figure>

cohesionless (while the boulders have significant cohesion). The case shown in a corresponds to the one shown in Fig. 1. The crater radius in the case of the homogeneous target is 8.82 m.
<figure><img src="image_9.png" /><figcaption>**Fig. 4** | Impact outcome as a function of impact size. Shown is the scaled crater radius as a function of gravity scaled impact size for different target cohesions. Shown are the results of the simulations using a nominal coefficient of friction of f = 0.55. The scaling relationships for coarse and perlite sand (H&H2011; ref. 44) are shown for comparison. The case with a cohesion of Y = 0.05 Pa leads to the best fit with the observed SCI crater. Even small increases of the cohesion lead to significantly smaller craters. The background gradient shows the transition between the strength and gravity dominated cratering regimes. Source data are provided as a Source Data file.</figcaption></figure>

between the strength and gravity regime (Fig. 4), where cohesive and gravitational forces determine the impact outcome concurrently.

Discussion

The cratering efficiency defined as $\Pi_R = R(\frac{\rho}{m})^{1/3} = D_{crat}/D_{proj} (\frac{3\rho}{4\pi\delta})^{1/3}$ (left term in Eq. (1)), where $D_{crat}$ is the crater diameter and $D_{proj}$ is the projectile diameter, largely determines the surface ages of asteroids based on crater counting.

As shown by the SCI experiment and confirmed by the results of the study presented here, small-scale impacts in low-cohesion targets have a very large cratering efficiency, with a correspondingly large ratio of $D_{crat}/D_{proj}$ of > 100 (Fig. 4; with $D_{crat}/D_{proj} = \Pi_R \times (\frac{3\rho}{4\pi\delta})^{-1/3} \approx \Pi_R \times 2.1$). This is in contrast with the assumption of a (constant) small ratio of $D_{crat}/D_{proj} = 10$ that has been suggested to be valid over a larger range of scales$^{22}$. Important consequences of such a high cratering efficiency are high resurfacing rates and young crater retention ages. For example, a ratio $D_{crat}/D_{proj}$ of about 100 leads to an about $10^3$ times higher crater formation rate (assuming $q = -4$ for the slope of the impactor population, see Methods, subsection Calculations of the impact frequencies) than a ratio of $D_{crat}/D_{proj} = 10$ and a correspondingly younger surface age. There is increasing evidence that the surfaces of both recently visited asteroids Ryugu and Bennu are very weak and our results suggest that their surface must therefore be very young (of the order of -1-10 Myr). These findings are consistent with most recent analysis, e.g., refs. 23-25, but in contrast to previous studies which proposed much older surface ages (-1 Gyr) based on low cratering efficiencies, corresponding to strength regime scaling, e.g., refs. 22, 26.

Assuming asteroid Ryugu's cohesion is negligible and crater formation is controlled by the gravity$^{23}$, showed that the absolute model ages of asteroid Ryugu's geologic units are in about 2-30 Myrs range. As also indicated by the observations, there are variations in surface properties (e.g., strength), which possibly lead to varying ages for the different units. Such differences in surface ages might indeed be explained by small variations in cohesion, which lead to considerable changes of the cratering efficiencies resulting from the crater scaling-laws used here. As an example, our best fit scaling-law (Fig. 4; Y = 0.05 Pa) leads to $D_{crat}/D_{proj}$ of about 180 for an impactor size of

$a = 0.1$ m, which in turns leads to a crater radius of about 18 m (assuming an impact velocity of 5 km/s). Using a slightly higher cohesion of 1 Pa, we obtain $D_{crat}/D_{proj}$ of about 110 and a crater radius of about 11 m. Following the procedure used in$^{27}$, this leads to an approximately 4 times higher rate of crater formation in this size range for the Y = 0.05 Pa surface compared to the Y = 1 Pa surface (see Methods section, subsection Calculations of the impact frequencies) and a correspondingly younger surface age. As illustrated by this example, even a small amount of cohesion leads to a significantly reduced cratering efficiency. The influence of cohesion becomes larger with decreasing size scale (Fig. 4) and leads to a crater reduction effect for small-scale impacts on low-strength asteroid surfaces. Recent studies have shown that on weak asteroids, resurfacing by seismic shaking is not efficient$^{28}$. On the other hand, armouring by surface boulders has been proposed to explain the paucity of craters smaller than 2-3 m on the asteroid Bennu$^{24}$. Here, we demonstrate that the presence of a small amount of cohesion can actually contribute significantly to the crater reduction effect and therefore to the depletion of small craters in the observed crater distribution on small asteroids e.g., refs. 28, 29.

Methods

Shock physics code model

We modelled the target using the Tillotson equation of state (EoS) for basalt$^{30}$, with a modified initial bulk modulus (see next section) for the simulations of the late stage evolution. For the target material response to shear deformation we applied a simple pressure-dependent strength model, typical of geological materials$^{8,9}$. The Lundborg strength model describes the yield strength as:

$$Y = Y_0 + \frac{fP}{1+fP/(Y_{dm} - Y_0)}, \quad (2)$$

where P is pressure, f is the coefficient of internal friction, $Y_{dm}$ is the limiting strength at high pressure and $Y_0$ is the cohesion.

We use a constant cohesion, $Y_0$, for the weak asteroid materials considered in this study and apply a strain-based weakening model that prevents artificial clumping (similar to the approach used in$^{31}$). Our model uses a linear relation between cohesion $Y_0$ and total strain $\epsilon_{tot}$, and it is assumed that for $\epsilon_{tot} \ge 1$ cohesion is lost. Instead of using a tensile fragmentation model, like it was implemented in previous studies e.g.,$^3$ the tensile strength is defined by extrapolating the yield strength (versus pressure) curve (Eq. (2)) to intersect the pressure axis. In addition, the pressure is limited to a minimum value $P_{min} \ge Y_0$.

To compute the coefficient of internal friction from the measured angle of response $\theta^{15,17}$, we follow the relation described in ref. 32:

$$f = \frac{2\sqrt{2\sin(\theta)}}{3 - \sin(\theta)} \quad (3)$$

The porosity was modelled using the P-$\alpha$ model$^4$ with a crush curve defined by the parameters $P_s$, $P_e$, $P_t$, $\alpha_t$, $n_1$, $n_2$ and initial distention, $\alpha_0$. $P_e$ is the pressure defining the transition between the elastic and plastic regimes, $P_s$ the pressure at which all pores are compacted (distention $\alpha = 1$) and the parameters $P_e < P_t < P_s$ and $1 < \alpha_t < \alpha_0$ indicate a transition between two regimes with different slopes ($n_1$ and $n_2$).

For both target and projectile materials, the Tillotson EOS is applied. For the relatively low SCI impact velocity, no significant vapourisation or melting is expected and the Tillotson EOS is appropriate for this impact regime$^{33}$. The relevant material parameters used in this study are summarised in Table 1.

Initial conditions

The target is modelled as a hemispherical domain with a radius of 15 m (Fig. 3) using approximately $10^7$ SPH particles. We use a spherical
<table><caption>**Table 1 | Material model parameters for impact simulations into Ryugu analogues**</caption><thead><tr><th>Description</th><th>Impactor</th><th>Homogeneous material</th><th>Boulders</th></tr></thead><tbody><tr><td>Material</td><td>Aluminium</td><td>Basalt</td><td>Basalt</td></tr><tr><td>Impact angle (°)</td><td>90</td><td>–</td><td>–</td></tr><tr><td>Impact speed (km/s)</td><td>2</td><td>–</td><td>–</td></tr><tr><td>Equation of state</td><td>Tillotson<sup>a</sup></td><td>Tillotson<sup>b</sup></td><td>Tillotson<sup>b</sup></td></tr><tr><td>Strength model</td><td>von Mises</td><td>LUND<sup>c</sup></td><td>LUND<sup>c</sup></td></tr><tr><td>LUND strength parameters<sup>c</sup></td><td>–</td><td>–</td><td>–</td></tr><tr><td>Damage strength at zero pressure, Y<sub>0</sub> (Pa)</td><td>–</td><td>0/0.05/0.1/0.5</td><td>1 × 10<sup>8</sup></td></tr><tr><td>Strength at infinite pressure, Y<sub>inf</sub> (GPa)</td><td>0.34</td><td>3.5</td><td>3.5</td></tr><tr><td>Internal friction coefficient (damaged), f</td><td>–</td><td>0.4/0.55/0.8/1.0</td><td>0.8</td></tr><tr><td>Porosity model parameters<sup>d</sup></td><td>–</td><td>–</td><td>–</td></tr><tr><td>Initial porosity, φ<sub>0</sub></td><td>–</td><td>50%</td><td>–</td></tr><tr><td>Initial distension, α<sub>0</sub></td><td>–</td><td>2</td><td>–</td></tr><tr><td>Transition distention, α<sub>t</sub></td><td>–</td><td>1.25</td><td>–</td></tr><tr><td>Pressure at full compaction, P<sub>s</sub> (MPa)</td><td>–</td><td>213</td><td>–</td></tr><tr><td>Transition pressure, P<sub>t</sub> (MPa)</td><td>–</td><td>68</td><td>–</td></tr><tr><td>Elastic pressure, P<sub>e</sub> (MPa)</td><td>–</td><td>1.0</td><td>–</td></tr><tr><td>Slope regime 1, n<sub>1</sub></td><td>–</td><td>12</td><td>–</td></tr><tr><td>Slope regime 2, n<sub>2</sub></td><td>–</td><td>3</td><td>–</td></tr></tbody></table>

projectile (Table 1) with a mass of 2 kg, impacting head-on. Recent studies³³ have shown that for low-strength targets with (with associated high cratering efficiencies), the projectile shape does not significantly influence the crater size. Because of the large simulation domain required to model the impacts with very high cratering efficiencies, the projectile shape cannot be resolved in our simulations. We approximate the SCI impactor using a homogeneous sphere with a the same diameter as the SCI at impact but with a correspondingly lower density (roughly corresponding to Aluminium), following the approach used in previous studies, e.g., refs. 33–35.

For the simulations with varying projectiles (Fig. 4), we use projectile masses of 0.05, 2, 100 and 2000 kg, respectively, and adapt the size of the computational domain and the computation times accordingly.

All models used a constant gravity acceleration of $1.4 \times 10^{-4} \text{ m/s}^2$.

## Modelling approach for the late stage evolution

Modelling the entire crater formation in the gravity-regime on small (-10–1000 m) asteroids is challenging because of the vastly different time-scales of the shock-wave propagation and the crater formation. The first time-scale is governed by the elastic wave velocity $c_s$, typically a few km/s for rocks. For example, to model a SCI-like impact and to ensure numerical stability, the so-called Courant criteria requires that the timestep, $dt$, is smaller than the simulation resolution divided by the sound speed in the target, $c_s$: $dt < \text{resolution}/c_s \approx 10^{-6} - 10^{-7} \text{ s}$. On the other hand, the crater formation time in the gravity regime can last up to few hundred seconds¹.

**Direct shock physics code modelling of late stage evolution.** To model the late-stage evolution of the crater grow we follow the approach introduced in⁶ and apply in the shock physics code

<figure><img src="image_10.png" /><figcaption>**Fig. 5 | Comparison of two numerical approaches.** Shown are the results of direct SPH (this work) and SPH-SSDEM simulations³⁶. The crater profile obtained by the hybrid SPH-SSDEM simulation (yellow line) is over-plotted on the results of the direct SPH simulation. Indicated in grey is the uncertainty due to the SPH kernel interpolation and associated smoothing length. The final crater size and morphology are in good agreement.</figcaption></figure>

calculation a transition to a low-sound-speed medium. After the initial shock has passed, the late-stage evolution is governed by low-velocity granular flow, which can be accurately modelled using a low bulk sound speed material (i.e., a material with a small $c_s$, allowing for a larger $dt$).

At this time $t_{transition}$ we switch to a low-sound-speed medium (i.e., a fast scheme), which allows us to model the late-stage evolution up to about a thousand seconds after the impact. In this step, we apply a simplified Tillotson EoS, in which all energy related terms are set to zero. The remaining leading term of the EoS is governed by the bulk modulus $P = A(1 - \rho/\rho_0)$, which also determines the magnitude of the sound speed. We use $A \approx 0.027 \text{ MPa}$ and also reduce the shear modulus proportionally.

This approach has been extensively tested and validated against laboratory experiments into homogeneous and heterogeneous targets⁷.

The transition time $t_{transition}$ depends on the scale of the impact. We use the following definition $t_{transition} \approx L/v_{sound} \times 10$ where $L$ is the length scale (chosen to be the size of the computational domain) and $v_{sound}$ the sound speed of the material. According to this definition, $t_{transition}$ corresponds to roughly 10 sound crossing times. For the SCI-scale impact, we use $t_{transition} = 0.1 \text{ s}$.

To test the robustness of the approach, we performed additional simulations with different transition times. For this test, we use the larger scale impactor (10 × SCI-size; Fig. 4) as the simulations are faster to perform at larger scales. We use transition times of $t_{transition} = 0.5, 0.8, 1.4$ and $2.0 \text{ s}$. For these different transition times, the resulting crater radii of 52.0, 52.8, 52.8 and 53.6 m agree very well with each other; the deviations are within the uncertainty of the crater size determination (± 1.1 m).

**SPH-SSDEM framework for late stage evolution.** An alternative approach to direct SPH simulations is to use a SPH-SSDEM framework³⁶, in which the late stage crater evolution is modelled using the soft-sphere discrete-element method (SSDEM). The SSDEM method can be used to model low-velocity collisions and the flow of granular materials and it is well suited to simulate the late-stage of the impact process. For this study, the parallelised hierarchical tree SSDEM code, pkdgrav, is used. To compute particle contact interaction and control the material shear strength, we used a granular physics model with four dissipation/friction components in the normal, tangential, rolling and twisting directions³⁷,³⁸.

The combination of the SPH and SSDEM methods has been used in previous studies to simulate small body impact processes e.g. ³⁹,⁴⁰. The initial shock propagation and fragmentation stage is simulated using SPH and the outcome is then transferred into the SSDEM code, which computes that late-stage evolution. The particle-based description of the two methods allows direct information transition at the particle level. However, there are several issues with this

<footnote><sup>a</sup>ref. 46.</footnote>

<footnote><sup>b</sup>ref. 30.</footnote>

<footnote><sup>c</sup>ref. 8.</footnote>

<footnote><sup>d</sup>ref. 4.</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 6 | Boulder positions and impact geometry.** The impact point is in the centre, marked by the red dot. The location and dimension of the boulders (1–4) are given in Table 2.</figcaption></figure>

direct particle-particle transition, for instance the handling of particle overlaps and maintaining linear and angular momentum constant⁴⁰. Additionally, because SPH simulations use very high particle resolution for properly solving the shock wave propagation and energy dissipation, the direct transition with comparable particle resolution in the subsequent SSDEM simulations to complete the crater growth modelling would require computational times that are not achievable⁴¹.

We use a shape-construction algorithm⁴⁰ for the velocity field SPH-SSDEM transition procedure to solve these issues. The procedure consists of four steps: (1) we prepare a granular bed with a predefined particle size distribution in which material parameters are consistent with the SPH simulation, and let the particles settle down in the SSDEM boundary geometry and gravity conditions, (2) based on the given SPH output, we use the α-shape-construction algorithm to construct a surface isolating the compact material from the fast-moving ejecta whose speeds exceed a given limit; (3) we then use the isolated surface to carve out the SSDEM granular bed surface, and conduct a nearest neighbour search to map the linear and angular momentum from the SPH compact material data to the SSDEM bed; (4) finally, we add fast-moving ejecta into the SSDEM simulation scenario as individual particles.

Initial simulations using the SPH-SSDEM framework show a very good agreement with the outcome of the direct SPH shock physics code calculation (Fig. 5).

## Boulder position and impact geometry

The positions of the boulders with respect to the impact point are shown in Fig. 6. The boulders have the same density as the matrix material; their dimensions are given in Table 2. In order to model the large subsurface Okamoto boulder (Fig. 2)²⁰, we use a rectangular block with the dimensions and position indicated in Table 2. The left side of this block (Fig. 6) has a rounded shape to avoid interference with the boundary of the overall (hemispherical) computational domain.

## Crater size computation

In order to determine the size of the final craters, we first compute the continuum density distribution using the SPH kernel

<table><caption>**Table 2 | Position and properties of boulders in the simulations with heterogeneous targets**</caption><thead><tr><td>Boulder number</td><td>Position (x, y, z) (cm)</td><td>Radius (cm)</td></tr></thead><tbody><tr><td>1</td><td>250, 0, 0</td><td>100</td></tr><tr><td>2</td><td>50, -200, 0</td><td>100</td></tr><tr><td>3</td><td>100, 350, 0</td><td>250</td></tr><tr><td>4</td><td>√x² + y² < 900</td><td>-</td></tr><tr><td>-</td><td>x < -200</td><td>-</td></tr><tr><td>-</td><td>|y| < 500</td><td>-</td></tr><tr><td>-</td><td>-500 < z < 50</td><td>-</td></tr></tbody></table>

interpolation. We then define the crater radius as the point where the density at the level of the pre-impact surface falls below a critical value $ρ_{crit} = ρ_{initial}/2$, where $ρ_{initial}$ is the initial bulk density of the homogeneous material. The final crater radius is determined as an average of the radii determined in this way along the x-axis and the y-axis.

## Crater scaling relationships

The crater scaling relationships⁴² relate the outcome of an impact to the impact properties (e.g., impactor radius $a$, impactor density $δ$ and impact speed $U$) by a so-called coupling parameter, which is given by $C ≈ aδU^μ$. Here, $μ$ and $ν$ are target material specific exponents. The density scaling exponent is often assumed to be $ν ≈ 0.4$⁴³. The velocity scaling exponent, $μ$, is one of the main scaling constants needed to extrapolate lab-scale and numerical results to other regimes of applicability e.g., refs. 34, 44, 45. $μ$ depends on the target material properties and its value lies between two theoretical limits: $μ = 1/3$ if the crater formation is influenced by the impactor momentum alone and $μ = 2/3$ if it is influenced by the impactor energy alone.

The most general form of the scaling relationship for crater radius, $R$, as a function of impactor properties and assuming vertical impact of a spherical impactor, is given by Eq. (1). According to this equation, the scaled crater radius is proportional to

$$R\left(\frac{\rho}{m}\right)^{1/3} = \frac{R}{a}\left(\frac{3\rho}{4\pi\delta}\right)^{1/3} = D_{crat}/D_{proj}\left(\frac{3\rho}{4\pi\delta}\right)^{1/3} \quad (4)$$

From Eq. (1) it can be seen that $D_{crat}/D_{proj}$ is constant (for varying $a$) only if $ga < < Y/ρ$, i.e., in the strength regime. This means that using a constant $D_{crat}/D_{proj} = 10^{22}$ implicitly assumes that the cratering outcome is given by purely strength regime scaling, independently of the scale of the event and the actual material strength.

## Calculations of the impact frequencies

The relative impact frequencies are computed by assuming an impactor population described by a differential size distribution with single slope $q^{27}$. We do not compute absolute numbers but rather the ratios resulting from different scaling-law parameters. Based on Eq. (1), the impactor size for a given crater size is obtained which is then used as input to compute the number of expected impacts, following²⁷.

## Data availability

The SPH and pkdgrav simulation data generated in this study are available from the corresponding author upon reasonable request. Source data are provided with this paper.

## Code availability

A compiled version of the Bern SPH code and the pkdgrav code, as well as the necessary input files are available from the corresponding author upon reasonable request. SPH data visualisation was produced using Vapor (www.vapor.ucar.edu).
## References

1. Housen, K. R., Sweet, W. J. & Holsapple, K. A. Impacts into porous asteroids. *Icarus* **300**, 72–96 (2018).

2. Arakawa, M. et al. An artificial impact on the asteroid (162173) ryugu formed a crater in the gravity-dominated regime. *Science* **368**, 67–71 (2020).

3. Benz, W. & Asphaug, E. Simulations of brittle solids using smooth particle hydrodynamics. *Comput. Phys. Commun.* **87**, 253–265 (1995).

4. Jutzi, M., Benz, W. & Michel, P. Numerical simulations of impacts involving porous bodies: I. Implementing sub-resolution porosity in a 3D SPH Hydrocode. *Icarus* **198**, 242–255 (2008).

5. Jutzi, M. SPH calculations of asteroid disruptions: the role of pressure dependent failure models. *Planet. Space Sci* **107**, 3–9 (2015).

6. Raducan, S. D. & Jutzi, M. Global-scale reshaping and resurfacing of asteroids by small-scale impacts, with applications to the DART and hera missions. *Planet. Sci. J.* **3**, 128 (2022).

7. Ormö, J. et al. Boulder exhumation and segregation by impacts on rubble-pile asteroids. *Earth Planet. Sci. Lett.* **594**, 117713 (2022).

8. Lundborg, N. The strength-size relation of granite. *Int. J. Rock Mech. Min. Sci. Geomech. Abstracts* **4**, 269–272 (1967).

9. Collins, G. S., Melosh, H. J. & Ivanov, B. A. Modeling damage and deformation in impact simulations. *Meteorit. Planet. Sci.* **39**, 217–231 (2004).

10. Scheeres, D. J., Hartzell, C. M., Sanchez, P. & Swift, M. Scaling forces to asteroid surfaces: the role of cohesion. *Icarus* **210**, 968–984 (2010).

11. Sugita, S. et al. The geomorphology, color, and thermal properties of Ryugu: implications for parent-body processes. *Science* **364**, 252–252 (2019).

12. Jutzi, M., Michel, P., Benz, W. & Richardson, D. C. Fragment properties at the catastrophic disruption threshold: The effect of the parent body's internal structure. *Icarus* **207**, 54 (2010).

13. Yada, T. et al. Preliminary analysis of the Hayabusa2 samples returned from C-type asteroid Ryugu. *Nat. Astron.* **6**, 214–220 (2021).

14. Lajeunesse, E., Monnier, J. B. & Homsy, G. M. Granular slumping on a horizontal surface. *Phys. Fluids* **17**, 103302–103302–15 (2005).

15. Chourey, S., Koschny, D., Rott, M. & Schmausser, C. Determining the momentum transfer in regolith-like targets using the TUM/LRT electro-thermal accelerator. *Planet. Space Sci.* **194**, 105112(2020).

16. Lube, G., Huppert, H. E., Sparks, R. S. J. & Hallworth, M. A. Axisymmetric collapses of granular columns. *J. Fluid Mech.* **508**, 175–199 (2004).

17. Mitchell, J. K. et al. Mechanical properties of lunar soil: density, porosity, cohesion and angle of internal friction. *Lunar Planet. Sci. Conf. Proc.* **3**, 3235 (1972).

18. Housen, K. R. & Holsapple, K. A. Impact cratering on porous asteroids. *Icarus* **163**, 102–119 (2003).

19. Holsapple, K. A. The scaling of impact processes in planetary sciences. *Annu. Rev. Earth Planet. Sci.* **21**, 333–373 (1993).

20. Honda, R. et al. Resurfacing processes on asteroid (162173) Ryugu caused by an artificial impact of Hayabusa2's Small Carry-on Impactor. *Icarus* **366**, 114530 (2021).

21. Raducan, S. D., Davison, T. M. & Collins, G. S. Morphological diversity of impact craters on asteroid (16) psyche: insight from numerical models. *J. Geophys. Res.: Planets* **125**, e2020JE006466 (2020).

22. Bottke, W. F. et al. Interpreting the cratering histories of bennu, ryugu, and other spacecraft-explored asteroids. *Astron. J.* **160**, 14 (2020).

23. Cho, Y. et al. Geologic history and crater morphology of asteroid (162173) ryugu crater morphology and chronology of asteroid ryugu. *J. Geophys. Res. Planets* **126**, e2020JE006572 (2021).

24. Bierhaus, E. B. et al. Crater population on asteroid (101955) Bennu indicates impact armouring and a young surface. *Nat. Geosci.* **15**, 440–446 (2022).

25. Perry, M. E. et al. Low surface strength of the asteroid Bennu inferred from impact ejecta deposit. *Nat. Geosci.* **15**, 447–452 (2022).

26. Walsh, K. J. et al. Craters, boulders and regolith of (101955) Bennu indicative of an old and dynamic surface. *Nat. Geosci.* **12**, 242 (2019).

27. Jutzi, M. The shape and structure of small asteroids as a result of sub-catastrophic collisions. *Planet. Space Sci.* **177**, 104695 (2019).

28. Takaki, N. et al. Resurfacing processes constrained by crater distribution on Ryugu. *Icarus* **377**, 114911 (2022).

29. Tatsumi, E. & Sugita, S. Cratering efficiency on coarse-grain targets: Implications for the dynamical evolution of asteroid 25143 Itokawa. *Icarus* **300**, 227–248 (2018).

30. Benz, W. & Asphaug, E. Catastrophic disruptions revisited. *Icarus* **142**, 5–20 (1999).

31. Collins, G. S., Kenkmann, T., Osinski, G. R. & Wünnemann, K. Mid-sized complex crater formation in mixed crystalline-sedimentary targets: Insight from modeling and observation. *Meteorit. Planet. Sci.* **43**, 1955–1977 (2008).

32. Luther, R. et al. Momentum enhancement during kinetic impacts in the low-intermediate-strength regime: benchmarking and validation of impact shock physics codes. *Planet. Sci. J.* **3**, 227 (2022).

33. Raducan, S. D. et al. Influence of the projectile geometry on the momentum transfer from a kinetic impactor and implications for the DART mission. *Int. J. Impact Eng.* **162**, 104147 (2022).

34. Raducan, S. D., Davison, T. M., Luther, R. & Collins, G. S. The role of asteroid strength, porosity and internal friction in impact momentum transfer. *Icarus* **329**, 282–295 (2019).

35. Jutzi, M. & Michel, P. Hypervelocity impacts on asteroids and momentum transfer I. Numerical simulations using porous targets. *Icarus* **229**, 247–253 (2014).

36. Zhang, Y., Jutzi, M., Michel, P., Raducan, S. D. & Arakawa, M. A hybrid sph-ssdem framework for end-to-end impact cratering modeling. In *52nd Annual Lunar and Planetary Science Conference*, 1974 (2021).

37. Schwartz, S., Richardson, D. & Michel, P. An implementation of the soft-sphere discrete element method in a high-performance parallel gravity tree-code. *Granul. Matter* **14**, 363–380 (2012).

38. Zhang, Y. et al. Creep stability of the proposed AIDA mission target 65803 Didymos: I. Discrete cohesionless granular physics model. *Icarus* **294**, 98–123 (2017).

39. Schwartz, S. R. et al. Catastrophic disruptions as the origin of bilobate comets. *Nat. Astron.* **2**, 379–382 (2018).

40. Ballouz, R.-L., Walsh, K. J., Richardson, D. C. & Michel, P. Using a geometrical algorithm to provide N-body initial conditions for the gravitational phase of asteroid family formation. *MNRAS* **485**, 697–707 (2019).

41. Schwartz, S. R., Yu, Y., Michel, P. & Jutzi, M. Small-body deflection techniques using spacecraft: Techniques in simulating the fate of ejecta. *Adv. Space Res.* **57**, 1832–1846 (2016).

42. Holsapple, K. A. & Schmidt, R. M. Point source solutions and coupling parameters in cratering mechanics. *J. Geophys. Res.: Solid Earth* **92**, 6350–6376 (1987).

43. Schmidt, R. M. *Meteor Crater: Energy of formation - Implications of centrifuge scaling*. Lunar and Planetary Science Conference. Vol. 11, p. 2099-2128 (Pergamon Press, 1980).

44. Housen, K. R. & Holsapple, K. A. Ejecta from impact craters. *Icarus* **211**, 856–875 (2011).

45. Prieur, N. C. et al. The effect of target properties on transient crater scaling for simple craters. *J. Geophys. Res. Planets* **122**, 1704–1726 (2017).
46. Tillotson, H. J. Metallic equations of state for hypervelocity impact. *Gen. Atomic Rep.* GA-3216, 141 (1962).

**Acknowledgements**

The authors acknowledge funding support from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 870377 (project NEO-MAPP). This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. P.M. Acknowledges funding support from the French space agency CNES and from the CNRS through the MITI inter-disciplinary programmes.

**Author contributions**

M.J. performed and analysed the numerical simulations and led the research. S.D.R. helped to design the numerical study and performed the scaling-law analysis. Y.Z. and P.M. performed and analysed the SSDEM simulations. M.A. provided the SCI impact data. All authors contributed to interpretation of the results and preparation of the manuscript.

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-022-34540-x.

**Correspondence** and requests for materials should be addressed to Martin Jutzi.

**Peer review information** *Nature Communications* thanks Daniel Durda, Daniel Scheeres and the other, anonymous, reviewer for their contribution to the peer review of this work.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022