Modified Granular Impact Force Laws For The OSIRIS-REx Touchdown On The Surface Of Asteroid (101955) Bennu | Monthly Notices Of The Royal Astronomical Society

R-L Ballouz, K J Walsh, P Sánchez, K A Holsapple, P Michel, D J Scheeres, Y Zhang, D C Richardson, O S Barnouin, M C Nolan, E B Bierhaus, H C Connolly, Jr., S R Schwartz, O Çelik, M Baba, D S Lauretta, Modified granular impact force laws for the OSIRIS-REx touchdown on the surface of asteroid (101955) Bennu, Monthly Notices of the Royal Astronomical Society, Volume 507, Issue 4, November 2021, Pages 5087–5105, https://doi.org/10.1093/mnras/stab2365

The OSIRIS-REx mission collected a sample from the surface of the asteroid (101955) Bennu in 2020 October. Here, we study the impact of the OSIRIS-REx Touch-and-Go Sampling Acquisition Mechanism (TAGSAM) interacting with the surface of an asteroid in the framework of granular physics. Traditional approaches to estimating the penetration depth of a projectile into a granular medium include force laws and scaling relationships formulated from laboratory experiments in terrestrial-gravity conditions. However, it is unclear that these formulations extend to the OSIRIS-REx scenario of a 1300-kg spacecraft interacting with regolith in a microgravity environment. We studied the TAGSAM interaction with Bennu through numerical simulations using two collisional codes, pkdgrav and gdc-i . We validated their accuracy by reproducing the results of laboratory impact experiments in terrestrial gravity. We then performed TAGSAM penetration simulations varying the following geotechnical properties of the regolith: packing fraction (P), bulk density, inter-particle cohesion (σc), and angle of friction (ϕ). We find that the outcome of a spacecraft-regolith impact has a non-linear dependence on packing fraction. Closely packed regolith (P ≳ 0.6) can effectively resist the penetration of TAGSAM if ϕ ≳ 28° and/or σc ≳ 50 Pa. For loosely packed regolith (P ≲ 0.5), the penetration depth is governed by a drag force that scales with impact velocity to the 4/3 power, consistent with energy conservation. We discuss the importance of low-speed impact studies for predicting and interpreting spacecraft–surface interactions. We show that these low-energy events also provide a framework for interpreting the burial depths of large boulders in asteroidal regolith.

Granular material, in the form of regolith (i.e. broken-up rock particles), is ubiquitous in the uppermost layer of the surface of airless Solar system bodies, such as the Moon and asteroids. The terminology of sedimentary deposits can be used to describe the range of sizes of these particles. The terms pebble, cobble, and boulder describe individual particles of sizes ranging from 4 mm to 6.4 cm, 6.4 cm to 2.6 m, and > 2.6 m, respectively (Williams et al. 2006). For simplicity, when we refer to regolith, we mean the population of particles on a planetary surface that are the size of a pebble or smaller.

The main process for regolith formation on a large planetary surface is thought to be the comminution of large boulders by impact (e.g. Hörz et al. 1975). On the Moon, the impact excavation of boulders by the formation of large craters is balanced by their subsequent comminution by the constant bombardment of small meteoroids (e.g. Basilevsky et al. 2014; Costello et al. 2018). Spacecraft exploration of asteroids (e.g. Cheng et al. 2002; Fujiwara et al. 2006; Lauretta & DellaGiustina et al. 2019; Sugita et al. 2019) showed that their surfaces were also covered by a layer of regolith. Based on images returned of the asteroids (433) Eros and (243) Ida, regolith production was thought to be the result of the formation of large craters (≳1 km) in the gravity regime, where ejecta fragments are expected to be produced in large quantities and retained by the asteroid; in contrast, smaller craters (≲1 km) formed by meteoroids are thought to be less efficient sources for ejecta retention as they form in the strength regime (Geissler et al. 1996; Robinson et al. 2002). For the small (300 m diameter) and regolith-covered asteroid Itokawa, Barnouin-Jha et al. ( 2008) suggested that the surface regolith originated from its parent body. However, recent observations of ejecta formation and retention on small asteroids (101955) Bennu and (162173) Ryugu (∼0.5 and 0.9 km diameters, respectively) by small impactors has challenged this conventional view of ejecta retention by small craters (Arakawa et al. 2020; Perry et al. 2021).

Alternative mechanisms of regolith formation have been proposed, including the fragmentation of boulders through thermal cycling (e.g. Dombard et al. 2010; Delbo et al. 2014). The first evidence for thermal fragmentation on asteroids was the observation of debris aprons around boulders on Eros (Dombard et al. 2010). Through laboratory experiments and numerical simulations, thermal fragmentation has been shown to be as effective as impact comminution for small asteroids (El Mir et al. 2019). Images returned from the near-Earth asteroid Bennu by the Origins, Spectral Interpretation, Resource Identification, and Security–Regolith Explorer (OSIRIS-REx) mission (Lauretta et al. 2017, 2021) showed evidence of both thermal fatigue (Molaro et al. 2020a) and impact breakdown (Ballouz et al. 2020) of metre-scale boulders.

Spacecraft interactions with the regolith surface of a small Solar System body have resulted in various outcomes. These range from soft landings [the first of which was performed by NEAR-Shoemaker at asteroid (433) Eros (Dunham et al. 2002)], to metre-high bouncing [Hayabusa at asteroid (25143) Itokawa (Yano et al. 2006); MASCOT at (162173) Ryugu (Jaumann et al. 2019)], to spacecraft displacements of hundreds of metres [Philae at comet 67P/Churyumov-Gerasimenko (Biele et al. 2015)]. Understanding the underlying mechanics of a regolith surface in a low-gravity environment is crucial for understanding the response of a spacecraft to landing on an asteroid. Here, we study this interaction between spacecraft and regolith in the context of the Touch-And-Go (TAG) maneuver of the OSIRIS-REx spacecraft on Bennu.

NASA's OSIRIS-REx mission arrived at Bennu, a primitive near-Earth asteroid, in 2018 December (Lauretta & DellaGiustina et al. 2019). After a detailed observation campaign, the Touch-And-Go-Sample Acquisition Mechanism (TAGSAM; Bierhaus et al. 2018) was deployed to collect material from the asteroid's surface on 2020 October 20 (Lauretta et al. 2021). The TAGSAM consists of a sampler head with an articulated pogo-stick-like arm. The sampling strategy was to release nitrogen gas into the regolith, thereby entraining some material into the sampling device. However, the interaction of a sampling mechanism with a regolith surface in the low-gravity environment of a small body is difficult to predict, even if the physical conditions on the surface of an asteroid are well known. This is because the response of granular material to low-speed impacts is poorly understood.

KD2013 suggest that the |$\mu $| -dependence in the speed-squared inertial drag force could correspond to an added-mass effect, where the volume of the particle flow field grows in proportion to |$\mu $|⁠ . Furthermore, the friction component now depends on the density of the impactor. This suggests that frictional contacts are loaded by the motion of the projectile, such that the medium is stronger than that set by lithostatic pressure and Coulomb friction. Although Katsuragi & Durian ( 2007, 2013) developed a theoretical model to describe the complex response of granular materials to impact, and the dependence on material properties of the impactor and target, they did not characterize the influence of variable gravity.

Motivated to understand secondary re-impacts of material ejected after the formation of primary impact craters on asteroid surfaces, Nakamura et al. ( 2013) built upon the work of Katsuragi & Durian ( 2007) by conducting laboratory experiments and numerical simulations of 6-mm plastic projectiles impacting 50- and 420-|$\mu\rm m$| glass beads. They studied the deceleration profiles of the impactors and isolated the effects of the individual terms in equation ( 2), and they included an additional term, |${F_{{\rm{viscous}}}} = \,\,6\,\,{\rm{\pi \eta }}av$|⁠ , where |${\rm{\eta }}$| is the viscosity of the granular material. Of the force terms that characterize a granular bed's resistance to penetration (i.e. |${F_{{\rm{drag}}}}$|⁠ , |${F_{{\rm{pressure}}}}$|⁠ , and |${F_{{\rm{viscous}}}})$|⁠ , only the pressure term has a gravity dependence. By measuring the contributions of each term to the deceleration of the impactor, Nakamura et al. ( 2013) found that the gravity-dependent lithostatic-pressure term has a very marginal effect on the deceleration compared to the other effects. The lithostatic-pressure term was measured to have a magnitude 1/40 of the inertial-drag term. This suggested that the local gravity has a negligible influence on the outcome of a cratering impact. Altshuler et al. ( 2014) conducted impact experiments at gravities ranging from 0.4 to 1.2 times Earth gravity and found a similar lack of dependence of the penetration depth of an impactor on the local gravity. Although these studies suggest little to no gravity dependence on the outcome of a granular impact, the range of gravities were limited to values greater than ∼0.1 times Earth gravity (∼4 orders of magnitude greater than that found on the asteroid Bennu, for example). It is unclear whether this description of granular impacts is still valid at very low gravities.

To better match the gravity environment of small Solar system bodies, impact experiments in reduced gravity (as little as 10–4 times Earth gravity) have been performed using Atwood machines (Murdoch et al. 2017), parabolic flights with aircraft, and the International Space Station (Brisset et al. 2018). While these studies provided qualitative insight into the behaviour of small impactors in low gravity, they require specialized equipment and environments and thus are limited in their ability to probe the wide parameter space of impactor and target properties needed to obtain a rigorous understanding of the range of expected behaviours for a spacecraft touching down on an asteroid surface.

Here, our goal is to identify the governing physical parameters of an asteroid surface that determine its response to a landed spacecraft. Many important questions, especially that of the role of the near-zero gravity, cannot easily be answered with physical experiments. Physically simulating such a wide range of possible interaction scenarios on Earth is difficult and expensive. Therefore, numerical simulations are an ideal tool to explore possible scenarios. Maurel et al. ( 2018) and Thuillet et al. ( 2018) both used the N-body Soft-Sphere Discrete Element Method (SSDEM) code pkdgrav to explore the impact dynamics of DLR-CNES's MASCOT lander when deployed by JAXA's Hayabusa2 spacecraft on to the surface of asteroid Ryugu. They correctly predicted that MASCOT would bounce off the surface before eventually settling. The surface operation of NASA's OSIRIS-REx spacecraft was explored in Ballouz ( 2017) and Sánchez et al. ( 2013). They presented preliminary results on the range of possible behaviors of the TAGSAM-regolith interaction and found that varying inter-particle friction and cohesion can control penetration depth, but did not provide an empirically derived force model for this behaviour. Celik et al. (2019) explored the trajectory and impact dynamics of a small deployable camera in the Phobos environment in anticipation of JAXA's future Phobos sample return mission, MMX. They found that sufficiently high impact angles can lead to the ricochet. Here, we build upon the results first presented in Ballouz ( 2017) by analysing spacecraft-regolith simulations in the context of the force laws developed by the granular physics community.

