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

Martin Jutzi  Space Research and Planetary Sciences, Physikalisches Institut,
University of Bern, Switzerland (∗E-mail<EMAIL_ADDRESS> Sabina D.
Raducan  Space Research and Planetary Sciences, Physikalisches Institut,
University of Bern, Switzerland (∗E-mail<EMAIL_ADDRESS> Yun Zhang
University of Maryland, College Park, USA.  Universite Côte d’Azur,
Observatoire de la Côte d’Azur, Laboratoire Lagrange, Nice, France.  Patrick
Michel  Universite Côte d’Azur, Observatoire de la Côte d’Azur, Laboratoire
Lagrange, Nice, France.  Masahiko Arakawa  Graduate School of Science, Kobe
University, Kobe 657-8501, Japan

Abstract

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.

## Introduction

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 [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 [2] and therefore the surface cohesion on Ryugu must
be very small ($<$ approximately 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 [2]), 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., 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, 4,
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, 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, 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, subsection Shock physics
code model). For our nominal simulations, we keep the target initial density
constant at 1300 kg/m3 (initial porosity $\phi_{0}$ = 50%), which is slightly
higher than the bulk density of Ryugu (1190 kg/m3) measured by Hayabusa2 [11].
We include 50%, microporosity (as described in [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 (Figure 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
about 0.2 Pa and a small friction coefficient ($f$ = 0.4).

Figure 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$\pm$0.8 m), with its associated error [2].
Source data are provided as a Source Data file.

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 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. 14, 15] and 30∘ for quartz sand [e.g.
16, 15]. 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. 17, 15]. 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, 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 approximately 0.2 Pa (Figure 1a).

Figure 2: End-to-end SPH simulation of the SCI impact. a) – c) Snapshots of
the simulation are shown at different times. At $t$ = 1200s the crater
evolution is finished. d) Observed SCI crater on Ryugu [2]. 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 Iijima boulder; [18]). 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; [18]). Figure 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 cohesionless (while the boulders have
significant cohesion). The case shown in a) corresponds to the one shown in
Figure 1. The crater radius in the case of the homogeneous target is 8.82 m.

However, in addition to friction and cohesion, the target porosity is also an
important parameter [e.g., 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 ($\lessapprox$ 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
($\lessapprox$ 1 Pa), while weak soils have compressive strength $>$ 104 Pa
[19]. 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. 20]. 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-0.5 Pa) and varying initial porosity ($\phi_{0}$ =
30-60%). Including these porosity variations, our results suggest an overall
upper limit for the cohesion on Ryugu’s surface of approximately 0.5 Pa
(Figure 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 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}$ = 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 [18] 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\approx$ 200-300 s
(Fig. 2), which is consistent with the observations [2].

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, $\delta$; mass,
$m$ of the impactor and the target properties: density, $\rho$; strength, $Y$,
and gravitational acceleration, $g$. It is given by [e.g. 20, 21]:

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

where the constants $\mu$, $\nu$, $K_{R1}$ and $K_{R2}$ are determined
empirically [e.g. 21].

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_{1}$ = 0.60, $K_{2}$ = 0.25, $\mu$ = 0.44 and
$Y=Y_{0}$ = 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 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\left(\frac{\rho}{m}\right)^{1/3}=D_{crat}/D_{proj}\left(\frac{3\rho}{4\pi\delta}\right)^{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\left(\frac{3\rho}{4\pi\delta}\right)^{-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 approximately 1–10 Myr).
These findings are consistent with most recent analysis [e.g. 23, 24, 25], but
in contrast to previous studies which proposed much older surface ages
(approximately 1 Gyr) based on low cratering efficiencies, corresponding to
strength regime scaling [e.g. 26, 22].

Figure 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; 27] 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.

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 [28], 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,
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 [29]. 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., 30, 29].

## Methods

### Shock physics code model

