2. General governing equations

Figure 2 represents the setting in the context of an initially axisymmetric body in figure 2(a) that is rotating at the angular velocity $\boldsymbol \omega$ . The body is taken to comprise grains that, at any given time, may be broadly distinguished into two regions: a solid central body where grains are consolidated and/or compacted and/or sintered, and a mobile layer of regolith that arises when landslides occur, as shown in figure 2(b). On the other hand, in between two landsliding events the entire body will be the central body. The shape and size of the central body is thus modified when regolith is mobilised, flows and stops. We incorporate the possibility of multiple landslides initiated, say, by impact events over the course of the body’s lifetime, and the inset in figure 2 shows how these landslides deposit over each other.

Figure 2.

Schematic of a rotating planetary body with angular velocity $\boldsymbol{\omega }$ . Three different coordinate systems are indicated and these are discussed in the main text. The differential volume element ${\text{d}} V$ of an axisymmetric regolith layer is shown. (b) The shape of an initially spherical body after multiple global axisymmetric landsliding events. The brown layer depicts an ongoing landslide (regolith flow). The inset shows how landslides are deposited over time over the central body.

We now proceed as follows. We will first gather equations that govern regolith flow on a body and couple it to the body’s rotational dynamics in a fairly general setting. The flow will be taken to be three-dimensional. We will allow the body’s angular velocity to be arbitrarily located, and not necessarily aligned with the body’s principal axes of inertia; thus the body may tumble in space. Furthermore, the landslides need not be axisymmetric and thus cause the centre of mass of the body to shift. The rotational dynamics of the system will incorporate effects of the system’s changing mass moment of inertia. No specific assumptions will be made on forces and torques acting on the system. We will subsequently introduce assumptions that will simplify the presentation, but will still be useful for realistic applications.

We will employ three coordinate systems shown in figure 2(a), viz. the inertial coordinate system $\{ \mathcal{P},P,\hat {\boldsymbol{E}}_i\}$ with origin at some space-fixed point $P$ and unit vectors $\hat {\boldsymbol E}_i$ ; the body-fixed coordinate system $\{\mathcal{E}_0,C,\hat {\boldsymbol{e}}_i\}$ , which is attached to the central body and is defined by its principal axes $\hat {\boldsymbol{e}}_i$ and $C$ lies at the central body’s centre of mass; and $\{ \mathcal{E},\, R,\, \hat {\boldsymbol{g}}_i\}$ is a curvilinear coordinate system defined by the natural unit basis $\hat {\boldsymbol g}_i$ associated with the body’s surface. Henceforth, we represent all vector quantities measured in the inertial frame $\mathcal{P}$ and the body-fixed coordinate system by, respectively, bold upper and lower case letters.

A common simplification in the analysis of landslides (Pudasaini & Hutter Reference Pudasaini and Hutter2007, Chap. 2.7.3) is to assume that the flowing regolith constitutes an incompressible continuum. This assumption rests upon several observations. Large terrestrial landslides have a nearly plug-like flow (Melosh Reference Melosh1986), so that shearing occurs only in a thin layer near the base. Consequently, the dilatation that accompanies granular flows is confined to a small region of the flow. Later, we will depth average our equations, so that ignoring these volume changes will not affect our results in a significant manner. The mass balance and linear momentum balance for an incompressible flow in the body-fixed coordinate system $ {\mathcal{E}_0}$ are, respectively,

(2.1) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0 \quad \text{and} \quad \rho \frac {\partial {\boldsymbol{u}}}{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{u}\otimes \boldsymbol{u})&=\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }+\rho (\boldsymbol{b}-2\boldsymbol{\omega }\times \boldsymbol{u}-\boldsymbol{r}\times \boldsymbol{\alpha }-\boldsymbol{A}_c), \end{align}

where $\boldsymbol{u}$ is the flow’s velocity relative to $\mathcal{E}_0$ , $\rho$ is the regolith’s mean density, $\boldsymbol{\sigma }$ is the stress tensor, $\boldsymbol \omega$ and $\boldsymbol{\alpha }$ are, respectively, the central body’s angular velocity and angular acceleration, $\boldsymbol{r}$ is a mass element’s position vector relative to $C$ , $\boldsymbol{b}=\boldsymbol{b}^0-\boldsymbol{\omega }\times (\boldsymbol{\omega }\times \boldsymbol{r})$ is the total effective gravity in terms of the central body’s gravity $\boldsymbol{b}^0$ and the centrifugal acceleration $\boldsymbol{\omega }\times (\boldsymbol{\omega }\times \boldsymbol{r})$ , $2\boldsymbol{\omega }\times \boldsymbol{u}$ is the Coriolis acceleration and $\boldsymbol{A}_C$ is the inertial acceleration of the central body’s centre of mass. To complete (2.1) we need to prescribe a rheology appropriate for the granular continuum that relates the stress tensor to the flow’s strain rate, and we return to this in the next section.

Next, we consider how the central body’s dynamics couples to the regolith flow. In the inertial coordinate system $\mathcal{P}$ , the linear momentum balance of the system comprising of the central body and the flowing regolith is

(2.2) \begin{equation} \frac {{\text{d}}}{{\text{d}} t}({M_c\boldsymbol{U}_{\!c}})+\frac {{\text{d}}}{{\text{d}} t}({M_r\boldsymbol{U}_{\!r}})=\boldsymbol{F}-\boldsymbol{F}_{\!s}, \end{equation}

where the subscripts ‘ $c$ ’ and ‘ $r$ ’ indicate the central body and the flowing regolith, respectively, $M$ is the mass, $\boldsymbol{U}$ is the centre of mass velocity, $\boldsymbol{F}$ represents external forces on the system and $\boldsymbol{F}_{\!s}$ is the rate of momentum carried away by the mass shed from the body’s surface. Due to erosion and deposition during landsliding, mass is exchanged between the stationary regolith (central body) and the flowing grains, so that $M_c$ and $M_r$ may be modified. Additionally, $M_r$ may also change due to surface mass shedding.

Finally, in the inertial coordinate system $\mathcal{P}$ , the global angular momentum balance of the system about the point $P$ is

(2.3) \begin{equation} \frac {{\text{d}} \boldsymbol{H}_{\!c}}{{\text{d}} t}+\frac {{\text{d}}\boldsymbol{H}_{\!r}}{{\text{d}} t}=\boldsymbol{T}- \boldsymbol{T}_{\!s}, \end{equation}

where the subscripts ‘ $c$ ’ and ‘ $r$ ’ are as above, $\boldsymbol{H}$ represents angular momentum, $\boldsymbol{T}$ is the external torque and $\boldsymbol{T}_{\!s}$ is the rate at which angular momentum is carried away by the shed mass. The above angular momentum balance is simplified further in Appendix A to obtain the following form:

(2.4) \begin{equation} \boldsymbol{\omega }\times (\boldsymbol{I}\boldsymbol{\cdot }\boldsymbol{\omega }+\boldsymbol{h}_r)+\boldsymbol{I}\boldsymbol{\cdot }\boldsymbol{\alpha }+\dot {\boldsymbol{h}}_r+\dot {\boldsymbol{I}}\boldsymbol{\cdot }\boldsymbol{\omega }+ \boldsymbol{R}_{O} \times M\!\boldsymbol{A}_c ={\boldsymbol{T}}-{\boldsymbol{T}}_s, \end{equation}