We begin in Section 2 by introducing the code updates that were implemented to model the collisional interaction of a non-spherical body with spherical particles. In Section 3, we compare granular impact simulations with laboratory experiments that were performed to calibrate our code. In Section 4, we present direct numerical simulations of the interaction of the OSIRIS-REx spacecraft with the surface of Bennu, studying the effect of packing fraction, bulk density, cohesive strength, and angle of friction. In Section 5, we summarize our study, discuss our results in the context of natural low-energy phenomena on asteroids, and briefly discuss the applicability of the methods developed here for future space missions that will interact with or deploy a lander to a small Solar system body's surface.

In this section, we first introduce the SSDEM N-body code pkdgrav . Then, we describe code updates to pkdgrav that enabled the simulation of spacecraft interaction with regolith. Following this section, we describe a code validation exercise performed with pkdgrav and another SSDEM code, gdc-i (see Sánchez & Scheeres 2011 for details of the implementation).

pkdgrav is a parallel N-body gravity tree code adapted for particle collisions (Richardson et al. 2000, 2009, 2011). Collisions between particles in pkdgrav are treated using SSDEM (Cundall & Strack 1979). A pkdgrav -based implementation of the SSDEM was introduced in Schwartz et al. ( 2012) and has since been enhanced with improved functionality for modelling rolling friction and inter-particle cohesion (Zhang et al. 2017, 2018). The code uses a second-order leapfrog integrator, with accelerations due to gravity and contact forces recomputed at each step.

The soft-sphere implementation in pkdgrav uses a spring-dashpot contact law to model the collisional forces between particles. In this model, a spherical particle overlapping with a neighbouring particle feels a reaction force in the normal and tangential directions determined by spring constants (⁠|${k_n}$| and |${k_t}$|⁠ ), with optional damping and effects that impose static and/or rolling friction. The damping parameters (⁠|${C_n}$| and |${C_t}$|⁠ ) are related to the conventional normal and tangential coefficients of restitution used in hard-sphere implementations, |${ \in _n}$| and |${ \in _t}$|⁠ . The static, twisting, and rolling friction components are parametrized by dimensionless coefficients |${\mu _s}$|⁠ , |${\mu _t}$|⁠ , and |${\mu _r}$|⁠ , respectively. Additionally, Zhang et al. ( 2017) introduced a shape factor, β.

The coefficient of static friction, |${\mu _s}$|⁠ , determines the maximum amount of tangential force that can be supported by the contact point of two particles. In the event that the tangential force exceeds the maximum value given by the static friction coefficient, the particles slip past one another. The coefficient of rolling friction is used to determine the resulting torque due to two particles rolling against one another. The twisting friction coefficient is used to determine the torque arising due to a difference in the rotation rate along the normal vector of a particle-particle or particle-wall contact. The parameter β codifies the bulk influence of particle shape. In practical terms, non-zero coefficients of friction and β mimic some asphericity and roughness in the grains, increasing the bulk resistance to shear.

In pkdgrav , inertial walls are geometric primitives (such as spheres, rectangles, and cylinders) with a finite mass that can react to impulses by spherical particles (Richardson et al. 2011). These walls can be combined to form an assembly of walls that collectively reacts to impulses as a single rigid body. An inertial wall has the following key properties that allow it to behave as a rigid body: a mass, a centre of gravity, a spin vector, a moment of inertia tensor, and an orientation matrix. The exact geometry of a spacecraft or lander can then be reconstructed in the pkdgrav environment, and we can thus model its interaction with a granular surface.

TAGSAM-specific dynamics, including a constant force spring, were included in the code and are detailed in Ballouz ( 2017). Here, we summarize the details of the TAGSAM design and its implementation in pkdgrav , but the interested reader may refer to Bierhaus et al. ( 2018) for more insight into the actual mechanism. TAGSAM contacts the surface of Bennu at an initial approach speed of approximately 10 cm s−1. Once contact is established, nitrogen gas is discharged and flows through an annular nozzle around the lower part of the collector head. The gas mobilizes the regolith underneath, then around, TAGSAM (Bierhaus et al. 2021). The gas phase of sampling is not modelled here. Rather, we only model the pre-gas collisional interaction of TAGSAM with a regolith surface. We attempted to replicate its geometry as close as possible by combining six geometric primitives (three cylinders and three disks). The TAGSAM collector head is made out of an outer hollow cylinder with a diameter of ∼30 cm and an inner hollow cylinder with a diameter of ∼ 20 cm. Interior to the main cylindrical housing is a smaller deflector cone that tapers downwards. The TAGSAM is attached to the spacecraft by a pogo-stick–like robotic arm. The total mass of the spacecraft, including TAGSAM, is ∼1300 kg. The TAGSAM head can move with six degrees of freedom (DOF) as a single rigid body. We note that the real TAGSAM does not rotate about the axis that is co-linear with the arm. However, as modeled here, the TAGSAM can rotate about three axes. The inclusion of an additional degree of freedom has no bearing on our results as off-axis collisions with the surface are not studied. The spacecraft and TAGSAM dynamics are coupled by the arm, which contains a constant-force spring that is designed to begin compressing when the TAGSAM head feels a force that exceeds the spring force, which here we model as 67 N, based on the as-built flight unit.

To establish the validity of pkdgrav and gdc-i , we compare them against a series of low-speed impact experiments done in the laboratory under Earth’s gravity (9.81 m s−2). The experiments included impacts of spherical projectiles (diameter of 2.54 cm) made of steel, aluminium, and nylon onto two types of targets: dry sand and glass beads. Relevant properties of these targets are documented in Table 1. Dry sand is a common material for testing of impact phenomena, especially for hypervelocity cases at several kilometres per second (Holsapple 1994). The glass beads are commonly used by experimenters because of their uniformity and relative ease of handling (e.g. Schwartz et al. 2013). The target sample was placed in a hemispherical container approximately 0.5 m in diameter. It was brought to a stable and compact state by manually shaking the container until any settlement ceased. Then a projectile of a chosen type and size was dropped from various heights to impact the surface. Drop heights from zero to 3 m provide a range of impact velocities from zero to 7 m s−1. Impact experiments with each unique combination of projectile type, target type, and impact speed were repeated four to seven times. The spherical impactor penetration was measured as the height to the top of that impactor from a reference height after the impact. Table A1 documents all the outcomes of the impact experiments included in this study, and Table 2 summarizes these outcomes by reporting the mean penetration depth (dexp) and its standard deviation (σd).

Impact experiment target materials and their properties: ρ is the bulk density, r is the range in radii of individual particles, c is the sound speed of the material, K is the bulk elastic modulus, ν is the Poisson’s ratio, E is the Young's modulus, ɸ is the angle of friction, and P is the packing fraction.

Low-speed impact experiment outcomes for simulations and experiments. The impact speed (U), projectile type (Proj.), and target type (Target) for each case are shown. The penetration depths of an impactor in the simulations are shown in the columns Sim: pkdgrav and Sim: gdc-i . The mean of the impact experiment penetration depth is shown in the column dexp, and the standard deviation is shown in the column σd.

The impact simulations were set up by filling a cylinder with a radius of 10 cm and height of 9.2 cm (with a covered bottom) with ∼400 000 spherical particles. The spheres have a normal distribution of radii with a mean of 1 mm and a standard deviation of 0.1 mm. The distribution was cut off at 1 standard deviation to prevent the particles from crystallizing or packing too closely together. A crystalized packing can lead to artificial outcomes, such as a Newtown's cradle effect, owing to unnatural symmetry lines and planes. The spring constant of each particle was set such that the typical sound speed of an impact matched that listed in Table 1 across the particle. This is consistent with previous impact simulations that used a soft-sphere approach (Wada et al. 2006).

We conducted simulations of low-speed impact cratering onto the two granular materials used in the laboratory experiments: glass beads and dry sand. For impact experiments onto dry sand, we used coefficients of restitution parameters similar to that of gravel (Yu et al. 2014): ϵn and ϵt ∼ 0.55. Based on simulated uni-axial compression tests and matching a ɸ = 35°, we used |${\mu _s}$| = 1.0, |${\mu _r}$| = 1.0, |${\mu _t}$| = 1.3, and β = 0.7. The simulated dry sand had a Young's modulus ∼ 200 MPa and a Poisson ratio of 0.23. For impact experiments onto glass, we used coefficients of restitution parameters derived from previous studies (Schwartz et al. 2013): ϵn and ϵt ∼ 0.9. Based on simulated uni-axial compression tests and matching a ɸ = 22°, we used |${\mu _s}$| = 0.4, |${\mu _r}$| = 0.1, and β = 0.5. The simulated glass beads had a Young’s modulus ∼ 150 MPa and a Poisson ratio of 0.31. Details of the uni-axial compression and angle of friction tests were presented in Ballouz ( 2017). The impactor had a diameter of 2.54 cm, as in the laboratory experiments. Its density was varied between that of steel, aluminium, and nylon (7.6, 2.7, and 1.15 g cm−3, respectively). The outcomes of the impact simulations using pkdgrav are shown in Table 2.

We also performed code benchmarking activities with a separate implementation of SSDEM (Sánchez & Scheeres 2011), gdc-i. gdc-i has also been used to study the TAGSAM interaction with a simulated microgravity surface (Sánchez et al. 2013). The version of gdc-i used in this study employs a linear spring-dashpot contact law, similar to pkdgrav . Furthermore, the particles interact through static and dynamic friction, and a rolling friction term was also implemented to mimic the behavior of non-spherical particles. The simulations were set up similarly to those described in Sec 3.2. We used 68 000 particles with diameters of 4.5–5.5 mm for targets that represent glass beads and dry sand. Other material parameters were matched including angle of repose and the density of the particles. The particles were contained in a cylindrical container with a diameter of 30 cm and a height of ∼15 cm. The diameter of the cylinder containing the target material was >5 times the size of the projectile, which gives some assurance that the border effects are minimal (Seguin et al. 2008). All targets had a packing fraction of approximately 0.62. The impact experiment outcomes using this SSDEM implementation are summarized in Table 2.