We modelled the target using the Tillotson equation of state (EoS) for basalt
[31], 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})},$ (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 [32]). Our model
uses a linear relation between cohesion $Y_{0}$ and total strain
$\epsilon_{tot}$, and it is assumed that for $\epsilon_{tot}\geq 1$ cohesion
is lost. Instead of using a tensile fragmentation model, like it was
implemented in previous studies [e.g., 5], 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}\geq-Y_{0}$.

To compute the coefficient of internal friction from the measured angle of
response $\theta$ [15, 17], we follow the relation described in [33]:

$f=\frac{2\sqrt{2}sin(\theta)}{3-sin(\theta)}$ (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 distension, $\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 vaporisation or melting
is expected and the Tillotson EOS is appropriate for this impact regime [34].
The relevant material parameters used in this study are summarised in Table 1.

Table 1: Material model parameters for impact simulations into Ryugu analogues. Description | Impactor | Homogeneous material | Boulders
---|---|---|---
Material | Aluminium | Basalt | Basalt
Impact angle (o) | 90 | – | –
Impact speed (km/s) | 2 | – | –
Equation of state |  Tillotsona |  Tillotsonb |  Tillotsonb
Strength model | von Mises |  LUNDc |  LUNDc
LUND strength parametersc |  |  |
Damage strength at zero pressure, Y0 (Pa) | – | 0/0.05/0.1/0.5 | $1\times 10^{8}$
Strength at infinite pressure, Yinf (GPa) | 0.34 | 3.5 | 3.5
Internal friction coefficient (damaged), $f$ | – | 0.4/0.55/0.8/1.0 | 0.8
Porosity model parametersd |  |  |
Initial porosity, $\phi_{0}$ | – | 50% | –
Initial distension, $\alpha_{0}$ | – | 2 | –
Transition distention, $\alpha_{t}$ | – | 1.25 | –
Pressure at full compaction, $P_{s}$ (MPa) | – | 213 | –
Transition pressure, $P_{t}$ (MPa) | – | 68 | –
Elastic pressure, $P_{e}$ (MPa) | – | 1.0 | –
Slope regime 1, $n_{1}$ | – | 12 | –
Slope regime 2, $n_{2}$ | – | 3 | –
a[35]; b[31]; c[8]; d[4]. |

### 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 projectile
(Table 1) with a mass of 2 kg, impacting head-on. Recent studies [34] 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. 36, 34, 37].

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}$ m/s2.

### Modelling approach for the late stage evolution

Modelling the entire crater formation in the gravity-regime on small (about
100 m – 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<\mathrm{resolution}/c_{s}\simeq 10^{-6}-10^{-7}$ s. On
the other hand, the crater formation time in the gravity regime can last up to
few hundred seconds [2].

#### Direct shock physics code modeling of late stage evolution

To model the late-stage evolution of the crater grow we follow the approach
introduced in [6] and apply in the shock physics code 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$ MPa and also reduce the shear modulus
proportionally.

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

The transition time $t_{transition}$ depends on the scale of the impact. We
use the following definition $t_{transition}\simeq 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 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 $\times$ 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 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 are within the uncertainty of the crater size determination
($\pm$ 1.1 m).

#### SPH-SSDEM framework for late stage evolution

An alternative approach to direct SPH simulations is to use a SPH-SSDEM
framework [38], 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 parallelized 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 [39, 40].

The combination of the SPH and SSDEM methods has been used in previous studies
to simulate small body impact processes [e.g. 41, 42]. 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 direct particle-particle transition, for instance the
handling of particle overlaps and maintaining linear and angular momentum
constant [42]. 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 [43].

We use a shape-construction algorithm [42] 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
$\alpha$-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).

Figure 5: Comparison of two numerical approaches. Shown are the results of
direct SPH (this work) and SPH-SSDEM simulations [38]. 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 gray is the uncertainty
due to the SPH kernel interpolation and associated smoothing length. The final
crater size and morphology are in good agreement.

### 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) [18], 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.

Figure 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. Table 2: Position and properties of boulders in the simulations with heterogeneous targets. Boulder number | Position (x, y, z) (cm) | Radius (cm)
---|---|---
1 | 250, 0, 0 | 100
2 | 50, -200, 0 | 100
3 | 100, 350, 0 | 250
4 | $\sqrt{x^{2}+y^{2}}<900$ | –
$x<-200$
$|y|<500$
$-500<z<50$