where $\boldsymbol{I}$ is the moment of inertia about $C$ of the system comprised, in general, a static part (central body) and flowing regolith, $\boldsymbol{h}_r$ is the angular momentum of the flowing regolith about $C$ in the body-fixed coordinate system $\mathcal{E}_0$ , $\dot {(\, )}$ represents time differentiation in the body-fixed coordinate system and $M$ and $\boldsymbol{R}_O$ are, respectively, the total system mass and the location of its centre of mass with respect to $C$ . The first term on the left-hand side of (2.4) arises from the tumbling motion of the central body, the second is the rotational inertia and the third and fourth terms are due to regolith flow and the final term is because the centre of mass of the system may not be at the origin $C$ , which introduces an extra contribution from the centre of mass’s acceleration.

Thus far, while developing our modelling framework, we have made no assumptions regarding forces and torques affecting our system, nor have we imposed any constraints on the kinematics of the central body or the nature of landslides. We may, in principle, implement this framework computationally. However, in the next section we introduce several physically motivated simplifications that make the presentation accessible, without sacrificing essential details or negating practical application.

Two remarks are in order before we close this section. First, when landslides occur the body consists of a solid central body and a flowing regolith layer. Thus, we cannot strictly ascribe to the body an angular velocity, at least in the sense of a rigid body. However, landsliding events are intermittent, with each landslide remaining active for a duration of time much shorter than the rotational dynamics time scale and, moreover, landslides are typically of depth much smaller in comparison with the body’s size. Thus, we may identify the central body’s rotational motion with that of the body itself. Second, as mentioned, we believe that landslides are most commonly initiated by seismic activity following an impact event. Here, we will not concern ourselves with how the location and amount of mobilised regolith is related to the details of the corresponding impact. Instead, in the sections to follow we will make simplifying assumptions and leave detailed modelling of this aspect to future work. This is done to retain focus in this work on how landslides affect, and are affected by, the shape and rotational motion of small rubble planetary bodies.

3. A simplified model

Figure 1 shows asteroids Bennu and Ryugu that were the targets of recent missions. Their shapes are observed to be fairly axisymmetric and the objects were found to be in pure spin about the symmetry axis. This suggests specialising our model to an axisymmetric situation, wherein the initial body and all subsequent landslides are axisymmetric about a space-fixed rotation axis. We further assume that the landslides are shallow, in that their depth normal to the body’s topography is much smaller than the extent of their spread over the body’s surface. Consequently, the momentum carried by the regolith flow and shed mass is much smaller than the system’s orbital momentum, which then balances the Sun’s gravitational force on the right in (2.2). The angular acceleration arising from landslide-induced changes in the body’s moment of inertia is nevertheless retained. This is because, for a spherical shell of thickness $H$ and diameter $D$ , the mass scales as $D^{2}H$ , whereas the moment of inertia scales as $D^{4}H$ . Indeed, Banik et al. (Reference Banik, Gaurav and Sharma2022) found that, although the change in the body’s rotation due to each landslide over its lifetime is small, the cumulative effect is significant. This is similar to the action of radiation torques discussed below, and makes the case for the present modelling that couples these torques with landsliding.

Now, the external forces $({\boldsymbol F})$ and moments $({\boldsymbol T})$ entering into the governing (2.2) and (2.4) arise from gravitational attraction of the Sun and other bodies, (thermal) radiation pressure and radiation torques. Radiation pressure – the Yarkovsky effect – influences orbital motion, while radiation torques – also called YORP torques, or simply YORP – change the body’s rotation, especially over long time periods; please see Vokrouhlický et al. (Reference Vokrouhlický, Bottke, Chesley, Scheeres and Statler2015) for more information. Here, we limit ourselves to small isolated objects. Thus, the Sun’s quadrupole torque – torque due to varying Solar gravitational forces experienced at different parts of the body – is negligible, as are tidal interactions. Furthermore, the orbital drift of sub-kilometre-sized asteroids caused by the Yarkovsky effect is typically of the order of $10^{-4}{-}10^{-3}$ AU/Myr, where AU is the astronomical unit and Myr represents a million years (Vokrouhlickỳ et al. Reference Vokrouhlickỳ, Milani and Chesley2000). Therefore, the orbital change – and consequently its influence on the YORP torque – may be neglected when simulating over a time scale of a few million years. Along with our assumption above regarding the balance of orbital inertia and the Sun’s gravitational force in (2.2), these considerations decouple the angular momentum balance (2.4) from (2.2). Thus, we may study landslides independently of the body’s orbital motion and external gravity forces.

The angular momentum balance (2.4) is further simplified for bodies that are nominally in a state of pure spin, i.e. bodies with a space-fixed rotation axis $\hat {\boldsymbol{e}}_3$ , which is generally their axis of greatest moment of inertia. This is true for many small planetary objects, including the asteroids in figure 1. As mentioned in the introduction, internal friction damps out tumbling motion in bodies in free space, driving them into a state of pure spin about their axis of greatest inertia (Burns et al. Reference Burns, Safronov and Gold1973; Sharma et al. Reference Sharma, Burns and Hui2005; Breiter et al. Reference Breiter, Rożek and Vokrouhlickỳ2012).

For a body in pure spin about $\hat {\boldsymbol e}_3$ we have $\boldsymbol{\omega }=\omega \hat {\boldsymbol{e}}_3$ and $\boldsymbol{\alpha }=\alpha \hat {\boldsymbol{e}}_3.$ Invoking this and the simplifications above, we arrive at the following coupled equations for shallow, incompressible regolith flow on a planetary body in pure spin:

(3.1) \begin{align} && &&\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}&=0, \\[-12pt]\nonumber\end{align}

(3.2) \begin{align} && && \rho \dot {\boldsymbol{u}}+\boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{u}\otimes \boldsymbol{u})&=\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }+\rho (\boldsymbol{b}-2\boldsymbol{\omega }\times \boldsymbol{u}-\boldsymbol{\alpha }\times \boldsymbol{r}),\\[-12pt]\nonumber\end{align}

and

(3.3) \begin{align} \hspace{1cm}\;\;\, I\alpha + \dot {h}_r + \dot {{I}}{\omega }&= T - {T}_s , && \end{align}

where $I$ and $h_r$ are, respectively, the system’s moment of inertia and the flowing regolith’s angular momentum about the rotation axis $\hat {\boldsymbol e}_3$ and $T-{T}_s = ({\boldsymbol T} - {\boldsymbol T}_s )\boldsymbol{\cdot }\hat {\boldsymbol e}_3$ . We recall that the system here consists of the central body and, when applicable, the flowing regolith. Note also that the angular momentum balance reduces to the scalar equation above when we take the body to be in a state of pure spin.

We also require a rheology of the flowing regolith and boundary conditions (BCs) to close the set (3.1) and (3.2). We introduce the former in § 3.3 below, but discuss the latter now. At the free top surface $F^s(\boldsymbol{r})=0$ of the flow, identified by the superscript $s$ , we impose the vanishing of the traction and the usual kinematic conditions, respectively,

(3.4a) \begin{align} \qquad \qquad\quad \boldsymbol{\sigma }^s\boldsymbol{\cdot }\boldsymbol{n}^s&=0, \end{align}

and

(3.4b) \begin{align} \frac {\partial F^s}{\partial t} +\left (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\right ) F^s &= 0, \end{align}

where $\boldsymbol{n}^s$ is the normal to the free surface.