We compare the outcomes of the simulations to the experiments, and to the experiment-derived empirical formulations from U2003 (equation 1) and KD2013 (equation 3). For both empirical formulations, we determine their predicted penetration depth for the target and projectile material properties and initial speed that were used in each corresponding experiment. For comparisons with the force law of KD2013, we numerically integrated equation ( 3) using the python package scipy and its odeint function, which solves a system of differential equations using the lsoda integrator (Virtanen et al. 2020), until the magnitude of the velocity reaches zero. In Figs 1(a)–(d), we show the cross-comparison of these five sources of penetration depth measurements or computations: two sets of simulation results (pkdgrav and gdc-i ), the outcome of the experiments, and two sets of computed results from the experimentally derived empirical formulas of U2003 and KD2013. We compare differences in outcome by showing the fractional difference in the penetration depth. We define the fractional difference in the measured penetration depth from a source A to source B as |${f_{{\rm{diff}}}} = | {{d_A} - {d_B}} |\,\,/{d_B}$|⁠ .

Cross-comparison of the five sources of penetration depth measurements or computations for each of the cases shown in Table 2. A value of zero on the vertical axis indicates perfect agreement. The vertical dashed black line separates cases where the targets were composed of glass beads (left side) from cases where the targets were composed of dry sand (right side). The comparison being made in each panel is indicated in in the top-left corner. (a) Experimental results compared to the two sets of simulation results: pkdgrav (open circles) and gdc-i (x's). (b) Experimental results compared to the two empirical formulations: KD2013 (open squares) and U2003 (closed red circles). (c) Simulation results compared to predictions from U2003. (d) Simulation results compared to predictions from KD2013.

The general trends in the simulation outcomes are similar to those in the laboratory experiments. As expected, higher impactor densities and higher impact speeds result in larger penetration depths (Table 2). In Fig. 1(a), we compare the simulation results to the experiments described in Section 3.1. Overall, we found good agreement, defined as fdiff ≲ 0.2, between the simulations and experiments for the cases of impacts into glass beads, but far more variability in simulation-to-experiment agreement for the dry sand cases. For impacts into glass beads (nine cases), pkdgrav and gdc-i simulations differed from experiments by the mean fdiff of 0.17 ± 0.12 and 0.15 ± 0.13, respectively, with the uncertainty indicating the standard deviation. For impacts into dry sand (seven cases), pkdgrav and gdc-i simulations differed from experiments by 0.26 ± 0.11 and 0.55 ± 0.22, respectively. In comparison, the standard deviation of the penetration depth for the experiments ranges from 0.02 to 0.14 of the mean for glass beads and 0.03 to 0.19 of the mean for dry sand. Therefore, some variation in the outcome is to be expected, but the simulations systematically overestimated the final penetration depth for the dry sand cases.

For the dry sand cases, a likely contributor to the large discrepancy between simulations and experiments is the difference in the radii of simulated particles compared to those in the experiments. The rougher nature of dry sand compared to glass beads may amplify differences between simulations and experiments, as particles are modelled as spheres. For well-packed systems, like the simulations and experiments achieved here, the contact area that an individual sand particle may experience can be significantly higher for irregularly shaped particles compared to more spherical particles. For smaller particle sizes with flat surfaces, the surface area for frictional contacts increases with decreasing particle size, providing more resistance to penetration. The influence of particle size on penetration dynamics has been previously studied (Lin & Wei 2012; Feng et al. 2019; Miyai et al. 2019). Here, the experiments used particles with radii rexp ∼ 0.35 mm, while pkdgrav used rpkd ∼ 1 mm particles, and gdc-i used rGDC ∼ 2.5 mm particles, such that rexp < rpkd < rGDC. For higher-speed cases (U > 1.4 m s−1; five of seven cases), the penetration depth outcomes systematically followed the same trend with experimental, pkdgrav , and gdc-i depths having dexp < dpkd < dGDC, respectively. Despite the relatively large fractional differences between simulations and experiments for dry sands, the simulations agree with the experiments for the single case of a nylon impactor (Case 14). For this particular case, it may be that the smaller pressure imparted by the lighter impactor does not place sufficient pressure on the sand particles to frictionally interlock, as was the case for denser impactors.

Panels (b) to (d) in Fig. 1 compare the simulation and experimental results to predictions from the empirical formulations of KD2013 (equation 3) and U2003 (equation 1). Overall, our simulations agree better with experimental results from the literature than they do with the experiments presented here, with variations in particle size as the likely cause. Fig. 1(b) shows that the empirical formulations also have difficulty in predicting the experimental outcomes for the dry sand case. For impacts into glass beads, the penetration depths predicted by KD2013 and U2003 differed from the experiments by a mean fdiff of 0.16 ± 0.13 and 0.15 ± 0.11, respectively. For impacts into dry sand, KD2013 and U2003 differed from experiments by 0.54 ± 0.26 and 0.40 ± 0.18, respectively. The disparity between the experimental results and the empirical formulations may also be attributed to differences in particle size. KD2013 and U2003 used a variety of granular materials for their experiments. The sizes of particles used by KD2013 in individual experiments ranged from r ∼ 0.15 mm (glass beads) to r = 3.5 mm on longest axis (popcorn). The particle size range used by U2003 spanned from r ∼ 0.15 mm (glass beads) to r = 4 mm on longest axis (rice). Although KD2013 and U2003 do not explicitly include a size dependency in their empirical formulations, their data points can diverge from the empirical curves by a factor of ∼ 2 (see fig. 5 of KD2013 and fig. 3 of U2003). This factor of ∼2 is similar to the maximum fractional differences we observe between the experiments, simulations, and predicted values from KD2013 and U2003. In support of this assessment of a particle size dependency, we show that the results of the simulations are well predicted by the empirical formulations of KD2013 (Fig. 1c) and U2003 (Fig. 1d). pkdgrav and gdc-i simulations differed from the empirical formulation of U2003 by a mean fractional difference of 0.21 ± 0.17 and 0.08 ± 0.06, respectively. The pkdgrav and gdc-i simulations differed from KD2013 by 0.17 ± 0.09 and 0.17 ± 0.12, respectively.

We take a closer look at the impact outcomes of the experiments and simulations into dry sand to assess whether the differences in the particle size can be accounted for within the existing framework of U2003. We choose the framework of U2003 as it provides a simple scaling formulation to compare with the combined results of simulations and experiments. This would be more difficult to accomplish with the force-balance framework of KD2013 as we were not able to record the detailed dynamical quantities of the impactor. Regardless, the well-established framework of U2003 also provides a common basis of comparison with KD2013, who showed that their experimental results were consistent with the penetration depth scaling of U2003 (equation 1).

(a) Comparison of the impacts into dry sand to the U2003 scaling law for penetration depth (equation 1). The gdc-i data (x’s) match the prediction from equation ( 1) (dashed line) best, followed by pkdgrav (open circles), and the experiments (filled triangles) have the worst match. (b) We introduce a particle-to-projectile-size ratio term to equation ( 1) and perform a linear least squares analysis to determine the dependency, α, of the penetration depth on this term, finding a best fit for α = 0.25. (c) We show penetration depth as a function of our new scaling law (equation 10) that introduces a dependency on the average particle size. The data are well modelled by equation ( 10), which is shown by the black dashed curve, giving R2 = 0.98.

We compare the values of penetration depths from numerical and laboratory results to equation ( 10) in Fig. 2(c). As expected, we find that all penetration depths, whether derived from simulations or experiments, collapse into a single trend line given by equation ( 10). This formulation of the penetration depth is likely valid for r/a < 1. When r/a > 1, the impactor is smaller than the mean target particle, leading to a bounce rather than penetration for low-speed impacts, when target and projectile have similar densities. We have demonstrated that the discrepancy between our simulations and experimental results can be explained by the difference in the size of the grains that make up the target. By understanding the factors that drive the systematic differences between simulation and experiments, we can more confidently move forward to use our validated codes to simulate conditions that are difficult to replicate in the lab. Furthermore, our results show that simulations of impacts into granular material can be sensitive to target grain size, producing differences up to a factor of 2 in the results presented here.

In this section, we describe simulations of the TAGSAM head penetrating a bed of granular material. As the properties of Bennu’s regolith surface were largely unknown prior to arrival, our strategy was to explore Bennu’s interaction with a regolith surface that had variable surface friction, cohesive properties, and packing fraction. Simulations such as these are valuable because they provide insight into the range of possible outcomes of a sample attempt by the spacecraft. Furthermore, data returned from the touchdown event itself will provide constraints on the geotechnical properties of the surface and sub-surface once compared to simulation results. Observations during the global surveys of Bennu provided measurements of some surface properties such as particle size, particle size frequency distribution, and mineralogy (e.g. Lauretta & DellaGiustina et al. 2019 and references therein). Given that these characteristics are fairly well constrained, the sampling attempt can also act as a useful and rare impact experiment onto particles in a microgravity setting, which is difficult to achieve experimentally on Earth (Brisset et al. 2018; Bierhaus et al. 2021; Murdoch et al. 2021). Based on the spacecraft’s reaction, and with the added context of the simulations described here, we may be able to better understand the dynamics of granular material in a microgravity environment.

We performed a series of simulations to separately explore the effects of (i) friction, (ii) the packing fraction of the granular bed, and (iii) cohesion. These characteristics of a granular surface are expected to largely govern its interaction with a spacecraft. We studied (i) and (ii) using the pkdgrav implementation described in Section 2. In the case of (iii), we explored the effect of inter-particle cohesion with gdc-i . The overall simulation strategy for modelling the complex 6+ DOF of the TAGSAM, the spacecraft arm, and the constant force spring have been described in Ballouz ( 2017) for pkdgrav and Sánchez et al. ( 2013) for gdc-i . We point the interested reader to those publications for further details, but elaborate briefly on the details here when appropriate.

We used pkdgrav to study the sensitivity of the TAGSAM outcome to ɸ. We set up simulations by filling a cylindrical container with a height of 60 cm and a radius of 60 cm with ∼ 150 000 particles and allowing them to settle in a gravity environment similar to Bennu's (∼ 7 |$\times $| 10–5 m s−2; Daly et al. 2020). Then, the actual touchdown simulation was conducted by placing the TAGSAM just above the surface of the granular bed and allowing it to move downward at a speed of 10 cm s–1 (Bierhaus et al. 2018). The total simulation time was ∼ 5 s. Each particle had a density of ∼ 2.8 g cm−3, which is typical of denser carbonaceous chondrites (Macke et al. 2011), the meteorites that are spectrally linked to carbonaceous asteroids such as Bennu. The regolith was composed of a continuous power-law size frequency distribution of spherical particles, with a differential power-law exponent of –2.5. Although the power law distribution of particles on the surface of Bennu is steeper (DellaGiustina et al. 2019; Burke et al. 2021), we selected this exponent because it balances the reality of the Bennu surface with computational limitations; steeper power laws require a computationally unfeasible number of small particles to be simulated within a space large enough to correctly capture the collisional dynamics of a spacecraft touching down on a regolith bed. We generated a continuous size distribution of particles with a minimum radius, rmin, of 0.5 cm, and a maximum radius, rmax, of 1.5 cm. Fig. 3 presents snapshots from one of the pkdgrav simulations, showing TAGSAM penetrating a bed of particles. The regolith bed had a packing fraction of 0.61, which is similar to that of random close-packing for monodisperse particles (∼ 0.63), and smaller than that of the hexagonal close-packing (∼0.74).

pkdgrav simulations of a model TAGSAM (visualized with slight transparency) penetrating a bed of approximately 150 000 particles (seen in cross-section) with an initial downward vertical speed of 10 cm s−1 in the microgravity environment of Bennu. The particles are coloured to represent distinct 6-cm layers. The vertical extent of the regolith bed is 60 cm, and the horizontal extent is 120 cm. The orange-red splotch is the shadow of the TAGSAM.

We explored the effect of variations in the coefficient of static, tangential, and rolling friction of the particles, and variations in the shape factor, β. However, we kept the normal and tangential coefficients of restitution constant at 0.55. Table 3 summarizes the friction and shape parameters used, their corresponding ɸ (measured through auxiliary angle of repose simulations), and the resulting peak force, Fpeak, felt by the spacecraft after touching down on the granular surface. The container size was sufficiently larger than the impactor such that there should be no boundary effects (Seguin et al. 2008). Furthermore, we performed a few test cases without a container, for friction end-members, and found that the peak forces were the same as cases where the container was included.

Summary of outcomes for pkdgrav simulations of the TAGSAM penetrating a regolith bed with a packing fraction of 0.61. Sim #, simulation number. μs is the coefficient of static friction, μr is the coefficient of rolling friction, μt is the coefficient of twisting friction, β is a shape factor, ɸ is the angle of friction of the granular bed, and Fpeak is the maximum force felt by the spacecraft.

The constant-force spring on the TAGSAM arm introduces a non-trivial dimension to the collisional dynamics of the systems. Therefore, we analyse one of the cases that demonstrates the dynamics of the constant force spring system in detail, in order to demonstrate the complex interaction between TAGSAM, the surface, and the spring. Fig. 4 illustrates the interaction between spacecraft and regolith surface for Sim 19. Fig. 4 shows how the force on the spacecraft is tracked (filled black circles, left vertical axis), and its effect on the compression of the constant force spring between the TAGSAM and the spacecraft (red dashed curve, right vertical axis). For this case, the regolith bed has ɸ ∼ 30°. The numbered regions in Fig. 4 (blue shading) indicate key events in the dynamical interaction between the spacecraft, the spring, and the TAGSAM as follows:

At t ∼ 0.5 s, the TAGSAM makes contact with the ground and experiences a sharp impulse of ∼140 N. The spring immediately engages as the force threshold of 67 N is exceeded.

The spring engagement causes a force balance between the downward-moving TAGSAM and the ground. The force by the ground on the TAGSAM is roughly at the spring threshold (67 N). The TAGSAM has a slight downward speed of 2 cm s−1. Due to the force balance, this remains roughly constant.

The spring compression reaches a maximum (t ∼ 1.7 s).

Once the spacecraft stops moving downwards, the spring begins to expand. This causes the spring force to drop to the 47–49 N range. Because the TAGSAM is forced downward by the spring, it continues to interact with the ground in this force range. The ground pushes back at the TAGSAM with an equal force. The TAGSAM remains in force balance. At the same time, the spring expands, as the spacecraft begins to move upwards. This continues until the spring reaches the 0-cm hard stop (t ∼ 3.3 s).

The spring is no longer engaged and the ground force is below the spring threshold. The spacecraft and TAGSAM continue moving upwards as a single rigid body. The spring force remains at zero until the end of the 5-s simulation.

Time evolution of the ground force on the TAGSAM (F, black filled circles, left vertical axis) and the spring compression (red dashed line, right vertical axis) for Sim 19 (Table 3). The spring engages after the TAGSAM feels a ground force that exceeds 67 N. The force by the spring depends on whether it is undergoing compression or expansion. This behavior is described by equations 3.11 and 3.12 in Ballouz ( 2017). The roman numerals indicate key events; see the text for details.

Fig. 4 illustrates how the spring dictates the dynamical interaction between the TAGSAM and the spacecraft, as well as modulating the total force that the TAGSAM feels during the sampling attempt. In the following section, we discuss how we model this interaction for different surface mechanical properties and how these properties can be inferred from the spacecraft's response to the sampling event (Lauretta et al. 2021). Details of how the spring compression is related to TAG dynamics were given in Ballouz ( 2017). Here, we focus on the dependence of the force response on the geotechnical properties of the regolith bed.

In our simulations, we find that ɸ can have a strong influence on a material’s response to an impact by a relatively massive spacecraft (Fig. 5). In Fig. 5(a), we show the penetration depth of the TAGSAM 1 s after contact, as this is exactly the time when the N2 gas from the sample head is blown into the regolith bed to mobilize material for sample collection (Bierhaus et al. 2018). For ɸ ≲ 28°, the granular bed offers almost no resistance, and the spacecraft-TAGSAM assembly is able to sink into the granular bed without triggering the constant force spring, or triggering it for only a fraction of a second. For ɸ ≳ 28°, the granular bed strongly impedes the downward motion of the TAGSAM, causing the constant force spring to trigger, as seen in Fig. 4. This limits the penetration depth to values < 4 cm. We performed a set of auxiliary simulations that tested the outcome of a touchdown without a constant force spring and found that for high ɸ cases (ɸ ≳ 28°), the spacecraft could bounce off the surface, reminiscent of the Philae lander on comet 67P/C-G (Biele et al. 2015).

Outcome of simulations of TAGSAM touching down on a regolith surface with variable angle of friction, ɸ. Simulations where particles have μS = 0.5, 1.0, 1.31 are represented by open black circles, open red squares, and blue x’s, respectively. (a) The penetration depth of TAGSAM after 1 s of contact with the surface depends on ɸ. (b) For these simulations, the peak force Fpeak felt by the spacecraft grows as ɸ to the power of 3/2.

As shown in Table 3, different combinations of the friction properties can lead to similar bulk ɸ (values of ɸ within 1° can be recreated using different combinations of μs, μr μt, and β). Despite this degeneracy in model parameters, we find that the final bulk properties of the regolith bed would lead to similar outcomes for penetration, as shown in Fig. 5(a). Fig. 5(a) shows that, for the regolith bed properties used in these simulations, ɸ is a good predictor of penetration depth and is independent of the specifics of the modelling parameters. The general shape of penetration depth appears similar to a negative sigmoid function. The saturation for low penetration depths (high ɸ) is due to the action of the constant force spring modulating the maximum force felt by the spacecraft. Sigmoid functions are widely studied in the field of machine learning as an activation function for the firing of artificial neurons. In an analogous sense, the sigmoid shape of the penetration depth emerges owing to the activation of the constant force spring. The saturation at high penetration depths (low ɸ) is likely due to a decreasing influence of friction in resisting penetration. Overall, the emergent qualities from the simulations suggest that increasing friction properties of the particles can control a transition from an unjammed state to a jammed state for the regolith bed. Fig. 5(b) shows that Fpeak ∝ ɸ3/2 once the constant force spring threshold is exceeded. For low ɸ (ɸ < 20°), it appears that Fpeak may grow even faster; however, it is difficult to say for certain because the number of cases explored for that region is small. The growth rate of Fpeak with ɸ is likely specific to the TAGSAM interaction as the constant force spring modulates the maximum force felt by the spacecraft.

To better understand why the penetration depth depends so sensitively on the friction properties of the particles, we can draw comparisons with studies into shear-thickening fluids (Mukhopadhyay et al. 2018). Mukhopadhyay et al. ( 2018) conducted studies on granular material (cornstarch) suspended in fluids. Such suspensions have the surprising property of resisting heavy loads when the impact speed of the projectile is sufficiently high. As a projectile impacts a fluid suspension, a mechanical wave travels through the suspension that causes the entire medium to dilate leading to the interlocking of nearby particles. It is possible that a granular material in low gravity may behave like a colloid in Earth gravity. Once the sampler makes contact with the surface, the material underneath the sampler shear-thickens. As the material underneath the sampler expands, the highly frictional particles form more particle–particle contacts. This creates a transition where the particles move from fluid-like flow to a solid-like response. In doing so, the strength of the surface increases, decreasing the penetration depth of the projectile. We provide examples of this shear-thickening behaviour in Fig. 6, which shows snapshots of the end of two simulations. The top row shows a case with high ɸ (Sim 19, ɸ = 30.6°), and the bottom row shows a case with low ɸ (Sim 1, ɸ = 17.8°). Figs 6(a) and (b) show the TAGSAM penetration depth at the end of the simulations, tend ∼ 4.5 s after initial contact. Figs 6(c) and (d) show the change in the column density of particles, for a 10-cm-thick slice through the centre of the simulation space, comparing the initial distribution of particles to that at tend. Blue colours correspond to underdensities with respect to the initial density distribution of particles. Redder colours correspond to overdensities. The large swaths of dark blue indicate regions where particles have been excavated. Fig. 6(c) shows the high ɸ case, where the excavation profile (blue) resembles a conventional crater, except for the region directly below TAGSAM. In Fig. 6(c), the region directly below TAGSAM appears to have dilated, suggesting that the particles interlocked and transitioned to a jammed state (Liu & Nagel 1988). Fig. 6(d) shows that the particles below the TAGSAM never develop into a jammed state. Rather, we see regions of local overdensities (red) as the particles are able to easily compress, filling voids and leaving behind new ones. The compressibility of the material allows the spacecraft to easily penetrate into the surface. The case shown in Fig. 6d would be analogous to a well-packed underdense smooth material, like hollow plastic balls in a children's ball pit, which allow denser objects to easily sink to the bottom. In this case, the TAGSAM represents a relatively small surface area where the entire spacecraft mass (∼ 1300 kg) interacts with the granular surface, resulting in a high contact pressure. Taking the spacecraft mass and dividing by the TAGSAM volume gives an impactor density at TAG, δTAG ∼ 54 g cm−3. This calculation is only valid if the spring is not engaged. Once the spring is engaged, the TAGSAM-regolith interaction is governed by the much smaller TAGSAM mass alone. However, as we will show in Section 5.1, the impactor density plays a small role in dictating the interaction in a microgravity environment.

Snapshots of simulations showing the difference in the penetration depth and granular dynamics for a high-friction case (top panels) and a low-friction case (bottom panels). (a) The final position of TAGSAM for a high-friction case (Sim 19). The particles are coloured to represent distinct 6-cm layers. (b) Same as (a) but for a low-friction case (Sim 1). (c) A map of the change in the column density (Δn) of the granular bed in (a), comparing the final to initial number density of particles through a 10-cm-thick slice centred in the middle of the granular bed. Bluer colours correspond to underdensities with respect to the initial density distribution of particles. Redder colours correspond to overdensities. The excavation profile resembles a conventional crater. The region directly below TAGSAM dilates as the particles interlock. (d) Same as (c) but for the low-friction case shown in (b). The excavation profile is thin and deep, resembling cratering into highly porous materials.

We further analyse the differences in penetration outcomes between low and high ɸ cases by considering the force chain networks that the granular bed develops as the TAGSAM touches down on its surface. A force chain consists of a linear string of rigid particles in point contact. A force chain network emerges when a granular material experiences an external load leading to the formation of interconnected force chains within the system, creating a jammed state where the material behaves like a solid (Cates et al. 1998). In this state, the granular system can support large applied loads in the direction of the original jamming force. However, individual chains may collapse if forces are applied in a different direction, and the system may subsequently lose its rigidity. Clark et al. ( 2012) showed that force chains exhibit intermittent behaviour, fluctuating rapidly in timescales approximately equal to the sound speed divided by the size of the impactor. Here, we only show snapshots of these force networks to illustrate how the networks characteristics can vary as a function of the geotechnical properties of the target granular material. In Fig. 7, we show the force chain networks that develop 0.2 s after the TAGSAM contacts the surface for a low and high ɸ. This time in the simulation co-incides with the peak force felt by the TAGSAM. As such, we define this time as t = tpeak. For regolith with low ɸ, the force chains are relatively shallow and weak (Fig. 7a). In Fig. 7(a), we see that the force chain network develops directly underneath the area of contact and only extends a few particle lengths below the surface. For regolith with high ɸ, the force chains are relatively deep and strong (Fig. 7b). In Fig. 7(b), we see that the force chain network extends from the contact area to deep within the sub-surface. The force chain network extends vertically down in the same direction as the applied load from TAGSAM, and some of its branches subsequently arc towards a perpendicular direction from the applied load. This allows a significant amount of mass in the subsurface to participate in resisting the penetration of the spacecraft.

Snapshot of simulations showing the force-chain network for two cases. The colour represents the normal force (N) that each particle experiences at t = tpeak. (a) For regolith with low ɸ, the force chains are relatively shallow and weak. (b) For regolith with high ɸ, the force chains are relatively deep and strong.

Simulations and experiments of granular systems show that for collections of particles, a critical packing fraction exists for which there are enough contacts per particle, known as the coordination number Nc, such that material can transition to a jammed state. Above the critical coordination number (Nc,crit), Nc and stress can increase rapidly as a function of P (O'Hern et al. 2002, Majmudar et al. 2007). For frictionless 3D spheres Nc,crit ∼ 6 (O'Hern et al. 2002). For the simulations described in this section, the median Nc = 3 < Nc,crit for frictionless spheres, which is consistent with our observation that TAGSAM can easily penetrate granular beds with low ɸ. However, it appears that for the high-friction cases, Nc,crit ∼ 3. In Figs 8 and 9, we illustrate the change in the volume and strength of the force chain networks in Sim 8 to 15, which span the end-members of friction cases that we studied. Figs 8 and 9 show the number and strength of particle–particle contacts in the normal and tangential direction, respectively, at t = tpeak. As ɸ increases, the magnitude of the force contacts increases, and the networks begin to branch perpendicular to the direction of loading, increasing the volume of particles that participate in resisting penetration. It is clear from Figs 8 and 9 that particles with sufficient friction contribute to the transition to a jammed state. Evidently, motion-loading of frictional contacts can enable a regolith bed to resist the penetration of a spacecraft, even in low gravity. This is consistent with laboratory experimental results that explore the role of friction in granular penetration ( KD2013) and counters conventional understanding that friction effects are negligible in the microgravity environment (Zacny et al. 2018).

Particles in the x–z plane that have relatively large contact forces (>0.5 N) at t = tpeak. The subpanels show the contact network for different values of the friction angle, ɸ, shown in the top-left corner. Each scatter-point represents the location of a contact, and its colour represents the magnitude of the normal force, Fn.

Same as Fig. 8, but for the tangential force, Ft, for each contact.

For regolith on the surface of an asteroid, the jamming transition is controlled by packing fraction and applied load (Liu & Nagel 1988). As the applied load by TAGSAM is known, we continue on to analysing the influence of the former property. The packing fraction can influence the manner in which the regolith can resist the load of a spacecraft by altering its bulk density, the average number of particle-to-particle contacts, and potentially ɸ (Metcalf 1966).

In this subsection, we investigate the effect of packing fraction, P, on the outcome of TAGSAM touching down on regolith. We implemented a method to generate packings with spherical particles for a user-specified P. Following Ringl et al. (2012), we build porous targets with a specified P through an algorithm that implements the following steps:

Set a particle at an arbitrary position in a specified volume that we wish to fill.

Calculate a temporary P for each particle by finding neighbours within a sphere of radius 3r. For particles on the faces, edges, and corners of the cubic volume, increase P by 1/2, 3/4, and 7/8, respectively.

Determine the particle with the smallest local packing fraction.

Attach a particle to it in a random direction.

Repeat from step 2 until the desired global packing fraction is reached.

This formulation enables the systematic study of the influence of P on low-speed impact outcomes. Our algorithm distributes the particles homogeneously, yielding a constant packing fraction per unit cell, within the specified volume (Fig. 10).

Initial conditions with variable user-defined packing fractions generated by the algorithm described in Section 4.2. The initial conditions are analysed using a 3D Voronoi tessellation to ensure spatial homogeneity of the packing fraction. The colors denote the mean packing fraction through a 5 cm × 5 cm cell in the x–z plane. For each case the median packing fraction, Pmed, of all the cells in the x–z plane is shown in the top left corner of the panel. (a) The spatial distribution of packing fraction for Pmed = 0.15. (b)–(d) Same as (a) but for Pmed = 0.2, 0.3, and 0.4, respectively.

To study the effect of P, we constructed cubic porous targets with lateral dimensions of 160 cm, vertical dimension of 100 cm, made up of particles with r = 1 cm, and packing fractions that ranged from 0.15 to 0.40. In Fig. 10, we verify the homogeneity of the packing by performing a 3D Voronoi tessellation analysis on the granular bed using the voro ++ package (Rycroft 2009). We chose to not gravitationally settle the particles as we solely want to study the behavior of non-gravity forces in the interaction between TAGSAM and the surface of a low-gravity planetary body for specific packings; thus, we set the net acceleration due to gravity in this subset of simulations to zero (in all dimensions). As shown by Nakamura et al. ( 2013), even under terrestrial gravity conditions, the drag force term shown in equation ( 3) (proportional to the square of the velocity) can dominate over the hydrostatic gravity term (proportional to the product of bulk density, gravity, and the depth of the impactor). One important caveat to this choice is that gravity, however little, may still influence the manner in which particles relax on the surface of an asteroid. This relaxation time is many orders of magnitude greater than the TAGSAM–regolith interaction; nevertheless, the manner in which particles settle could influence how their inherent material properties like friction can control the low-speed intrusion of the TAGSAM. We studied the interaction for two ɸ end-members, ɸ = 18.8 and 35.2, with the same static, rolling, and twisting friction parameters as described by Sim 8 and Sim 13 in Table 3.

Further, as we would like to characterize the influence of packing fraction on the governing velocity-dependent force, we also vary the impact speed U, from 0.1 to 0.5 m s−1, and the particle density, ρg = 2.2 or 3 g cm−3. We summarize the simulations completed for the packing fraction analysis and the resulting peak force, Fpeak, in Table 4.

Summary of outcomes for simulations of the TAGSAM with variable impact speed penetrating a regolith bed with variable packing fraction and particle density. P is the packing fraction of the bed, ρg is the density of individual particles, and Fpeak is the peak force felt by the TAGSAM.

To analyse the dependence of the force on packing fraction, we compare Fpeak to a velocity-squared dependent function as shown in the second term of equation ( 3): |${F_{{\rm{peak}}}} \propto K\mu \rho A{U^2}$|⁠ , where K is a fitting constant related to the shape of the impactor (Nakamura et al. 2013), |$\rho \,\, = \,\,P{\rho _{\rm{g}}}$| is the regolith bulk density (accounting for porosity), A is the TAGSAM contact area, U is the impact speed, and = tan(ɸ) is the coefficient of friction. In the limited range of ɸ values studied here, we found that in these low P regimes, the peak force is not sensitive to ɸ. In Fig. 11(a), we plot Fpeak, for every simulation described in Table 4, against velocity-squared force expression, finding a best-fitting value of K = 8 and an intercept of 15.7 N (see horizontal axis). The black dashed line is a 1:1 correlation that demonstrates how the model (horizontal axis) describes the data. Although the model fits the data fairly well (correlation coefficient R2 = 0.92), the intercept value does not have a valid physical significance. This would imply that a residual force would exist when U = 0. Furthermore, values of low Fpeak ≲ 25 N appear to diverge from the model. Therefore, although the U2 model provides a passable fit to the simulation data, there are some physical inconsistences to the model.

Peak forces, Fpeak, as a function of empirically derived models that describe the behaviour of TAGSAM penetrating loosely packed particles without the influence of gravity. The open circles are simulation results described in Table 4. (a) Fpeak is compared to the best-fit force model that has a velocity-squared term. Although the R2 value is relatively high, the model does not adequately capture the behavior for small values of Fpeak and requires an unphysically large value for the y-intercept. (b) The data are modelled with a velocity taken to the 4/3 power, following U2003. The data are better described by this model than that shown in (a).

The discrepancy between a velocity-squared term as suggested by KD2013 and our findings may also potentially be resolved by including a viscous term, which has a force term directly proportional to velocity, as suggested by Nakamura et al. ( 2013). However, this term was found to be a fraction of the velocity-squared term in that study. This result implies that the leading coefficient in the second term in equation ( 14) has a unit equivalent to dim[U2/3]. There are limitations to the applicability of this equation. For example, as shown in Section 4.1, well-packed systems exhibit a non-linear dependence on ɸ. Furthermore, energy conservation may not be applicable for higher-speed impacts where a non-negligible fraction of the energy may go into heating material or propagate to far distances as acoustic waves.

Our finding of a sensitivity of the penetration dynamics with packing fraction is consistent with previous studies (e.g. Umbanhowar & Goldman 2010; Royer et al. 2011), who varied ɸ over a narrow range of 0.51–0.62. Umbanhowar & Goldman ( 2010) identified a critical packing state, ɸc ∼ 0.59, where sheared grains neither dilate nor consolidate. Royer et al. ( 2011) find an impacting object can generate a pronounced compaction front in loose regolith (ɸ < ɸc) or dilate compacted regolith (ɸ > ɸc). Royer et al. ( 2011) noted that for loose beds, in particular, the behaviour is similar to that observed when a plow is pushed into loose soil or snow. Finally, inter-particle cohesion may play a significant role in resisting penetration if the particles can create sufficiently strong bonds (Sánchez & Scheeres 2014). We explore the role of cohesion in the next section.

This parameter is the strength provided by cohesive particles to cement larger ones together. This formulation would allow us to remove the smallest particles in the distribution and replace them with a net particle–particle tensile strength between the larger ones. This tensile strength is directly related to the material properties of the binding particles. With this implementation, the number of particles to be simulated remains manageable with the available DEM codes, while still including the effects of large numbers of fine particles. Sánchez et al. ( 2013) presented TAGSAM simulations performed with gdc-i on cohesive regolith with ɸ = 25°, varying the spring constant and gravity. Here, we present key updates to modelling interactions with cohesive regolith using gdc-i . For the gdc-i simulations, a bowl-shaped container, with diameter = 90 cm, was filled with 210 000 spherical particles with diameters between 0.8 and 1.2 cm, ɸ = 35°, P = 0.64, ρg = 2.5 g cm−3, and local gravity = 1 |$\mu$| g. The particles are settled using the procedure described in Sánchez & Scheeres ( 2011). The TAGSAM arm has a length of 1 m and can contract down to 70 cm. The spring in the arm engages when the force on the TAGSAM head is ≧67 N. The restoring force of the spring has the form 67 + 15(1 – 1.2c), where c is the compression of the arm in metres. Before the spring engages and when the arm has been fully contracted, the system moves as a solid object, but the TAGSAM is still free to rotate in two dimensions. The spacecraft has the inertia tensor of a cube of volume 1 m3 with a mass of 1292 kg. There is no hinge between the spacecraft and the arm, and these two entities are always aligned. The TAGSAM head has the inertia tensor of a solid cylinder with a mass of 8 kg. The system has 10 total DOF. The spacecraft has six DOF, the TAGSAM has three rotational DOF, and one additional DOF from the compressible arm.

We explored the effect of variable σc ∈ [25 Pa, 300 Pa] to assess the range of outcomes for the interaction and whether we would be able to discriminate between cohesive and strengthless regolith from spacecraft telemetry. Again, σc is the particle–particle cohesive strength and not the cohesive strength of the granular bed. For a regolith bed with ɸ = 35°, gdc-i simulations of the stability of rubble-pile asteroids against rotational disruption showed that the bulk cohesive strength is a factor of 0.2 of σc (Sánchez & Scheeres 2016, 2018). Fig. 12 shows snapshots of two simulations with gdc-i . Fig. 12(a) shows that a TAGSAM impact into a regolith with σc = 25 Pa leads to deep penetration, with almost no resistance from the regolith. Fig. 12(b) shows that increasing the cohesive strength slightly to σc = 125 Pa leads to a much different outcome, with the surface having sufficient strength to resist penetration.

Using gdc-i , we simulated TAGSAM impacts into regolith with variable interparticle cohesive strength, σc. These snapshots show the simulation state 2.7 s after contact. (a) Regolith with σc = 25 Pa provides little resistance to TAGSAM. (b) Regolith with σc = 125 Pa can effectively limit TAGSAM penetration.

In Fig. 13, we demonstrate how the cohesive strength of surface regolith may be measured using a TAGSAM-like system. As σc increases, the duration of spring compression increases. For σc = 50, 150, and 300 Pa, the total magnitude of spring compression was 1, 2.5, and 8.2 cm, respectively, and the total duration was 0.8, 2.0, and 3.6 s, respectively. The increase in cohesive strength of the regolith decreases its compressibility, leading to a higher and longer compressions of the spring. This ability of the spring to measure the effective strength of the surface was also demonstrated for the high-P, high-ɸ cases in Ballouz ( 2017). Here, the spring may be able to measure cohesive strengths σc > 25 Pa. For smaller values, the spacecraft would experience minimal resistance (Fig. 12a); however, further work is required to determine whether the relatively small ratio of TAGSAM diameter to bowl diameter influences the dynamics in the low-cohesion cases.

Force profiles from the gdc-i simulations for TAGSAM touching down on cohesive regolith with variable cohesive strength, σc. The force felt by the body, centre of mass (CM), and TAGSAM head are represented by red plus signs (which blur together into solid lines), the green dashed line, and the blue dashed line, respectively. (a)–(c) σc = 50 Pa, 150 Pa, and 300 Pa, respectively. Larger values of σc lead to larger arm compressions.

The effective compressive stress for a case where the constant force spring is just able to trigger is simply the force threshold divided by the surface area of TAGSAM: ∼500 Pa. This is approximately an order of magnitude greater than the minimum cohesive strength required for the spring to trigger, when also considering contributions from the frictional forces. This is consistent with observations of porous materials that have compressive strengths larger by a factor of 5–10 than their tensile strengths (e.g. Nakamura 2017).

In Section 4.3, we introduced a modified force law for impacts into granular material, equation ( 14), that arises from consideration of final depth scaling (equation 1; U2003) and energy balance (equation 13). Here, we demonstrate its validity by comparing the instantaneous penetration dynamics of a spherical aluminium projectile into dry sand (Section 3) to the behavior predicted by the modified granular force law (equation 14). These simulations present an adequate test bed as the full penetration dynamics were simulated (from impact to projectile stopping), and the dynamics were not influenced by contributions from application-specific mechanisms (like a constant force spring). In Fig. 14, we plot the instantaneous projectile speed, U(t), and penetration depth, Z(t), for Cases 10, 11, and 12 (see Table 2). The behaviour of the model was obtained by numerically integrating equation ( 14). Furthermore, we use a modified drag coefficient for spherical impactors based on KD2013. In Fig. 14(a), we see that U (t) in the simulations is well approximated by equation ( 14), particularly at low instantaneous speeds. In Fig. 14(b), we show that the simulated and analytically derived value of Z(t) also have a good match. For the comparison of the simulated and analytical Z(t), we have used our findings in Section 3.4.3 to modify the instantaneous depth by the impactor-to-target dimensionless parameter: (r/a)1/4. Overall, the modified force law (equation 14) is able to reproduce key characteristics of the penetration dynamics exhibited by the simulations, such as the stopping time and final penetration depth.

(a) Comparison of simulation data of impacts of aluminium spheres into dry sand. We show the speed of the projectile as a function of time for Cases 10 (red crosses), 11 (green triangles), and 12 (blue circles). The black dashed lines is the modeled instantaneous projectile speed obtained by numerically integrating Eq. ( 14). (b) The same as (a), but we compare the instantaneous penetration depth of the projectile.

In this study, we showed that the force laws that govern the interaction between TAGSAM and the surface of an asteroid can vary depending on the geotechnical properties of the regolith. Here, we have shown that (i) packing fraction, P, (ii) bulk density, ρ, (iii) interparticle cohesion, σc, and (iv) angle of friction, ɸ, all govern whether a slowly impacting object rests on the surface, penetrates into the interior, or ricochets. In particular, the packing fraction directly influences the latter three geotechnical properties. As such, packing fraction has a non-linear effect in a regolith's ability to resist penetration. For surfaces that have P ≲ 0.4, the TAGSAM-regolith interaction can be described by Eq. ( 14), and the accelerations felt by the spacecraft would be mostly determined by the regolith's P and ρ. For surface that have P ≳ 0.6, the TAGSAM-regolith interaction is most sensitive to regolith strength properties, σc and ɸ. In this regime, the constant force spring is triggered for ɸ ≳ 28° or σc ≳ 50 Pa, and we find that the spring prevents the immediate bouncing of the spacecraft from the surface, as happened for the Philae lander on comet 67P/C-G (Biele et al. 2015) and for the Hayabusa spacecraft on asteroid Itokawa (Yano et al. 2006). In the case of cohesionless material, high frictional static and rolling contacts are established through motion loading from the spacecraft, allowing the regolith to transition to a jammed state. In this scenario, a spacecraft without a constant force spring would ricochet off the surface.

For a nominal spacecraft approach speed U = 10 cm s−1, a bulk density and angle of repose derived from observations of Bennu, ρ = 1.2 g cm−3 and μ = tan(40°) (Barnouin et al. 2019, Scheeres et al. 2019; Barnouin et al. 2021). For an elastic medium, γ would be the inverse of the Poisson's ratio, and γ = 10 would be equivalent to a Poisson's ratio of 0.1, which is well within the range of possibilities for terrestrial materials. Based on this assessment and the discussion above regarding the relative strengths of Lunar and Martian soil, we assume γ = 10, giving σcrit = 24.3 Pa, which is the same order of magnitude as the value derived from our simulations. This calculation signifies that the cohesion on the surface of Bennu needs to be larger than 24.3 Pa for the TAGSAM penetration dynamics to be dominated by the cohesive interaction of the surface particles rather than their frictional interaction. This cohesive strength is larger than the value estimated for bulk Bennu (Barnouin et al. 2019; Scheeres et al. 2019) and surface stability assessments (Barnouin et al. 2021). Estimates of the strength of Bennu and Ryugu's deeper sub-surface have also been obtained based on crater-scaling relationships. Perry et al. ( 2021) find that cratering strength <100 Pa on Bennu's surface based on the presence of crater ejecta. Furthermore, Arakawa et al. ( 2020) find that the strength of the first subsurface layer of Ryugu (more than 1 m below the surface) is ∼100 Pa based on results from the small carry-on impactor (SCI). We note that these values of asteroid strength reflect the cratering strength, which can be a combination of tensile, shear, and compressive strength.

There is a depth dependence of U as the impactor penetrates into the asteroid surface and decelerates. If one simply considers U = 10 cm s−1, and assuming δ = 54 g cm−3 and g = 7e-5 m s−2 based on the surface gravity at the OSIRIS-REx sampling site (Daly et al. 2020, Lauretta et al. 2021), then zcrit = 287 m, which is larger than Bennu's radius. Instead, we numerically integrate the penetration dynamics for a cohesionless regolith, keeping track of the relative magnitudes of the second and third terms of equation ( 17). When these two terms are equivalent, we find that zcrit ∼ 0.6 m, when U = 0.3 cm s−1, signifying that gravity contributes negligibly to slowing down the spacecraft, as most of the energy has been dissipated by the time gravity has a significant effect.

Historically, knowledge of the response of planetary surfaces to an impact has been driven by a desire to understand the cratering (e.g. Holsapple 1993) and catastrophic disruption (e.g. Michel & Ballouz et al. 2020) processes. Compared to impact cratering, landing or touching down on an asteroid is an extremely low-energy event. For example, a 50-m-radius asteroid striking the Moon at 10 km s−1 has a kinetic energy of ∼1017 J. In contrast, the OSIRIS-REx spacecraft touching down on Bennu has a kinetic energy that was lower by 16 orders of magnitude. Therefore, the relevant forces that govern the outcomes of low-energy events are different than that of impact cratering, and our ability to predict their outcomes is, in contrast, immature.

Low-speed interactions with regolith surfaces have been studied mainly through laboratory experiments of impacts into granular matter in Earth gravity ( U2003, KD2013) and in environments that can simulate low-gravity interactions for a limited time (Brisset et al. 2018; Murdoch et al. 2021). As we have demonstrated, numerical simulations of granular material are an invaluable tool for exploring these low-energy impact regimes, as we are able to simulate the low-gravity environment of small asteroids easily and recreate the detailed design characteristics of spacecraft devices or landers (Maurel et al. 2018; Thuillet et al. 2018; Çelik et al. 2019, 2021). Combined with observations of the surface and telemetry from spacecraft proximity operations, numerical simulations can provide new insights into the behavior of granular materials in low gravity.

The exploration of the small near-Earth asteroids Bennu and Ryugu by the OSIRIS-REx and Hayabusa2 missions, respectively, has revealed natural forms of low-energy interactions on asteroid surfaces, underscoring the value of extracting geotechnical properties from spacecraft-to-surface interactions. Over the course of the OSIRIS-REx proximity operations at Bennu, particles were observed to spontaneously eject from the surface of Bennu with minimum energies between 8 and 270 mJ, corresponding to objects with diameters between millimetres and 10 cm and speeds of about 0.05–3 m s−1 (Lauretta & Hergenrother et al. 2019). One particle was observed to ricochet off the surface of the asteroid, similar to simulations and observations of small landers interacting with regolith in low-gravity environments (Çelik et al. 2019; Thuillet et al. 2021), which enabled a coefficient of restitution measurement on the Bennu surface of 0.57 ± 0.01 (Chesley et al. 2020). Plausible mechanisms for the ejection events are meteoroid impacts (Bottke et al. 2020), thermal fatigue (Molaro et al. 2020b), and phyllosilicate dehydration (Lauretta & Hergenrother et al. 2019). Combined with direct and geologic evidence of the production of low-speed ejecta from cratering events (Arakawa et al. 2020; Perry et al. 2021), it is evident that low-energy events contribute to the production of regolith and the physical transformation of the asteroid surface.

Because the approach speed of a spacecraft undergoing a controlled descent to the surface of a low-gravity asteroid is similar to that asteroid’s escape speed, the impact dynamics of crater ejecta should be similar to that of spacecraft and landers. In particular, low-speed impact dynamics have been used to estimate the burial depth of boulders on the surfaces of asteroids (Nakamura et al. 2013, Wright et al. 2020). The global values of the four geotechnical properties mentioned in Sec. 5.1 (P, ρ, σc, and ɸ) may be inferred if the geopotential properties of the asteroid are known (e.g. Barnouin et al. 2019; Scheeres et al. 2019; Hirabayashi et al. 2020, 2020); however, heterogeneity in the asteroid interior and surface can lead to large variations in these values, and gravity inversion techniques may lead to non-unique solutions (Scheeres et al. 2020). As such, one manner in which the geotechnical properties at a landing site of interest may be estimated, prior to touchdown, is by considering the burial depth of boulders in its vicinity and using equation ( 17) to constrain the possible values of P, ρ, σc, and ɸ. For example, for a cohesionless regolith where ρg = 2.2 g cm−3 and ɸ = 40°, a 1-m-diameter boulder impacting the Bennu surface vertically with U = 20 cm s−1 would be buried at a depth d = 0.15, 0.25, and 0.75 m for P = 0.5, 0.3, and 0.1, respectively. Ejecta are unlikely to impact the surface vertically, but this example illustrates the sensitivity of the burial depth of boulders to packing fraction. Furthermore, our code validation work, in Section 3, showed that an additional contributor to penetration dynamics of a boulder is the ratio of impactor size to particle size at the impact location. Though the sensitivity to penetration depth is fairly weak, it may play an important role for the largest boulders on asteroids, preventing them from effectively sinking below the surface if ejected at speeds below the escape speed. This framework may explain the observation of boulders perching on the surface of Itokawa (Miyamoto et al 2007; Wright et al. 2020) and the Hayabusa’s spacecraft bouncing interaction with the surface (Yano et al. 2006). On Bennu, there appear to be boulders both perched and lying underneath the local surface (Walsh et al. 2019; Daly et al. 2020; Jawin et al. 2020). Some of these boulders have been interpreted to be both buried and exhumed by surface mass-movement processes (Daly et al. 2020; Jawin et al. 2020); however, sites of potentially recent secondary impacts have also been identified (Perry et al. 2021). As such, the variable burial depths of boulders on Bennu potentially suggest that touching down on different parts of its surface would lead to different outcomes.

The granular force law we provide in equation ( 17) can be combined with spacecraft telemetry of the touchdown event that occurred in October 2021 to infer the properties of Bennu’s surface and near sub-surface. Furthermore, these equations can be utilized for future missions to small regolith-covered bodies. For future missions that will explore the surfaces of small Solar system bodies, such as the Martian Moons Exploration mission to Phobos (Kawakatsu et al. 2017) and the Hera mission to the Didymos binary system (Michel et al. 2018), particular attention to the burial depth of boulders may provide additional information on the local geotechnical properties of the regolith.

This material is based upon work supported by NASA under Contract NNM10AA11C issued through the New Frontiers Program. KJW was supported in part through the NASA Solar System Exploration Research Virtual Institute node Project ESPRESSO, cooperative agreement number 80ARC0M0008. P. Michel and Y. Zhang acknowledge support from CNES, Doeblin Federation, the UCA IDEX JEDI ‘Individual grants for young researchers’ programme and the European Union’s Horizon 2020 research and innovation program under grant agreement no. 870377 (project NEO-MAPP). We are grateful to the entire OSIRIS-REx Team for making the encounter with Bennu possible.

Simulation output data and platform-specific versions of the pkdgrav and gdc-i codes will be made available by the lead author, R-LB, at reasonable request. The pkdgrav implementation is fully described in Ballouz ( 2017). The gdc-i implementation is fully described in Sánchez et al. ( 2013).

Altshuler E. et al. , 2014, Geophys. Res. Lett., 41, 3032 10.1002/2014GL059229

Arakawa M. et al. 2020, Science, 368, 67 10.1126/science.aaz1701

Ballouz R.-L. , 2017, PhD thesis, Univ. Maryland, College Park

Ballouz R.-L. et al. 2020, Nature, 587, 205 10.1038/s41586-020-2846-z

Barnouin O.S. et al. 2019, Nat. Geosci., 12, 247 10.1038/s41561-019-0330-x

Barnouin O.S. et al. , 2021, J. Geophys. Res. Planets, in review

Barnouin-Jha O.S. et al. , 2008, Icarus, 198, 108 10.1016/j.icarus.2008.05.026

Basilevsky A.T. et al. , 2014, Planet. Space Sci., 117, 312 10.1016/j.pss.2015.07.003

Biele J. et al. , 2009, Acta Astronaut., 65, 1168 10.1016/j.actaastro.2009.03.041

Biele J. et al. 2015, Science, 349, aaa9816–1-6.

Bierhaus E.B. et al. 2018, Space Sci. Rev., 214, 107 10.1007/s11214-018-0521-6

Bierhaus E.B. et al. , 2021, Icarus, 355, 114142

Bottke W.F. et al. , 2020, J. Geophys. Res. Planets, 125, e2019JE006282–1-16. 10.1029/2019JE006282

Brisset J. et al. 2018, Progr. Earth Planetary Sci., 5, 73 10.1186/s40645-018-0222-5

Burke K.N. et al. , 2021, Remote Sensing, 13, 1315 10.3390/rs13071315

Cates M.E. et al. 1998, Phys. Rev. Lett., 81, 1841 10.1103/PhysRevLett.81.1841

Çelik O. et al. , 2019, Planet. Space Sci., 178, 104693

Cheng A.F. , 2002, in Bottke W. F. Jr. , Cellino A. , Paolicchi P. , Binzel R. P. , eds., Asteroids III. Univ. Arizona Press, Tucson, p. 351

Chesley S.R. et al. , 2020, J. Geophys. Res. Planets, 125, e2019JE006363–1-33. 10.1029/2019JE006363

Clark A.H. et al. , 2012, Phys. Rev. Lett., 109, 238302 10.1103/PhysRevLett.109.238302

Colwell J. , 2003, Icarus, 164, 188 10.1016/S0019-1035(03)00083-6

Costello E.S. et al. , 2018, Icarus, 314, 327 10.1016/j.icarus.2018.05.023

Cundall P.A. , Strack O.D.L. , 1979, Geotechnique, 29, 47

Daly M.G. et al. 2020, Sci. Adv., 6, eabd3649–1-11.

Delbo M. et al. 2014, Nature, 508, 233 10.1038/nature13153

DellaGiustina D.N. et al. 2019, Nature Astron., 3, 341 10.1038/s41550-019-0731-1

Dunham D.W. et al. , 2002, Icarus, 159, 433 10.1006/icar.2002.6911

Dombard A.J. et al. 2010, Icarus, 210, 713 10.1016/j.icarus.2010.07.006

El Mir C. et al. , 2019, Icarus, 333, 356 10.1016/j.icarus.2019.06.001

Feng Y. et al. , 2019, Soft Matter, 15, 3008 10.1039/C8SM02480D

Fujiwara A. et al. 2006, Science, 312, 1330 10.1126/science.1125841

Geissler P. et al. , 1996, Icarus, 120, 140 10.1006/icar.1996.0042

Hirabayashi M. et al. 2020, Icarus, 352:

Heiken G.H. , Vaniman D.T. , Bevan M.F. (eds) 1991, The Lunar Source Book. Cambridge Univ. Press, Cambridge

Holsapple K.A. , 1993, Annu. Rev. Earth Planet. Sci., 21, 333 10.1146/annurev.ea.21.050193.002001

Holsapple K.A. , 1994, Planet. Space Sci., 42, 1067 10.1016/0032-0633(94)90007-8

Housen K.R. , Holsapple K.A. 2011, Icarus, 211, 856 10.1016/j.icarus.2010.09.017

Hörz F. et al. , 1975, Moon, 13, 235 10.1007/BF00567517

Jaumann R. et al. 2019, Science, 365, 817 10.1126/science.aaw8627

Jawin E. R. et al. 2020, J. Geophys. Res. Planets, 125, e2020JE006475–1-21. 10.1029/2020JE006475

Katsuragi H. , Durian D. J. , 2007, Nat. Phys., 3, 420

Katsuragi H. , Durian D.J. , 2013, Phys. Rev. E, 87, 052208–1-5. (KD2013)

Kawakatsu Y. et al. , 2017, In 68th International Astronautical Congress, Vol. 17. International Astronautical Federation, Australia, p. A3.3A.5

Lauretta D. S. et al. , 2017, Space Sci. Rev., 212, 925

Lauretta D.S. , DellaGiustina D.N. et al. , 2019, Nature, 568, 55 10.1038/s41586-019-1033-6

Lauretta D.S. , Hergenrother C.W. et al. , 2019, Science, 366, eaay3544–1-10.

Lauretta D. S. et al. , 2021, In Longobardo A. , ed., Sample Return Missions. Elsevier, Amsterdam

Lin J. , Wei W. , 2012, Philos. Mag., 92, 3474 10.1080/14786435.2012.706373

Liu A. , Nagel S.R. , 1988, Nature, 396, 21 10.1038/23819

Macke R.J. , Consolmagno G.J. , Britt D.T. , 2011, Meteorit. Planet. Sci., 46, 1842 10.1111/j.1945-5100.2011.01298.x

Majmudar T.S. et al. , 2007, Phys. Rev. Lett., 98, 058001 10.1103/PhysRevLett.98.058001

Maurel C. et al. , 2018, Adv. Space Res., 62, 2099 10.1016/j.asr.2017.05.029

Metcalf J.R. , 1966, Int. J. Rock Mech. Mining Sci. Geomech. Abstracts, 3, 155

Michel P. et al. 2018, Adv. Space Res., 62, 2261 10.1016/j.asr.2017.12.020

Michel P. , Ballouz R.-L. et al. 2020, Nat. Commun., 11, 2655–1-11. 10.1038/s41467-020-16433-z

Miyamoto H. et al. 2007, Science, 316, 1011 10.1126/science.1134390

Miyai S. et al. 2019, Granular Matter, 21, 105

Molaro J. L. et al. , 2020a, Nat. Commun., 11, 1 10.1038/s41467-020-16528-7

Molaro J. L. et al. , 2020b, J. Geophys. Res. Planets, 125, e2019JE006325 10.1029/2019JE006325

Mukhopadhyay S. et al. , 2018, Phys. Rev. E, 97, 052602 10.1103/PhysRevE.97.052602

Murdoch N. et al. , 2017, MNRAS, 468, 1259 10.1093/mnras/stw3391

Murdoch N. , 2021, Mon. Not. R. Astron. Soc., 503, 3460 10.1093/mnras/stab624

Nakamura A. M. et al. , 2013, Icarus, 223, 222 10.1016/j.icarus.2012.11.038

Nakamura A.M. , 2017, Planet. Space Sci., 149, 5 10.1016/j.pss.2017.09.001

O'Hern C.S. et al. , 2002, Phys. Rev. Lett., 88, 075507 10.1103/PhysRevLett.88.075507

Perry M. et al. , 2021, Nat. Geosci., in review

Perko H.A. et al. , 2001, J. Geotech. Geoenviron. Eng., 127, 371

Richardson D. C. et al. , 2000, Icarus, 143, 45 10.1006/icar.1999.6243

Richardson D. C. et al. , 2009, Planetary Space Sci., 57, 183 10.1016/j.pss.2008.04.015

Richardson D. C. et al. 2011, Icarus, 212, 427 10.1016/j.icarus.2010.11.030

Robinson M.S. et al. 2002, Meteorit. Planet. Sci., 37, 1651 10.1111/j.1945-5100.2002.tb01157.x

Royer J.R. et al. , 2011, EPL (Europhys. Lett.), 93, 28008 10.1209/0295-5075/93/28008

Richter L.O. , 2005, Fall Meeting 2005. American Geophysical Union, San Francisco, CA, P21A-0134

Rycroft C. H. , 2009, Chaos: An Interdisciplinary J. Nonlinear Sci., 19, 041111

Sánchez P. , Scheeres D.J. , 2011, ApJ, 727, 120 10.1088/0004-637X/727/2/120

Sánchez P. et al. , 2013, 44th Lunar and Planetary Science Conference. Lunar and Planetary Institute, Texas

Sánchez P. , Scheeres D.J. , 2014, Meteorit. Planet. Sci., 49, 788 10.1111/maps.12293

Sánchez P. , Scheeres D.J. , 2016, Icarus, 271, 453 10.1016/j.icarus.2016.01.016

Sánchez P. , Scheeres D.J. , 2018, Planet. Space Sci., 157, 39 10.1016/j.pss.2018.04.001

Scheeres D.J. et al. 2019, Nature Astron., 3, 352 10.1038/s41550-019-0721-3

Scheeres D.J. et al. 2020, Sci. Adv., 6, eabc3350–1-16.

Schwartz S. et al. , 2012, Granular Matter, 14, 363

Schwartz S. R. et al. , 2013, Icarus, 226, 67 10.1016/j.icarus.2013.05.007

Seguin A. et al. 2008, Phys. Rev., E78, 010301 10.1103/PhysRevE.78.010301

Sugita S. et al. , 2019, Science, 364, 252 10.1126/science.aaw0422

Thuillet F. et al. , 2018, A&A, 615, A41 10.1051/0004-6361/201832779

Thuillet F. et al. , 2021, A&A, 648, A56 10.1051/0004-6361/201936128

Uehara J. S. et al. 2003, Phys. Rev. Lett., 90, 194301 (U2013) 10.1103/PhysRevLett.90.194301

Umbanhowar P. , Goldman D.I. , 2010, Phys. Rev. E, 82, 010301 10.1103/PhysRevE.82.010301

Virtanen P. et al. , 2020, Nat. Methods, 17, 261

Wada K. et al. , 2006, Icarus, 180, 528 10.1016/j.icarus.2005.10.002

Walsh K.J. et al. , 2019, Nat. Geosci., 12, 242 10.1038/s41561-019-0326-6

Williams J. et al. , 2006, Open-File Report 2006-1195. United States Geological Survey, Virginia

Wright E. et al. 2020, Icarus, 351, 113963

Yano H. et al. , 2006, Science, 312, 1350 10.1126/science.1126164

Yu Y. et al. , 2014, Icarus, 242, 82 10.1016/j.icarus.2014.07.027

Zacny K. et al. , 2018, in Abreu N. , ed. Primitive Meteorites and Asteroids. Elsevier, Amsterdam, p. 439

Zhang Y. et al. 2017, Icarus, 294, 98 10.1016/j.icarus.2017.04.027

Zhang Y. et al. , 2018, ApJ, 857, 15 10.3847/1538-4365/aa9f28

Oxford University Press is a department of the University of Oxford. It furthers the University's objective of excellence in research, scholarship, and education by publishing worldwide

This PDF is available to Subscribers Only

For full access to this pdf, sign in to an existing account, or purchase an annual subscription.

Your Trustworthy Specialist And Solution Provider

Excellent Performance, Accurate Dimension, Fast Delivery

100% On-Time Shipment Protection

Modified granular impact force laws for the OSIRIS-REx touchdown on the surface of asteroid (101955) Bennu | Monthly Notices of the Royal Astronomical Society | Oxford Academic

Featured Products

News & Blog

Flow Wrap Machine Market Size, Scope, Growth, Competitive Analysis – Bosch, FUJI MACHINERY, Coesia Group, PFM Packaging Machinery – The Sabre

North America Holds Major Market Share in Global Flow Wrap Machines Industry - EIN Presswire

Tray Packing Machine Market 2022 Business Outlook, SWOT Analysis, Key Business Strategies, Industry Players and Forecast 2026 – KHS GmbH, IMA S.p.A., OYSTAR Holding GmbH – ChattTenn Sports

Newsletters

Sintratec aims to take 3D printing software "to the next level" with new nesting update - 3D Printing Industry

Newsletters

Whitesnake Discussing New Album to Follow Farewell Tour

How to Pack Like a Seasoned Pro for Super-Long Trips - InsideHook

13 Chocolate Egg Hacks For Easter

Best Chip Brands Ranked