### Crater size computation

In order to determine the size of the final craters, we first compute the
continuum density distribution using the SPH kernel 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
$\rho_{crit}=\rho_{initial}/2$, where $\rho_{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 [44] relate the outcome of an impact to the
impact properties (e.g., impactor radius $a$, impactor density $\delta$ and
impact speed $U$) by a so-called coupling parameter, which is given by
$C\approx a\delta^{\nu}U^{\mu}$. Here, $\nu$ and $\mu$ are target material
specific exponents. The density scaling exponent is often assumed to be
$\nu\approx 0.4$ [45]. The velocity scaling exponent, $\mu$, is one of the
main scaling constants needed to extrapolate lab-scale and numerical results
to other regimes of applicability [e.g., 27, 46, 36]. $\mu$ depends on the
target material properties and its value lies between two theoretical limits:
$\mu=1/3$ if the crater formation is influenced by the impactor momentum alone
and $\mu$ = 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}$
(4)

From Eq. (1) it can be seen that $D_{crat}/D_{proj}$ is constant (for varying
$a$) only if $ga<<Y/\rho$, 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$
[28]. 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 cater size is obtained which is then used as input to compute the number
of expected impacts, following [28].

## 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.  _Computer Physics Communications_ 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.  _The Planetary Science Journal_ 3, 128 (2022).  URL https://iopscience.iop.org/article/10.3847/PSJ/ac67a7.
  * [7] Ormö, J. _et al._ Boulder exhumation and segregation by impacts on rubble-pile asteroids.  _Earth and Planetary Science Letters_ 594, 117713 (2022).
  * [8] Lundborg, N.  The strength-size relation of granite.  _International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts_ 4, 269–272 (1967).
  * [9] Collins, G. S., Melosh, H. J. & Ivanov, B. A.  Modeling damage and deformation in impact simulations.  _Meteoritics & Planetary Science_ 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.  _Nature Astronomy_ (2021).
  * [14] Lajeunesse, E., Monnier, J. B. & Homsy, G. M.  Granular slumping on a horizontal surface.  _Physics of 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.  _Planetary and Space Science_ 194, 105112 (2020).
  * [16] Lube, G., Huppert, H. E., Sparks, R. S. J. & Hallworth, M. A.  Axisymmetric collapses of granular columns.  _Journal of Fluid Mechanics_ 508, 175–199 (2004).
  * [17] Mitchell, J. K. _et al._ Mechanical properties of lunar soil: Density, porosity, cohesion and angle of internal friction.  _Lunar and Planetary Science Conference Proceedings_ 3, 3235 (1972).
  * [18] 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).
  * [19] Housen, K. R. & Holsapple, K. A.  Impact cratering on porous asteroids.  _Icarus_ 163, 102–119 (2003).
  * [20] Holsapple, K. A.  The Scaling of Impact Processes in Planetary Sciences.  _Annual Review of Earth and Planetary Sciences_ 21, 333–373 (1993).
  * [21] Raducan, S. D., Davison, T. M. & Collins, G. S.  Morphological Diversity of Impact Craters on Asteroid (16) Psyche: Insight From Numerical Models.  _Journal of Geophysical Research: Planets_ 125 (2020).  URL https://onlinelibrary.wiley.com/doi/10.1029/2020JE006466.
  * [22] Bottke, W. F. _et al._ Interpreting the Cratering Histories of Bennu, Ryugu, and Other Spacecraft-explored Asteroids.  _The Astronomical Journal_ 160, 14 (2020).
  * [23] Cho, Y. _et al._ Geologic History and Crater Morphology of Asteroid (162173) Ryugu Crater morphology and chronology of asteroid Ryugu.  _Journal of Geophysical Research: Planets_ n/a, e2020JE006572 (2021).
  * [24] Bierhaus, E. B. _et al._ Crater population on asteroid (101955) Bennu indicates impact armouring and a young surface.  _Nature Geoscience_ 1–7 (2022).  URL https://www.nature.com/articles/s41561-022-00914-5.  Publisher: Nature Publishing Group.
  * [25] Perry, M. E. _et al._ Low surface strength of the asteroid Bennu inferred from impact ejecta deposit.  _Nature Geoscience_ 1–6 (2022).  URL https://www.nature.com/articles/s41561-022-00937-y.  Publisher: Nature Publishing Group.
  * [26] Walsh, K. J. _et al._ Craters, boulders and regolith of (101955) Bennu indicative of an old and dynamic surface.  _Nature Geoscience_ 12, 242–246 (2019).
  * [27] Housen, K. R. & Holsapple, K. A.  Ejecta from impact craters.  _Icarus_ 211, 856–875 (2011).
  * [28] Jutzi, M.  The shape and structure of small asteroids as a result of sub-catastrophic collisions.  _Planetary and Space Science_ 177, 104695 (2019).
  * [29] Takaki, N. _et al._ Resurfacing processes constrained by crater distribution on Ryugu.  _Icarus_ 377 (2022).
  * [30] Tatsumi, E. & Sugita, S.  Cratering efficiency on coarse-grain targets: Implications for the dynamical evolution of asteroid 25143 Itokawa.  _Icarus_ 300, 227–248 (2018).
  * [31] Benz, W. & Asphaug, E.  Catastrophic Disruptions Revisited.  _Icarus_ 142, 5–20 (1999).
  * [32] 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.  _Meteoritics & Planetary Science_ 43, 1955–1977 (2008).
  * [33] Luther, R. _et al._ Momentum Enhancement during Kinetic Impacts in the Low-intermediate-strength Regime: Benchmarking and Validation of Impact Shock Physics Codes.  _The Planetary Science Journal_ 3, 227 (2022).  URL https://iopscience.iop.org/article/10.3847/PSJ/ac8b89.
  * [34] Raducan, S. D. _et al._ Influence of the projectile geometry on the momentum transfer from a kinetic impactor and implications for the DART mission.  _International Journal of Impact Engineering_ 162, 104147 (2022).  URL https://www.sciencedirect.com/science/article/pii/S0734743X21003225.
  * [35] Tillotson, H. J.  Metallic Equations of State for Hypervelocity Impact.  _General Atomic Report_ GA-3216, 141 (1962).
  * [36] 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).
  * [37] Jutzi, M. & Michel, P.  Hypervelocity impacts on asteroids and momentum transfer I. Numerical simulations using porous targets.  _Icarus_ 229, 247–253 (2014).
  * [38] 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).
  * [39] Schwartz, S., Richardson, D. & Michel, P.  An implementation of the soft-sphere discrete element method in a high-performance parallel gravity tree-code.  _Granular Matter_ 14, 363–380 (2012).
  * [40] 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).
  * [41] Schwartz, S. R. _et al._ Catastrophic disruptions as the origin of bilobate comets.  _Nature Astronomy_ 2, 379–382 (2018).
  * [42] 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).
  * [43] Schwartz, S. R., Yu, Y., Michel, P. & Jutzi, M.  Small-body deflection techniques using spacecraft: Techniques in simulating the fate of ejecta.  _Advances in Space Research_ 57, 1832–1846 (2016).
  * [44] Holsapple, K. A. & Schmidt, R. M.  Point source solutions and coupling parameters in cratering mechanics.  _Journal of Geophysical Research: Solid Earth_ 92, 6350–6376 (1987).
  * [45] Schmidt, R. M.  Meteor Crater: Energy of formation - Implications of centrifuge scaling.  In _Lunar and Planetary Science Conference_ , vol. 11, 2099–2128 (1980).
  * [46] Prieur, N. C. _et al._ The effect of target properties on transient crater scaling for simple craters.  _Journal of Geophysical Research: Planets_ 122, 1704–1726 (2017).

## 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
interdisciplinary programs.

## Author contribution

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 interest

The authors declare no competing interests.