We assume no erosion of material from the central body during a landslide, i.e. mass is not added from the central body to the flowing regolith across the flow’s base. This is a simplifying but not limiting assumption, and may be relaxed in the future. The flow adheres to the basal topography and the surface traction there is assumed to be governed by Coulomb’s friction law. The mathematical statements of these BCs at the flow’s base are, respectively,

(3.5a) \begin{align} \boldsymbol{u}^b\boldsymbol{\cdot }\boldsymbol{n}^b&=0, \\[-12pt]\nonumber \end{align}

and

(3.5b) \begin{align} \big ({\boldsymbol 1} - \boldsymbol{n}^b\otimes \boldsymbol{n}^b\big )\boldsymbol{\cdot }\boldsymbol{\sigma }^b\boldsymbol{\cdot }\boldsymbol{n}^b &=-\mu \boldsymbol{n}^b\boldsymbol{\cdot }\boldsymbol{\sigma }^b \boldsymbol{\cdot }\boldsymbol{n}^b\frac {\boldsymbol{u}^b}{|\boldsymbol{u}^b|}, \end{align}

where $\mu$ is the basal friction coefficient. Landslides flow over the central body, which comprises the same type of grains, so that we set $\tan ^{-1}\mu = \zeta$ , the regolith’s internal friction angle. Note that (3.5a ) is equivalent to the usual no-penetration condition. Finally, mass shedding is assumed to commence when the pressure at the flow’s base, denoted by the superscript $b$ , becomes tensile in nature, i.e. when

(3.6) \begin{equation} \boldsymbol{n}^b\boldsymbol{\cdot }\boldsymbol{\sigma }^b\boldsymbol{\cdot }\boldsymbol{n}^b\gt 0. \end{equation}

This criterion is analogous to that employed by Scheeres (Reference Scheeres2015) and Sánchez & Scheeres (Reference Sánchez and Scheeres2020), except that they consider the loss of a single grain from the body’s surface, which happens when the normal reaction vanishes.

Now, limiting ourselves to shallow axisymmetric landslides on an initially axisymmetric body allows us to reduce the spatial dimensionality of the governing equations to one, because derivatives in the longitudinal direction vanish and we may depth average our equations in a direction normal to the base. The latter is equivalent to ignoring variations in the regolith flow across its depth. Furthermore, because the regolith flow’s thickness is much smaller compared with the central body’s mean radius $R_{\textit{CB}}$ , we may now introduce the small parameter

(3.7) \begin{equation} \varepsilon = {H}/{R_{\textit{CB}}}\ll 1, \end{equation}

where $H$ is the average thickness of the flow when it is initiated. This allows us to approximate various quantities by ignoring variations that are $\mathcal{O}(\varepsilon )$ or smaller. These steps simplify our mathematical description considerably, as we now see.

3.1. Axisymmetric coordinate system

We will now depth integrate the linear momentum balance (3.2) in a manner similar to Gray et al. (Reference Gray, Wieland and Hutter1999) and Luca et al. (Reference Luca, Tai and Kuo2016). To facilitate this we first introduce the natural axisymmetric coordinate system associated with the body’s shape; see, for example, the schematic in figure 2(b). We will express our governing equations in this coordinate system.

The central body’s axisymmetric surface is defined by its distance $r_{\textit{CB}}=R(\theta )$ from $C$ in terms of the co-latitude $\theta$ . To describe the regolith flow over the central body’s surface, we introduce the contravariant coordinate system $(x^1, x^2, x^3) = (n,\theta ,\phi )$ , where $n$ is the normal distance from the central body’s surface and $\phi$ is the longitude. A material point within the regolith flow is then located by

(3.8a) \begin{align} \boldsymbol{r} = r_{\textit{CB}}\hat {\boldsymbol{e}}_r +n\hat {\boldsymbol n}=R(\theta )\hat {\boldsymbol{e}}_r +n\hat {\boldsymbol n}, \end{align}

where

