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 craters1. 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 regime2 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 m2), 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 code3,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 in6 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 cs, allowing for a larger time step dt). This approach has been extensively tested and validated against laboratory experiments into homogeneous and heterogeneous targets7. 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 model8,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 another10. 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 (Y0 = 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/m3 (initial porosity ϕ0 = 50%), which is slightly higher than the bulk density of Ryugu (1190 kg/m3) measured by Hayabusa211. We include 50%, microporosity (as described in refs. 4, 12), which is consistent with recent analysis of Hayabusa2 samples returned from Ryugu13.

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 (Y0 = 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).

Fig. 1: Crater radius as a function of material properties.
figure 1

a Crater radius as a function of cohesion Y0 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 (Y0 = 0). The horizontal line shows the measured SCI crater size at the pre-impact level (7.25 ± 0.8 m), with its associated error2. Source data are provided as a Source Data file.

Full size image

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. 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 (⪅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 target1). 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 compaction1. The crater would form in the compaction cratering regime only if the compressive strength was of the same order as the cohesion (⪅1 Pa), while weak soils have compressive strength > 104 Pa18. 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 (Y0 = 0–0.5 Pa) and varying initial porosity (ϕ0 = 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 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 (Y0 = 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 point20 on Ryugu.

Fig. 2: End-to-end SPH simulation of the SCI impact.
figure 2

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 Ryugu2. 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;20). 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).

Full size image

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 observations2.

Fig. 3: Comparison of SCI impact simulations with varying target surface properties.
figure 3

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 Fig. 1. The crater radius in the case of the homogeneous target is 8.82 m.

Full size image

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\nu -2-\mu )}{3\mu }}\right.\\ \ \ {\left.+{\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}}\right\}}^{\frac{-\mu }{2+\mu }}$$
(1)

where the constants μ, ν, KR1 and KR2 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 K1 = 0.60, K2 = 0.25, μ = 0.44 and Y = Y0 = 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.

Fig. 4: Impact outcome as a function of impact size.
figure 4

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.

Full size image

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 Dcrat is the crater diameter and Dproj 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 Dcrat/Dproj 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 Dcrat/Dproj = 10 that has been suggested to be valid over a larger range of scales22. Important consequences of such a high cratering efficiency are high resurfacing rates and young crater retention ages. For example, a ratio Dcrat/Dproj of about 100 leads to an about 103 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 Dcrat/Dproj = 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,24,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 gravity23, 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 Dcrat/Dproj 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 Dcrat/Dproj of about 110 and a crater radius of about 11 m. Following the procedure used in27, 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 efficient28. 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 Bennu24. 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 basalt30, 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 materials8,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, Ydm is the limiting strength at high pressure and Y0 is the cohesion.

We use a constant cohesion, Y0, 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 in31). Our model uses a linear relation between cohesion Y0 and total strain ϵtot, and it is assumed that for ϵtot ≥ 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 Pmin ≥ − Y0.

To compute the coefficient of internal friction from the measured angle of response θ15,17, we follow the relation described in ref. 32:

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

The porosity was modelled using the P − α model4 with a crush curve defined by the parameters Ps, Pe, Pt, αt, n1, n2 and initial distension, α0. Pe is the pressure defining the transition between the elastic and plastic regimes, Ps the pressure at which all pores are compacted (distention α = 1) and the parameters Pe < Pt < Ps and 1 < αt < α0 indicate a transition between two regimes with different slopes (n1 and n2).

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 regime33. The relevant material parameters used in this study are summarised in Table 1.

Table 1 Material model parameters for impact simulations into Ryugu analogues
Full size table

Initial conditions

The target is modelled as a hemispherical domain with a radius of 15 m (Fig. 3) using approximately 107 SPH particles. We use a spherical projectile (Table 1) with a mass of 2 kg, impacting head-on. Recent studies33 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,34,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 × 10−4 m/s2.

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 cs, 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, cs: dt < resolution/cs ≃ 10−6 − 10−7 s. On the other hand, the crater formation time in the gravity regime can last up to few hundred seconds2.

Direct shock physics code modelling of late stage evolution

To model the late-stage evolution of the crater grow we follow the approach introduced in6 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 cs, allowing for a larger dt).

At this time ttransition 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 − ρ/ρ0), which also determines the magnitude of the sound speed. We use A ≈ 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 targets7.

The transition time ttransition depends on the scale of the impact. We use the following definition ttransition ≃ L/vsound × 10 where L is the length scale (chosen to be the size of the computational domain) and vsound the sound speed of the material. According to this definition, ttransition corresponds to roughly 10 sound crossing times. For the SCI-scale impact, we use ttransition = 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 × SCI-size; Fig. 4) as the simulations are faster to perform at larger scales. We use transition times of ttransition = 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 ( ± 1.1 m).

SPH-SSDEM framework for late stage evolution

An alternative approach to direct SPH simulations is to use a SPH-SSDEM framework36, 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 directions37,38.

The combination of the SPH and SSDEM methods has been used in previous studies to simulate small body impact processes e.g. 39,40. 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 constant40. 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 achievable41.

We use a shape-construction algorithm40 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).

Fig. 5: Comparison of two numerical approaches.
figure 5

Shown are the results of direct SPH (this work) and SPH-SSDEM simulations36. 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.

Full size image

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)20, 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.

Fig. 6: Boulder positions and impact geometry.
figure 6

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.

Full size image
Table 2 Position and properties of boulders in the simulations with heterogeneous targets
Full size table

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 ρ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 relationships42 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.443. 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}$$
(4)

From Eq. (1) it can be seen that Dcrat/Dproj is constant (for varying a) only if ga < < Y/ρ, i.e., in the strength regime. This means that using a constant Dcrat/Dproj = 1022 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 q27. 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, following27.