(3.8b) \begin{align} \hat {\boldsymbol{ n}} = {\left (R\hat {\boldsymbol{e}}_r - R'\hat {\boldsymbol{e}}_\theta \right )}N^{-1}, \end{align}

is the outward unit surface normal, $N=\sqrt {R^2 + {R^{\prime }}^2}$ , $\hat {\boldsymbol{e}}_i$ are the unit vectors in spherical coordinates and prime $(\prime)$ denotes differentiation with respect to $\theta$ .

The natural basis vectors $\boldsymbol{g}_i$ are obtained as $ \boldsymbol{g}_i={\partial \boldsymbol{r} }/{\partial x^i}$

(3.9) \begin{align} \boldsymbol{g}_1 = \hat {\boldsymbol{ n}}, \hspace {0.5cm} \boldsymbol{g}_2 = \frac {\kappa ^b_\theta }{\kappa _\theta }(R'\hat {\boldsymbol{e}}_r+R\hat {\boldsymbol{e}}_\theta )\hspace {0.5cm} \text{and} \hspace {0.5cm} \boldsymbol{g}_3 = \frac {\hat {\boldsymbol{e}}_\phi }{\kappa _\phi }, \hspace {0.5cm} \end{align}

where $\kappa _{\theta }$ and $\kappa _{\phi }=\sqrt {\kappa ^2_{\phi n}+\kappa ^2_{\phi g}}$ are the curvatures of $\theta -$ and $\phi -$ coordinate curves, with $\kappa _{\phi g}$ and $\kappa _{\phi n}$ being, respectively, the geodesic and normal curvatures of the $\phi$ -curves, and the superscript $b$ indicates that the quantity is evaluated at the central body’s surface. The formulae for these curvatures are

(3.10) \begin{align} \kappa ^b_{\theta } = \frac {R^2+2{R^{\prime }}^2-RR''}{N^3}, \hspace {0.25cm} \kappa ^b_\phi = \frac {1}{R\sin \theta },\hspace {0.25cm}\kappa _{\theta } = \frac {\kappa ^b_\theta }{1+\kappa ^b_\theta n}\hspace {0.25cm}\text{and}\hspace {0.25cm} \kappa _{\phi } = \frac {\kappa ^b_\phi }{ 1+ \kappa ^b_{\phi n}n}, \\[-12pt]\nonumber \end{align}

and

(3.11) \begin{align} \kappa ^b_{\phi n} = -\frac {(R\cos {\theta })'}{NR \sin \theta },\hspace {0.25cm} \kappa ^b_{\phi g} = \frac {(R\sin \theta )'}{ N R\sin \theta },\hspace {0.25cm} \kappa _{\phi n} =\frac {\kappa ^b_{\phi n}}{1+\kappa ^b_{\phi n}n} \hspace {0.25cm}\text{and}\hspace {0.25cm} \kappa _{\phi g} =\frac {\kappa ^b_{\phi g}}{1+\kappa ^b_{\phi n}n}. \end{align}

Finally, the radius of curvature for the $\theta$ - and $\phi$ -curves are, respectively,

(3.12) \begin{align} \mathcal{R}_\theta =1/\kappa _\theta \quad \text{and} \quad \mathcal{R}_\phi = 1/\kappa _\phi . \end{align}

The basis vectors ${\boldsymbol g}_i$ calculated in (3.9) are orthogonal. In such bases, employing index notation where possible, the divergence of vector $({\boldsymbol v})$ and tensor $({\boldsymbol T})$ fields have the compact representations

(3.13) \begin{align} &\hspace {-4 cm}\!\!\! \displaystyle {\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{v} =\displaystyle \sum _{i=1}^3{\cfrac {1}{ {g}}}\,{\frac {\partial }{\partial x^{i}}}\left (\frac {v_{i}\,{ {g}}}{\sqrt {g_{\textit{ii}}}}\right )\!, \end{align}

and

(3.14) \begin{align} &\hspace {1.5 cm}\ \ \ \ \ \ \, \displaystyle {\boldsymbol{\nabla }}\boldsymbol{\cdot }\boldsymbol{T} =\displaystyle \sum _{i,j=1}^3{\sqrt {g_{\textit{ii}}}}\left ( \frac {1}{{g}}\frac {\partial }{\partial x^j}\left ( \frac {{g T_{\textit{ji}}}}{\sqrt {g_{\textit{ii}}}\sqrt {g_{\textit{jj}}}}\right ) + \displaystyle \sum _{k=1}^3\varGamma ^{i}_{\textit{jk}}\frac {{T_{\textit{kj}}}}{\sqrt {g_{\textit{jj}}}\sqrt {g_{\textit{kk}}}}\right )\hat {\boldsymbol{g}}_i , \end{align}

where $g_{\textit{ij}}=\boldsymbol{g}_i\boldsymbol{\cdot }\boldsymbol{g}_{\textit{j}}$ is the metric tensor, $g=\sqrt {\text{det} [g_{\textit{ij}} ]}$ (no sum) is its norm, $\hat {\boldsymbol{g}}_i={\boldsymbol{g}_i}/{\sqrt {g_{\textit{ii}}}}$ (no sum) are unit vectors, $v_i=\boldsymbol{v}\boldsymbol{\cdot }\hat {\boldsymbol{g}}_i$ and $T_{\textit{ij}}=\hat {\boldsymbol{g}}_i\boldsymbol{\cdot }\boldsymbol{T}\boldsymbol{\cdot }\hat {\boldsymbol{g}}_{\textit{j}}$ are the covariant components of the vector and tensor fields, respectively, and $\varGamma ^{i}_{\textit{jk}}=\boldsymbol{g}^i\boldsymbol{\cdot }({\partial \boldsymbol{g}_{\textit{j}}}/{\partial x^k} )$ are the Christoffel symbols. For our coordinate system we obtain

(3.15) \begin{align} g_{11} =1, \hspace {0.25cm} g_{22} = \left (\displaystyle \frac {\kappa ^b_\theta }{\kappa _\theta }N\right )^2, \hspace {0.25cm} g_{33} = \displaystyle \frac {1}{\kappa ^2_\phi }, \hspace {0.25cm} g_{\textit{ij}} =0\ \ \forall \ \ i\neq j, \hspace {0.25cm} {g} =\frac {N\kappa ^b_\theta }{\kappa _\theta \kappa _\phi } \hspace {0.2 cm} \text{and}\ \hspace {0.2cm} J=g, \end{align}

where $J$ is the Jacobian of the transformation from Cartesian to the axisymmetric coordinate system.

3.2. Non-dimensionalisation

Before proceeding, we non-dimensionalise our equations as per the following definitions:

(3.16) \begin{equation} \begin{array}{cc} \left \{n,t\right \}=\left \{H n^{*},\sqrt {\displaystyle {R_{\textit{CB}}}/{b}}\,t^{*}\right \}\!; \hspace {0.5cm} \left \{u_n,u_\theta ,u_\phi \right \}=\sqrt {\textit{bR}_{\textit{CB}}}\big \{\varepsilon {u}^*_n,{u}^*_\theta ,{u}^*_\phi \big \};\\[7pt] \left \{\omega ,T\right \}=\left \{\sqrt {b/R_{\textit{CB}}}\,{\omega }^*,\rho b {R^4_{\textit{CB}}}\,{T}^*\right \}\!; \hspace {2mm} \text{and} \hspace {2mm} \sigma _{\textit{ij}}=\rho b H\,{\sigma }^*_{\textit{ij}}, \end{array} \end{equation}

where $H$ is the initial average thickness of the landslide, $b$ is the gravitational acceleration on the surface of a sphere of radius $R_{\textit{CB}}$ and superscript $^*$ identifies non-dimensional quantities; the $\,{}^*\,$ is subsequently dropped for a cleaner notation.

3.2.1. Governing equations

The continuity and linear momentum balance equations, upon using (3.9)–(3.16) and assuming axisymmetry, simplify to, respectively,

(3.17) \begin{align} &\frac {\partial }{\partial n}\left (Ju_n \right )+\frac {\partial }{\partial \theta }\big (\mathcal{R}_\phi u_\theta \big )=0, \\[-12pt]\nonumber \end{align}

(3.18) \begin{align} &\varepsilon \frac {\partial }{\partial t}\left (\textit{Ju}_n \right )+ \frac {\partial }{\partial n}\left (\textit{JA}_{nn}\right )+ \varepsilon \frac {\partial }{\partial \theta }\big ( \mathcal{R}_\phi A_{n\theta }\big )=J\left ( \kappa _\theta A_{\theta \theta }+ \kappa _{\phi n} A_{\phi \phi }+ f_n\right )\!, \\[-12pt]\nonumber \end{align}

(3.19) \begin{align} &\frac {\partial }{\partial t}\big (J^2\kappa _\phi u_\theta \big )+ \frac {\partial }{\partial n}\big (J^2\kappa _\phi A_{ n \theta } \big )+ \frac {\partial }{\partial \theta }\left ( {\textit{JA}_{\theta \theta }}\right ) = \frac {\partial }{\partial \theta }\big (J\kappa _\phi \big )\mathcal{R}_\phi A_{\theta \theta } \, \, . . . \nonumber &&\\ &\hspace {7cm}. . . \, + J^2\kappa _\phi \big (\kappa _{\phi g} A_{\phi \phi } + f_\theta \big ) \\[-12pt]\nonumber \end{align}

and

(3.20) \begin{align} \hspace {-1 cm} \frac {\partial }{\partial t}\left ( J\mathcal{R}_\phi u_\phi \right )+ \frac {\partial }{\partial n}\left ( {J\mathcal{R}_\phi A_{ n \phi }}\right )+ \frac {\partial }{\partial \theta }\big ( \mathcal{R}^2_\phi {A_{\theta \phi }}\big ) = {J \mathcal{R}_\phi f_\phi },&& \end{align}

where

(3.21) \begin{align} A_{nn} = \varepsilon u^2_n -\sigma _{nn}, \hspace {0.25cm} A_{ni} = u_n u_i - \sigma _{ni} \hspace {0.25cm} \text{and} \hspace {0.25cm} A_{\textit{ij}} = u_iu_{\textit{j}} - \varepsilon \sigma _{\textit{ij}}, \end{align}

with $i$ and $j$ representing either $\theta$ or $\phi$ and where $f_i$ are the body force terms

(3.22a) \begin{align} &f_n=b^0_n+\mathcal{R}_\phi {\kappa _{\phi n}}\big (\mathcal{R}_\phi {\omega ^2} +2\omega u_\phi \big ), \\[-12pt]\nonumber \end{align}

(3.22b) \begin{align} &f_\theta =b^0_\theta +\mathcal{R}_\phi {\kappa _{\phi g}}\big (\mathcal{R}_\phi {\omega ^2} +2\omega u_\phi \big ) \\[-12pt]\nonumber \end{align}

and

(3.22c) \begin{align} \hspace {2 cm} \qquad &f_\phi =-{\mathcal{R}_\phi }\left (2\omega (\varepsilon \kappa _{\phi n}u_n+\kappa _{\phi g}u_\theta )+\alpha \right ),&& \end{align}

in terms of $b^0_n = {\boldsymbol b}^0\boldsymbol{\cdot }\hat {\boldsymbol g}_1$ and $b^0_\theta = {\boldsymbol b}^0\boldsymbol{\cdot }\hat {\boldsymbol g}_2$ .

3.2.2. Boundary conditions

The non-dimensionalised kinematic BC at the top and bottom surfaces are given by, respectively,

(3.23a) \begin{align} \hspace {-2 cm}\,\,\,\,\, J^s u^s_n-\frac {u^s_\theta }{\kappa ^s_\phi }\frac {\partial h}{\partial \theta }&=J^s\frac {\partial h}{\partial t} \\[-12pt]\nonumber \end{align}

and

(3.23b) \begin{align}\hspace {-3 cm} u^b_n&=0, \end{align}

where $s$ and $b$ represents top and bottom surfaces. The normal to the top surface is

(3.24) \begin{align} \boldsymbol{n}^s&=\frac {\boldsymbol{\nabla }\!F}{|\boldsymbol{\nabla }\!F|}=\left [J^2+\displaystyle \left ({\mathcal{R}_\phi }\frac {\partial h}{\partial \theta }\right )^2\right ]^{-\frac {1}{2}}\left ({J\hat {\boldsymbol{g}}_1-{\mathcal{R}_\phi }\frac {\partial h}{\partial \theta }\hat {\boldsymbol{g}}_2}\right )\!. \end{align}

Employing (3.24), we obtain the non-dimensionalised dynamic BCs at the top surface as

(3.25a) \begin{align} \sigma ^s_{nn}-\varepsilon \frac { \sigma ^s_{n\theta }}{J^s\kappa ^s_{\phi }}\frac {\partial h}{\partial \theta }&=0, \\[-12pt]\nonumber \end{align}

(3.25b) \begin{align} \sigma ^s_{n\theta }-\varepsilon \frac { \sigma ^s_{\theta \theta }}{J^s\kappa ^s_{\phi }}\frac {\partial h}{\partial \theta }&=0 \end{align}

and

(3.25c) \begin{align} \qquad \sigma ^s_{n\phi }-\varepsilon \frac { \sigma ^s_{\theta \phi }}{J^s\kappa ^s_{\phi }}\frac {\partial h}{\partial \theta }&=0. \end{align}

At the bottom surface, with $i$ representing either $\theta$ or $\phi$ , we have

(3.26a) \begin{align} \sigma ^b_{ni} = -\mu \hat {u}^b_i\sigma ^b_{nn}, \end{align}

where, we have introduced the scaled velocity

(3.26b) \begin{align} \hat {u}^b_i =\frac {u^b_i}{\sqrt {(u^b_\theta )^2+(u^b_\phi )^2}}. \end{align}

3.3. Depth averaging

We now define the depth-averaged quantities, identified by $\bar {(\boldsymbol{\cdot })}$

(3.27a) \begin{align} \,\,\, \left \{\overline {v}(\theta ,t), \overline {w}(\theta ,t)\right \} &= \frac {1}{h}\int ^{h}_{0}\left \{ v(n,\theta ,t), w(n,\theta ,t)\right \} \, {\text{d}} n \end{align}

and

(3.27b) \begin{align} \hspace {-1 cm} \overline {vw} &= \beta _{vw}(\theta ,t) \overline {v} \,\overline {w}, \end{align}

where $h$ is the flow’s depth scaled by its initial average thickness $H$ and the depth-variation parameter $\beta _{vw}$ is introduced to express the average of the product of velocities as the product of their individual averages. A plug-like velocity profile is a fairly good assumption for granular flows (Melosh Reference Melosh1986; Pudasaini & Hutter Reference Pudasaini and Hutter2007). Consequently, we approximate $\beta _{vw}$ as unity, as did Gray et al. (Reference Gray, Wieland and Hutter1999). This allows us to equate the depth-averaged velocity to the basal velocity, i.e.

(3.28) \begin{equation} \bar {\boldsymbol u}= {\boldsymbol u}^b. \end{equation}

To close our mathematical description, we need a rheological description of the shallow regolith flow. We express the stress tensor as

(3.29) \begin{align} \boldsymbol{\sigma } =-p\boldsymbol{1} +\boldsymbol{\tau }, \end{align}

where $p$ is the pressure and $\boldsymbol \tau$ is the deviatoric stress tensor. Gray & Edwards (Reference Gray and Edwards2014) and Baker, Barker & Gray (Reference Baker, Barker and Gray2016) proposed a depth-averaged constitutive law of the form

(3.30) \begin{equation} h\bar {\boldsymbol{\tau }} =\rho \nu h^{3/2} \bar {\boldsymbol{D}}, \end{equation}

where $\nu$ and $\bar {\boldsymbol{D}}$ are the depth-averaged viscosity and two-dimensional strain rate tensor, respectively. The depth-averaged viscosity $\nu$ may be expressed in terms of material parameters and slope of the basal surface; see Baker et al. (Reference Baker, Barker and Gray2016, (2.27)). However, the contribution from $\bar {\boldsymbol \tau }$ in the depth-averaged linear momentum balance is generally rather small (Baker et al. Reference Baker, Barker and Gray2016; Rocha, Johnson & Gray Reference Rocha, Johnson and Gray2019). Given this, for simplicity, here we ignore $\bar {\boldsymbol \tau }$ and set

(3.31) \begin{align} \overline {\sigma }_{\theta \theta }=\overline {\sigma }_{\phi \phi } =\overline {\sigma }_{nn}=-\overline {p}. \end{align}

This is consistent with the observation (Melosh Reference Melosh1986) that downslope granular flows are plug-like, with shearing confined near the base, so that shear stresses are generally small through the flow’s depth. The only resistance to the flow is thus from basal friction.

We now list some further approximations that we make while depth averaging. The curvatures may be written as $\kappa _i = \kappa _i^b +O(\varepsilon )$ , because we assume curvatures to be of $O(1/R_{\textit{CB}})$ . Therefore, we set curvatures to be constant while depth averaging, which introduces an error of $O(\varepsilon )$ . We also neglect the contribution of Coriolis force arising from the normal component $u_n$ of the flow velocity, as it is also of $O(\varepsilon )$ . Considering only the partial contribution from the Coriolis term is termed as the traditional approximation in the geophysical literature (Eckart Reference Eckart1960; Gerkema et al. Reference Gerkema, Zimmerman, Maas and Van Haren2008). We introduce the above approximations in the sections below. However, for completeness, in Appendix B we provide equations for regolith flow that include the constitutive law (3.30), the effect of curvature variations and which do not make the traditional approximation.

The depth-averaged momentum equation in the $n$ -direction now becomes

(3.32a) \begin{align}\hspace {-1 cm} p &= \psi (h-n) +O(\varepsilon ), \end{align}

where

(3.32b) \begin{align} \hspace {3 cm}\,\,\, \psi &= -(f_n +\kappa _{\theta }{{\bar u}_\theta ^2} + \kappa _{\phi n}{\bar u}^2_{\phi })&& \end{align}

is defined in terms of the body force term

(3.32c) \begin{align} f_n&=b^0_n+\mathcal{R}_\phi {\kappa _{\phi n}}\big (\mathcal{R}_\phi {\omega ^2} +2\omega {\bar u}_\phi \big ). \end{align}

The form (3.32a ) and (3.32b ) of the pressure field is similar to the one obtained by employing the quasi-hydrostatic approximation in geophysical fluid dynamics (White et al. Reference White, Hoskins, Roulstone and Staniforth2005). In (3.32b ) and (3.32c ), we note the presence of contributions of centripetal acceleration through terms quadratic in velocity components $u_\theta$ and $u_\phi$ and angular velocity $\omega$ . Coriolis acceleration results in the last term in $f_n$ . While centripetal terms always lower the basal pressure $p$ , the Coriolis effect may raise or lower $p$ depending upon the signs of $\omega$ , the curvature $\kappa _{\phi n}$ and the azimuthal velocity ${\bar u}_{\phi }$ . As mentioned, mass shedding commences when basal pressure becomes negative which, from (3.32a ), is equivalent to the requirement that

(3.33) \begin{equation} \psi \lt 0. \end{equation}

Incorporating the above, our final depth-averaged continuity and linear momentum balance equations governing shallow regolith flow when mass shedding is absent, i.e. ${\psi }\gt 0$ , become, respectively,

(3.34) \begin{align} \frac {\partial }{\partial t}( J h) + \frac {\partial }{\partial \theta }\left (\mathcal{R}_\phi {{\bar u}_\theta } h\right ) &=0, \\[-12pt]\nonumber \end{align}

(3.35) \begin{align} \frac {\partial }{\partial t}( J {\bar u}_\theta h) + \frac {\partial }{\partial \theta }\left \{\left ({{\bar u}_\theta ^2} +\frac {\varepsilon }{2}\psi h\right )\mathcal{R}_\phi h\right \} &=\Bigg [ -{\mu \psi \hat {{\bar u}}_\theta } + b^0_\theta +\, \kappa _{ \phi g}\left \{{\bar U}^2_\phi + \varepsilon {\psi }\frac {h}{2}\right \} \Bigg ]\textit{Jh} \end{align}

and

(3.36) \begin{align} \hspace {-3 cm} \frac {\partial }{\partial t}\left ( {J \mathcal{R}_\phi {\bar U}_\phi } h \right ) + \frac {\partial }{\partial \theta }\big (\mathcal{R}^2_\phi {\bar u}_\theta {\bar U_\phi } h\big ) &= -{\mu J\mathcal{R}_\phi {\psi }{\hat {{\bar u}}}_\phi } h, \end{align}

where ${\hat {{\bar u}}}_\theta$ and ${\hat {{\bar u}}}_\phi$ are defined as in (3.26b ) and

(3.37) \begin{align} {\bar U}_\phi = {\bar u}_\phi +\mathcal{R}_\phi {\omega } \end{align}

is the regolith’s azimuthal velocity in the inertial frame. The equations above have to be solved for the flow rate $\left \{{\bar u}_\theta , {\bar u}_\phi \right \}$ and depth $h$ of the landslide. Finally, mass is shed from regions where ${\psi }\lt 0$ . In those regions, because we employ a depth-averaged formulation, mass shedding necessarily implies that material across the full flow depth is removed. Thus, we set

(3.38) \begin{equation} h={\bar u}_\theta ={\bar u}_\phi =0. \end{equation}

In Appendix D we provide details of how we have validated our landslide model above.

We now make two remarks regarding the regolith’s total energy and the angular momentum. First, the energy in this system is dissipated only due to basal sliding. Indeed, the energy equation is shown in Appendix C to be

(3.39) \begin{align} \frac {\partial }{\partial t}\left \{\textit{Jh}\left (\frac {|\bar {U}|^2}{2} +V(\theta ) - \varepsilon b^0_n\frac {h}{2}\right )\right \} + \frac {\partial }{\partial \theta }\left \{\left (\frac {|\bar {U}|^2}{2} +V(\theta ) - \varepsilon b^0_n h\right )\mathcal{R}_\phi {\bar u}_\theta h\right \} \nonumber \\ = -\mu J\psi h \big (|{\bar u}|+ {\mathcal{R}_\phi \omega \hat {{\bar u}}_\phi }\big ) + O(\varepsilon ), \end{align}

where $|\bar {U}|^2 = \sqrt {{\bar U}^2_\theta + {\bar U}^2_\phi }, \ \ |\bar {u}|^2 = \sqrt {{\bar u}^2_\theta + {\bar u}^2_\phi }$ and $V(\theta )$ is the gravitational potential at the surface the of central body, so that

(3.40) \begin{equation} b^0_\theta = - \frac {\mathcal{R}_\phi }{J}\frac {\partial V}{\partial \theta }. \end{equation}

In the next section we show that (3.36) represents the angular momentum balance of a ring-shaped volume element ${\text{d}} V$ of the axisymmetric regolith flow – sketched in brown in figure 2(a). Thus, (3.36) reveals that the angular momentum of this ring element is advected across latitudes by the velocity ${\bar u}_\theta$ and exchanged with the central body by basal friction.

3.4. Rotational dynamics

It remains to follow the evolution of the body’s rotation rate $\omega$ . The angular momentum balance (3.3) after non-dimensionalisation and invoking axisymmetry becomes

(3.41a) \begin{align} \hspace {2 cm}\,\,\,\, \alpha = \left \{T-{2\pi \varepsilon }\int _{\psi \geq 0} {J}\mathcal{R}_\phi \frac {\partial }{\partial t}\left (h{{\bar U}}_\phi \right ) \,{\text{d}}\theta \right \}I_c^{-1} +O(\varepsilon ^2), \end{align}

where

(3.41b) \begin{align} \hspace {-7 cm} I_c = {2\pi }\int ^\pi _0\frac {1}{5}\mathcal{R}_\phi ^3R^2\,{\text{d}} \theta \end{align}

is the central body’s moment of inertia about symmetry axis. The integral in (3.41a ) is over those regions where the basal pressure does not vanish, i.e. $\psi \geqslant 0$ . We find that the term inside the integral in (3.41a ) is exactly the first term in (3.36) which may, hence, be interpreted as the angular momentum balance for the ring-shaped granular element ${\text{d}} V$ in figure 2(a). We note that employing (3.36) we may also write the expression for angular acceleration in the simpler form

(3.42) \begin{equation} \alpha = I_c^{-1}\left \{T+{2\pi \varepsilon }\int _{\psi \geq 0} \mu J\mathcal{R}_\phi \psi \hat {\bar {u}}_\phi h \,{\text{d}}\theta \right \} +O(\varepsilon ^2). \end{equation}

From (3.41a ) or (3.42) we may determine the angular acceleration $\alpha$ at any time in terms of the average azimuthal velocity $\bar {u}_\phi$ and the depth $h$ of the landslide. Time integration of $\alpha$ yields $\omega$ as a function of time.

4. Applications

We now investigate two different applications in turn.

4.1. Isolated landslides

We investigate regolith flows on rotating bodies in two limits, viz. the rotation rate being zero and very fast. This is done to understand gravity- and rotation-driven flows, respectively.

Figure 3.

Final regolith deposited by landsliding of an initially uniform regolith layer on stationary and rotating oblate and prolate ellipsoids. The regolith height is magnified two times for clarity.

We initiate the flow by distributing regolith uniformly over the body’s surface with depth $H=\varepsilon R_{\textit{CB}}$ , where we arbitrarily set $\varepsilon = 0.01$ . Thus, the initial conditions are

(4.1) \begin{align} &\text{at}\quad t=0: \qquad h=1 \quad \text{and} \quad \bar {u}_\theta =\bar {u}_\phi =0 \quad \text{for all} \quad \theta \in (0,\pi ), & \end{align}

where $h$ is the flow’s scaled depth. The flow of the ensuing landslide is then found by numerically solving (3.34)–(3.36) utilising the non-oscillatory central schemes of Kurganov & Tadmor (Reference Kurganov and Tadmor2000). The gravity field of the axisymmetric body is found by slicing the body into infinitesimal discs along the symmetry axis and then adding their individual contributions (Huré Reference Huré2012).

Figure 3 presents the final deposition profile of an initially uniformly distributed regolith layer on an oblate and a prolate body that were either stationary or rotating fast. When rotation is absent, the gravitational field is too weak to overcome the basal resistance that keeps the regolith immobile; thus results are shown for a small friction angle $\zeta =10^\circ$ . We find that landslides transport grains to the poles on an oblate body and to the equator for a prolate body when they are not rotating. This is because the tangential component of the gravity is always towards the poles for an oblate body, but towards the equator for a prolate body. In contrast, rotational effects always drive grains to the equator, as the tangential component of the centrifugal force along the body’s surface is always directed towards the equator. Thus, for a prolate body, rotation assists the gravity-driven motion augmenting grain accumulation at the equator. On an oblate body, rotation and gravity compete, and we may find a comparable amount of regolith at both poles and the equator. We remark that rotation about the long axis of a prolate body is generally unstable. Thus, that example may be treated as illustrative.

We consider now applications to small planetary objects. For this we create axisymmetric bodies from shape models of actual, nominally axisymmetric, asteroids, such as Bennu and Ryugu from figure 1, and Itokawa. This is done by averaging the surface profiles of these bodies across longitudes at a fixed latitude.

Figure 4 shows final deposits of initially uniformly spread regolith on asteroids Bennu, Ryugu and Itokawa. Were Bennu and Ryugu not rotating, then we would find grains accumulated at mid-latitudes. This is because the equilibria at the surface of such top-shaped bodies are located at their mid-latitudes (Banik et al. Reference Banik, Gaurav and Sharma2022). Interestingly, we find grains deposited around the poles of Ryugu due to the presence of a local topographical depression. These are the remnants of the initially placed regolith layer. Because the tangential component of both centrifugal and gravitational force is weak at the poles, the grains in their vicinity did not move away. New grains are not transported to the poles. At a higher rotation rate, we find that regolith accumulates at the equator, as in the examples in figure 3.

Finally, while Itokawa may be somewhat crudely approximated as an axisymmetric body, it is observed to rotate very slowly with a time period of 12 h about its short axis. Thus, to fit into our modelling we here investigate landsliding on a non-rotating Itokawa. We find that, while regolith migrates towards Itokawa’s poles, a significant fraction is captured in the depression around the ‘neck’ region where Itokawa’s ‘head’ and ‘body’ are joined.

Figure 4.

Final regolith deposited by landsliding of an initially uniform regolith layer on axisymmetrically modelled asteroids Bennu, Ryugu and Itokawa. The regolith thickness is magnified two times for clarity.

4.2. Shape and spin evolution of rubble piles

We now explore two instances of the evolution of the shape and spin state of an initially spherical rubble asteroid over its lifetime due to impact-induced landsliding and the action of YORP torque. Before doing that we discuss in some detail the parameters we employ in our simulation, the manner in which they are chosen and how initial conditions are specified.

4.2.1. Simulation details

As already mentioned in the introduction and at the end of § 2, we do not model the impact process here. We simply assume that over time the asteroid experiences a fraction of impacts that are energetic enough to bring about global landsliding. Given the time $t_{\!s}$ over which we wish to investigate the asteroid’s evolution, we follow the procedure of Banik et al. (Reference Banik, Gaurav and Sharma2022) to calculate the number of such energetic impact events. We first estimate the minimum size $d_{\textit{min} }$ of the impactor that may cause global landslides (Richardson et al. Reference Richardson, Melosh, Greenberg and O’Brien2005). We then compute the number of impactors with size bigger than $d_{\textit{min} }$ during time $t_{\!s}$ (Holsapple Reference Holsapple2022). These impacts may occur stochastically over a time scale $t_{\!s}$ , but here we report results corresponding to a uniform temporal distribution. We then integrate the equations of motion, employing the parameters listed in table 1. Whenever an impact happens, we initiate a landslide by instantaneously mobilising a constant thickness of layer spread over the body. The thickness of the mobilised regolith is set in the manner of Banik et al. (Reference Banik, Gaurav and Sharma2022), who took a thickness-to-slant-height ratio of $\varepsilon _{ls}=0.005$ for a top-shaped asteroid. For our case, this corresponds to a thickness-to-radius ratio of $\varepsilon =0.0075$ . The initial condition for the regolith flow are thus given by

(4.2) \begin{align} &\text{at}\quad t=0: \qquad h=1, \quad \bar {u}_\theta =\bar {u}_\phi =0 \quad \text{for all} \quad \theta \in (0,\pi ) \quad \text{and} \quad \omega = \omega _0, \end{align}

where $\omega _0$ is the primordial rotation rate of the body.

Table 1.

Parameters used for the simulations.

When an impact mobilises a landslide, we lower the basal friction angle to incorporate the effects of post-impact seismic reverberations. This is similar to acoustic fluidisation in large landslides (Melosh Reference Melosh1986). We set the basal friction angle to $10^{\circ }$ during the time over which seismic shaking persists, following the Luca et al. (Reference Luca, Tai and Kuo2016) estimate. The basal frictional angle recovers to its original value of $30^\circ$ once seismic shaking ceases; note that $30^\circ$ is typical of geological materials in terrestrial conditions (Nedderman Reference Nedderman2005). The duration for seismic shaking on the Moon was observed to be around 1 h. by the Apollo missions (Latham et al. Reference Latham, Ewing, Press, Sutton, Dorman, Nakamura, Toksöz, Wiggins, Derr and Duennebier1970). In the absence of more precise constraints, we allow seismic effects to persist for this length of time in our simulations, which corresponds to $t\approx 2$ non-dimensional time units.

Once the landslide stops, grains are deposited on the body changing its overall shape. Subsequently, any new landslides mobilises a uniform regolith layer over this new shape. The body’s rotational dynamics also incorporates its modified shape. We correspondingly update the Jacobian, curvatures and gravity for the body. All through, YORP torque continues to act. We take the YORP torque to be equal to that experienced by the asteroid Bennu were it is located in the main belt. To achieve this, the YORP torque was scaled by multiplying the torque experienced by Bennu (Hergenrother et al. Reference Hergenrother, Maleszewski, Nolan, Li, Drouet d’Aubigny, Howell, Kareta, Izawa and Barucci2019) by $(a_0/a)^2$ , where $a_0=1.1$ AU is Bennu’s original semi-major axis and $a=2.3$ AU corresponds to that of an object located in the main belt (Bottke et al. Reference Bottke, Vokrouhlickỳ, Walsh, Delbo, Michel, Lauretta, Campins, Connolly, Scheeres and Chelsey2015).

Before closing this section, we make the following observation. The YORP torque is very sensitive to variations in surface features (Statler Reference Statler2009). Consequently, we should ideally account for spin evolution driven by YORP torque along with the modification of this torque due to reshaping, as done in the DE simulations of Cotto-Figueroa et al. (Reference Cotto-Figueroa, Statler, Richardson and Tanga2015). Here, however, we focus on landslide-induced shape evolution and treat the YORP torque as constant. This is done to understand the importance of landslide-driven processes. Fully coupled models will be developed in the future.

4.2.2. Results and discussion

Figure 5 shows the evolution of the rotation rate and shape of an initially spherical asteroid over 4 Myr under the combined effects of a constant YORP torque and impact-induced global landsliding. Following the methodology described above, we estimate that the asteroid experiences 58 impacts large enough to initiate global landslides during this period. As noted earlier, these impacts are distributed uniformly over 4 Myr, and their occurrences are indicated by the spikes in the red curve in figure 5. The asteroid evolves into a top-shaped configuration, similar to the bodies shown in figure 1. Its shape is progressively modified by the deposition of multiple landslides, as illustrated schematically in figure 2(b).

Figure 5.

The evolution of the shape and rotation rate of an initially spherical rubble asteroid due to a constant YORP torque and impact-induced landslides over 4.0 Myr. The blue dashed curve tracks changes in the body’s rotation rate when only the YORP torque is active.

The asteroid’s rotation rate in figure 5 increases continuously under the YORP torque, while each landsliding event reduces this rate by redistributing regolith toward the equator, thereby increasing the body’s moment of inertia. This feedback prevents the rotation rate from exceeding a critical limit. In contrast, the rotation rate under the action of the YORP torque alone – shown by the blue dashed curve in figure 5 – rises steadily. Along this curve the body’s initial spherical shape is kept fixed. This continuous increase in the rotation rate will eventually lead to rotational fission, i.e. surface mass loss driven by centrifugal forces, when the rotation period drops below 3 h for our assumed density of 1250 $\textrm {kg}/\textrm {m}^3$ (Scheeres Reference Scheeres2015). However, when landsliding effects are included, rotational fission is suppressed; indeed, the rotation period along the red curve in figure 5 remains above 3.75 h. This aspect is further probed in figure 7. We recall from (3.32a ) that the basal pressure $p=\psi h$ must vanish for mass shedding. This is equivalent to the requirement (3.33) that $\psi \lt 0$ . Figure 7 plots $\psi$ at the equator at times when landsliding is absent, so that $\bar {\boldsymbol u} = {\boldsymbol 0}$ in (3.32b ). Under these conditions $\psi$ functions as an effective gravity, comprising gravitational and centrifugal accelerations. We find that $\psi$ at the equator remains positive during the body’s evolution, confirming that mass shedding does not take place.

Figure 6.

Same as figure 5 but the initial rotation period is decreased to 3.6 h. The dashed curve shows the spin evolution of figure 5.

Figure 7.

The change in $\psi$ given by (3.32b ) at the equator during the evolution of initially slow and fast rotating rubble asteroids reported in, respectively, figures 5 and 6.

This result is significant, suggesting that rotational fission may be far less common – if it occurs at all – during the evolution of small rubble-pile bodies than is currently assumed. Finally, figure 5 shows that the asteroid’s top-shaped form develops within a few Myr, which may account for the observed abundance of top-shaped asteroids in near-Earth orbits.

Figure 6 reproduces figure 5, except that the primordial body is initially rotating faster. For comparison, the rotation rate curve from the simulation in figure 5 is also included. We now find that the decrease in rotation rate caused by landsliding events dominates over the increase due to the YORP torque, resulting in a net lengthening of the rotation period. From figure 6, we infer that the shape changes induced by landslides – and hence their influence on the body’s rotational dynamics – are more pronounced at higher rotation rates. Furthermore, figure 7 shows that, although the spin rate decreases in figure 6, $\psi$ at the equator continues to decline. This behaviour is attributed to changes in the gravitational field resulting from the progressive elongation of the body’s shape.

Comparing figures 5 and 6, we note that the final shapes are qualitatively similar. Moreover, the rotation rate curves in figure 6 appear to converge to nearly the same final value. We attribute this to the use of the same uniform distribution of impacts, so that landslides occur at identical times during the body’s evolution. The outcome would differ if the impacts occurred randomly.

The evolution of $\psi$ – the effective gravity – at the equator in figure 7 suggests that a continued increase in the angular velocity driven by the YORP effect and/or further bulging at the equator due to landsliding may eventually result in mass shedding from the equatorial region. Allowing for landsliding and the associated reshaping delays this mass shedding, which may be identified as ‘rotational fission’. At the same time, we note that our model assumes that post-impact regolith mobilisation occurs uniformly over the surface. In reality, impacts are expected to preferentially mobilise larger volumes of regolith at higher elevations, which would act to suppress the equatorial bulge and, thereby, further delay the onset of mass shedding. We defer a quantitative investigation of this effect to future work.

Two technical points are worth noting. First, in the simulations of figures 5 and 6, the curvature ceases to remain small after a few Myr, violating the approximation in § 3.3 that curvatures are small. To address this, we employed the more accurate – and more complex – formulation given by (B6)–(B7), which retains all terms at $O(\varepsilon )$ . Second, the simulations were terminated once the curvature near the equator became excessively large and mass shedding began.

We now connect our results with past work, focusing on studies that investigated surface mass motion. Harris et al. (Reference Harris, Fahnestock and Pravec2009) examined shape evolution driven by the movement of surface grains by considering the equilibrium of a single particle within the effective geopotential of a spinning body and its relocation once a local slope-stability criterion was violated. From this analysis, they derived axisymmetric equilibrium shapes that were successfully fitted to top-shaped asteroids. This work was extended by Scheeres (Reference Scheeres2015) to include spin changes due to the YORP torque and a more detailed treatment of post-failure regolith redistribution. It is interesting that some of the shapes observed during the evolution in figures 5 and 6 are similar to the ones reported in Scheeres (Reference Scheeres2015) figures 9, 10, and this lends support to both our approaches. Hirabayashi (Reference Hirabayashi2015), Hirabayashi et al. (Reference Hirabayashi, Sánchez and Scheeres2015) and Sánchez & Scheeres (Reference Sánchez and Scheeres2020) further developed this line of investigation to include cohesive grains, employing a combination of linear elasticity, and FE and DE simulations. An extension to irregular shapes was presented by Hirabayashi & Scheeres (Reference Hirabayashi and Scheeres2019). Collectively, these studies predicted the onset of failure due to increasing rotation rates from YORP torques and revealed the possibility of bulk failure preceding surface mobilisation when the interior is weak.

The above studies were essentially equilibrium analyses: they predicted the onset of failure, but did not follow the body’s dynamical evolution after failure. Because of computational limitations, perhaps, they did not follow the long-term evolution of a rubble body’s shape and spin state and, hence, could not comment on the time scales over which the body may achieve its final shape. They also focused on regolith flow triggered by centrifugal effects alone which, for a body with a strong interior, would mobilise only a small amount of surface material. In contrast, processes such as impact cratering may mobilise much larger volumes of regolith flows that, as shown here, will exert a significant influence on the body’s rotational dynamics – to the extent that rotation rates may be suppressed at values well below those at which centrifugal effects can drive landslides.