An expression for the angle of repose of dry cohesive granular materials on Earth and in planetary environments

Fig. 1 displays snapshots of simulations obtained with $d=30\hspace{0.17em}\mathrm{\mu }$m and $d=1$ mm under terrestrial gravity conditions. To measure the angle of repose, the projection of the heap on two orthogonal planes ($x$ and $y$) is computed, which yields two side-view images of the heap separated by a rotation of $90°$. An image-processing algorithm is then applied to compute the upper boundary of the granular material, i.e., the heap surface, as well as the total surface area associated with the two-dimensional heap projection. Following ref. 18, we define the angle of repose as the base angle of the isosceles triangle, which has the same surface area as the heap (Fig. 1). Therefore, two values of this angle are obtained for each heap, from which the mean ${\theta }_{\mathrm{r}}$ and SD ${\sigma }_{{\theta }_{\text{r}}}$ are computed. We have found no considerable changes in the angle of repose by including computations of this angle from additional projections separated by smaller rotation angles (18).

Moreover, we have found that heaps formed by particles of diameter smaller than a value $d_{\mathrm{c}}$ ($\approx 50\mathrm{\mu }$m under terrestrial gravity) display a strongly irregular shape, with considerable variations in local slope throughout their surface (Fig. 1A). In this regime of strongly irregular heap morphology, which we call here regime I, particle avalanches are largely inhibited by the formation of large, irregular particle agglomerates along the surface of the heap. By contrast, heaps formed with particles of size $d\gtrsim d_{\mathrm{c}}$ (regime II) display rather regular morphology and nearly constant slope along their surface (Fig. 1B).

We remark that it becomes increasingly difficult to produce heaps in our simulations as $d$ decreases down to values smaller than $d_{\mathrm{c}}$. In particular, we could not obtain any heap for $d$ smaller than about $50\%$ of $d_{\mathrm{c}}$, since particles’ sticky behavior for such small particle sizes causes the heap to reach the upper limit of the simulation domain. The conclusions presented in our article are thus applicable to heaps in regime II, i.e., provided particles are not smaller than about $d_{\mathrm{c}}$.

The solid squares in Fig. 2 denote the values of ${\theta }_{\mathrm{r}}\left(d\right)$ obtained from DEM simulations under Earth’s gravity. Moreover, Fig. 2 also displays a comprehensive compilation of experimental data for the angle of repose as a function of particle size, taken from measurements of ${\theta }_{\mathrm{r}}$ from granular heaps produced with different materials (see legend and key in Fig. 2). As we can see in Fig. 2, our numerical predictions follow remarkably well the trend of the experimental data.

Furthermore, Fig. 3A shows the values of ${\theta }_{\mathrm{r}}$ obtained for different values of gravity, ranging from Pluto gravity ($g=0.62$ m$/{\mathrm{s}}^2$) to 100 times Earth’s gravity. We see that the same qualitative behavior of ${\theta }_{\mathrm{r}}$ as a function of $d$ is found for all values of gravity investigated. Moreover, for a given particle size, a trend of increasing ${\theta }_{\mathrm{r}}$ with decreasing gravity is observed. In the following section we present a mathematical expression that will shed light on our numerical predictions, as well as on the experimental observations of the angle of repose as a function of particle size.

To obtain a mathematical expression that fits our DEM predictions for ${\theta }_{\mathrm{r}}$ as a function of particle size, we first note that the angle of repose is described via ${\theta }_{\text{r}}\approx {\tan }^{-1}\left[{\mu }_{\mathrm{e}\mathrm{ff}}\right]$, where ${\mu }_{\mathrm{e}\mathrm{ff}}$ is an effective coefficient of friction describing the maximal friction resistance of the granular material against sliding through avalanche (23). This coefficient encodes information on material properties, as well as on macroscopic characteristics related to the heap structure, such as packing fraction and surface roughness.

However, these characteristics are affected by the nature of interparticle interaction forces. As particle size decreases, attractive particle interaction forces become increasingly more relevant compared to gravitational forces, thereby leading to an increase in the bulk porosity and affecting the dynamics of stress distribution within the heap. Experiments and DEM simulations revealed that the solid fraction of powders follows, approximately, a power-law relation with particle size (15). Taking inspiration from this insight, we approximately describe ${\theta }_{\mathrm{r}}$ (in radians) using the expressionwhere the effective friction coefficient ${\mu }_{\mathrm{e}\mathrm{ff},\infty }$ and the characteristic length-scale $D$ must be determined from the best fit to observations of the angle of repose as a function of particle size. The best fit to our simulation data (blue solid squares in Fig. 2) using Eq. 1 gives ${\mu }_{\mathrm{e}\mathrm{ff},\infty }\approx 0.447$ and $D\approx 87\hspace{0.17em}\mathrm{\mu }$m, with coefficient of correlation $R^2\approx 99\%$. This best fit is denoted by the solid line in Fig. 2.

Moreover, the solid lines in Fig. 3A denote the best fits to the DEM results associated with various values of g level; i.e., $\tilde{g}=g/g_{\mathrm{E}\mathrm{a}\mathrm{r}\mathrm{t}\mathrm{h}}$, where $g_{\mathrm{E}\mathrm{a}\mathrm{r}\mathrm{t}\mathrm{h}}=9.81$ m$/$${\text{s}}^2$, using Eq. 1 and assuming that ${\mu }_{\mathrm{e}\mathrm{ff},\infty }\approx 0.447$ is independent of gravity. The dotted line in Fig. 3A corresponds to the asymptotic value of the angle of repose, ${\theta }_{\mathrm{r},\infty }\equiv {\tan }^{-1}\left[{\mu }_{\mathrm{e}\mathrm{ff},\infty }\right]$. Indeed, we find that the angle of repose ${\theta }_{\mathrm{r},\mathrm{n}\mathrm{c}}$ obtained from simulations without attractive particle interaction forces (i.e., assuming noncohesive materials) is approximately equal to ${\theta }_{\mathrm{r},\infty }$, independent of particle size (Fig. 3B). There is a small trend of increasing ${\theta }_{\mathrm{r},\mathrm{n}\mathrm{c}}$ with gravity—Fig. 3B suggests an increase of the order of $10\%$ in ${\theta }_{\mathrm{r},\mathrm{n}\mathrm{c}}$ with an increase in $\tilde{g}$ from 0.06 to 100—which we attribute to a slightly greater tendency of particles to interlock upon an increase in gravitational stress, but further work (not in the scope of this article) would be needed to elucidate this small effect. Fig. 3 A and B thus shows that particle size controls on ${\theta }_{\mathrm{r}}$ are dictated mainly by cohesive particle–particle interactions.

Fig. 3C displays the values of $D$ obtained from the best fits in Fig. 3A, as a function of the $g$ level. As we can see from Fig. 3C, for the range of $\tilde{g}$ considered in the present work, the characteristic length-scale $D$ can be approximately described by the equationwhere $D_{\mathrm{E}}$ denotes the terrestrial value of $D$, i.e., $87\mathrm{\mu }$m. We remark that the solid red line in Fig. 3C is no fit to the simulation data (solid circles). Instead, the red line denotes Eq. 2, i.e., our conjecture that $D$ is equal to $D_{\mathrm{E}}/\sqrt{\tilde{g}}$. We have found that fitting the DEM results using the equation $D=D_{\mathrm{0}}/{\tilde{g}}^{\alpha }$ leads to $D_0\approx 76\hspace{0.17em}\mathrm{\mu }$m and $\alpha \approx 0.45$ with coefficient of correlation larger than $99\%$ (dashed line in Fig. 3C), which provides, thus, an exponent close to 0.5 as assumed in Eq. 2.

Using Eq. 2 we can rewrite Eq. 1 aswhich provides a means to estimate ${\theta }_{\mathrm{r}}$ (in radians) for given $g$ and $d$, provided $D_{\mathrm{E}}$ and ${\mu }_{\mathrm{e}\mathrm{ff},\infty }$ for the granular material are known, e.g., from experiments under terrestrial gravity.

From Eq. 3, we see that, in the limit $\sqrt{\tilde{g}}d/D_{\mathrm{E}}\to \infty$, the angle of repose asymptotically approaches the value ${\theta }_{\mathrm{r},\infty }\equiv {\tan }^{-1}\left[{\mu }_{\mathrm{e}\mathrm{ff},\infty }\right]$, independent of gravity and particle size. For the particulate material investigated here, ${\theta }_{\mathrm{r},\infty }\approx 24°$ (Figs. 2 and 3). Therefore, an approximate value for ${\mu }_{\mathrm{e}\mathrm{ff},\infty }$ can be estimated from experiments using particles that are so large that cohesive effects on the angle of repose can be neglected, whereupon ${\mu }_{\mathrm{e}\mathrm{ff},\infty }$ can be computed using the equation ${\mu }_{\mathrm{e}\mathrm{ff},\infty }=\tan \left({\theta }_{\mathrm{r},\infty }\right)$. Subsequently, experiments under Earth’s gravity using different values of $d$ can be performed to estimate $D_{\mathrm{E}}$ from the best fit to Eq. 3.

To shed light on Eq. 3, we note that granular materials display increasing relative cohesiveness the smaller the particle size is (15). We follow refs. 24 and 25 and characterize this cohesiveness by computing the Bond number (Bo), i.e., the ratio of the maximum cohesive force acting on a particle to the gravitational force. For the system investigated here, we obtain $\mathrm{B}\mathrm{o}=12\gamma /\left(d^{2}g{\rho }_{\mathrm{p}}\right)$ (using Eqs. 13 and 14), where ${\rho }_{\mathrm{p}}$ is the material density and $\gamma$ is the surface energy density, which scales the attractive interparticle (van der Waals) force. Given this expression for the Bond number and the scaling of the angle of repose with the particle size in Eq. 3, we propose the modelwhere ${\mu }_{\mathrm{e}\mathrm{ff},\infty }$ and $\beta$ are empirical parameters. As shown in Fig. 4, using this model, we find an excellent data collapse of our simulation results for all values of particle size and gravity considered, with ${\mu }_{\mathrm{e}\mathrm{ff},\infty }\approx 0.458$ and $\beta \approx 0.014$. Furthermore, from the definition of Bo, Eqs. 1 and 4 lead towhich relates the characteristic length-scale $D$ in Eq. 1 to the material parameters.

We note that our model is valid in the regime of regular heap morphology (regime II), i.e., when $d\gtrsim d_{\mathrm{c}}$. In particular, Eq. 1 predicts ${\theta }_{\mathrm{r}}\to 90°$ when $d\to 0$, but experiments revealed a nearly constant angle of repose for particle sizes of the order of a few micrometers, i.e., in regime I of irregular heap morphology (26). Our simulations show that $D\approx 2d_{\mathrm{c}}$ and that the heap morphology changes from regular (regime II) to irregular (regime I) when the Bond number exceeds approximately $10^4$, thus providing a means to approximately estimate $d_{\mathrm{c}}$ from the material properties.

We have also performed some simulations with heaps $150\%$ larger than the ones considered in the present discussion, but we could not find a significant change in the values of ${\theta }_{\mathrm{r}}$. We believe that, qualitatively, even larger heaps should display the same dependence of ${\theta }_{\mathrm{r}}$ on $d$ observed in the experiments quoted in Fig. 2 and in the DEM simulations in this article. However, it is uncertain to which extent our model parameters should be adapted to assess the angle of repose of much larger heaps quantitatively. In particular, the dynamics of much larger systems could be affected by particle–air interactions, which must be incorporated in a future extension of our model. Furthermore, since the packing dynamics are affected by the preparation history of the packing (see the discussion in ref. 15), the pouring height and the size of the insertion region chosen for reproducing experimental conditions could exert an important control on the value of ${\theta }_{\mathrm{r}}$. We have indeed noticed a trend of increasing avalanche sizes by increasing the width of the particle insertion region in our simulations of larger heaps. Therefore, the quantitative behavior of ${\theta }_{\mathrm{r}}$ for much larger heaps will constitute topic of future research.

Reduced-gravity experiments of granular heap formation (28) in the Bremen drop tower, using basalt beads of diameter $500\mathrm{\mu }$m, revealed a slight trend of increasing angle of repose with decreasing gravity for $0.03\le \tilde{g}\le 1$. Although the experiments were performed using a Hele–Shaw cell to obtain quasi-2D heaps (28), qualitatively, the behavior of ${\theta }_{\mathrm{r}}$ with gravity for these heaps is consistent with the prediction from our model that granular materials behave more cohesively and produce larger angles of repose under lower gravity.

Furthermore, parabolic flight experiments with slowly rotating drums, operating in the discrete avalanche regime at $g$ levels ranging from 0.1 to 1, revealed a slight trend of increasing angle of repose ${\theta }_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}}$—the slope above which the granular surface becomes unstable and starts to avalanche (29)—with decreasing gravity (30, 31). The increase of ${\theta }_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}}$ with decreasing $g$ observed in ref. 30 suggests enhanced cohesive behavior of granular materials under reduced gravity (30). Microgravity experiments with rotating drums operating at $g$ levels varying from 0.01 to 1.9 and at high rotational speeds, in the so-called continuous avalanching regime—in which the granular surface undergoes avalanche flow permanently—also revealed increasing cohesive behavior with decreasing $g$ (32, 33). In these experiments, a scaling of the slope angle of the granular surface with $\sqrt{1/g}$ was found (32, 33). To compare, experiments in the centrifuge at $g$ levels from 1 to 25 revealed insignificant correlation between this angle and gravity (34), as also predicted from our model.

We note that particle interaction forces can be affected by still poorly understood geochemical processes (35), while granular flows under low atmospheric pressure conditions can be further influenced by the presence of interstitial gas. Moreover, for values of Bond number smaller than unity, irregular particle geometries, inherent to natural granular systems, can cause statistical fluctuations comparable to the difference between the angles of repose predicted for spherical particles under $g$ levels of the same order. For instance, parabolic flight experiments using Toyoura sand ($d\approx 250\hspace{0.17em}\mathrm{\mu }$m) revealed no significant dependence of the angle of repose on gravity for the range of investigated $g$ levels ($1/6$ to 2) (36). We conclude, based on our results, that the angle of repose of natural sand under Martian gravity should not differ much from its terrestrial counterpart ($30°\text{to}35°$), which is consistent with observations from the slope of Martian dunes at their avalanche side (37).

Characterization of average particle size of strongly polydisperse systems can be made using different parameters, such as the Sauter mean particle size $d_{\mathrm{S}}$, i.e., the diameter of the particle that has the same volume to surface area ratio as the entire ensemble. Two comprehensive compilations of experimental datasets on ${\theta }_{\mathrm{r}}$ as a function of $d_{\mathrm{S}}$, presented in refs. 26 and 38, reveal the same characteristics as in Fig. 2, i.e., a noticeable, power-law–like decrease with $d_{\mathrm{S}}$ within an intermediate range of particle size, followed by slower decay in the limit of larger $d_{\mathrm{S}}$, with ${\theta }_{\mathrm{r}}$ asymptotically approaching a value ${\theta }_{\mathrm{r},\infty }$ between $20°$ and $30°$ that is roughly independent of particle size. We thus propose that Eq. 1 should be valid without regard of the parameter employed to characterize particle size, since the corresponding value of $D$ can be estimated from the best fit to ${\theta }_{\mathrm{r}}\left(d\right)$ using Eq. 1.

We further remark that Eq. 1 should break down when $d$ becomes much smaller than a certain value ($d_{\mathrm{c}}$), which, for the system considered in our DEM simulations, is of the order of $50\hspace{0.17em}\mathrm{\mu }$m. Indeed, Eq. 1 predicts that ${\theta }_{\mathrm{r}}\to 90°$ as $d\to 0$. However, measurements of ${\theta }_{\mathrm{r}}$ of strongly polydisperse dust systems (using the Sauter mean diameter $d_{\mathrm{S}}$ to characterize particle size) revealed a breakdown in the trend of increasing ${\theta }_{\mathrm{r}}$ with decreasing particle size, when particles become sufficiently small (26). As explained by Lanzerstorfer (26), the reason for this behavior is the formation of small agglomerates of dust-sized particles, which roll down the heap slope, thus limiting the angle of repose. An upper limit for ${\theta }_{\mathrm{r}}$ of about $55°$ was reported in the experiments in ref. 26.

Although further work will be required to characterize the dynamics of particle agglomerates in DEM simulations and to investigate the factors controlling the value of $d_{\mathrm{c}}$, we note that the highest values of ${\theta }_{\mathrm{r}}$ in most experimental datasets quoted in Fig. 2 are consistent with the upper bound reported by Lanzerstorfer (26). The only exception is the dataset of Lumay et al. (18) (open squares in Fig. 2), which displays ${\theta }_{\mathrm{r}}$ exceeding $70°$. A possible explanation for the higher ${\theta }_{\mathrm{r}}$ values in this dataset is that the silicon carbide abrasives adopted in ref. 18 consist of a broad distribution of complex particle geometric shapes, including elongated and sharp-edged particles. DEM simulations (39) have shown that cohesive rods form packings of lower solid fraction than that of cohesive beads; i.e., complex particle shapes can enhance the cohesive behavior of granular materials. By virtue of these observations and in view of our Eqs. 4 and 7, we thus expect that elongated particles, compared to beads, produce steeper heaps, while sharp edges could further enhance friction and heap stability. The present work could be thus continued by elucidating this behavior in DEM simulations.

However, future work should now focus on deepening the understanding of our findings from a theoretical analysis of granular heap stability, which poses a long-standing and challenging research topic in the field of granular matter physics (40).

The equations for the translational and rotational motion of the $i$th particle in the system readwhere ${\mathbf{r}}_i$ and ${\mathit{\omega }}_i$ denote the position and angular velocity of particle $i$, respectively; $\mathbf{g}$ is gravity; while ${\mathbf{F}}_{ij}$ and ${\mathbf{M}}_{ij}$ denote the interaction force and torque, respectively, exerted by particle $j$ on particle $i$. Furthermore, $n$ is the total number of particles in the system, each particle having diameter $d$, density ${\rho }_{\mathrm{p}}$, mass $m=\pi {\rho }_{\mathrm{p}}{d}^3/6$, and moment of inertia $I=m{d}^2/10$. In Eq. 8, ${\mathbf{F}}_{ij}$ denotes the sum of the contact force ${\mathbf{F}}_{ij}^{\mathrm{c}}$ and the nonbonded attractive interaction force ${\mathbf{F}}_{ij}^{\mathrm{a}}$ exerted by particle $j$ on particle $i$. These forces are modeled as discussed next.

For the normal component of the contact force, ${\mathbf{F}}_{ij,\mathrm{n}}^{\mathrm{c}}$, viscoelastic interaction is assumed (11), while a modified Cundall–Strack model (10) is applied to compute the tangential force, ${\mathbf{F}}_{ij,\mathrm{t}}^{\mathrm{c}}$. The equations for ${\mathbf{F}}_{ij,\mathrm{n}}^{\mathrm{c}}$ and ${\mathbf{F}}_{ij,\mathrm{t}}^{\mathrm{c}}$ readwhere ${\delta }_{ij,\mathrm{n}}=d-\left|{\mathbf{r}}_j-{\mathbf{r}}_i\right|$; ${\mathbf{e}}_{ij,\mathrm{n}}=\frac{{\mathbf{r}}_j-{\mathbf{r}}_i}{\left|{\mathbf{r}}_j-{\mathbf{r}}_i\right|}$; ${\mathbf{r}}_i$ and ${\mathbf{r}}_j$ are the positions of particles $i$ and $j$, respectively; $v_{ij,\mathrm{t}}=\left|{\mathbf{v}}_{ij,\mathrm{t}}\right|$, where ${\mathbf{v}}_{ij,\mathrm{t}}={\mathbf{v}}_{ij}-\left({\mathbf{v}}_{ij}\cdot {\mathbf{e}}_{ij,\mathrm{n}}\right)\cdot {\mathbf{e}}_{ij,\mathrm{n}}$ denotes the relative tangential velocity at the point of contact, with ${\mathbf{v}}_{ij}={\dot{\mathbf{r}}}_j-{\dot{\mathbf{r}}}_i$; and ${\mathbf{e}}_{ij,\mathrm{t}}={\mathbf{v}}_{ij,\mathrm{t}}/v_{ij,\mathrm{t}}$ is the tangential unit vector; while the contact forces in Eqs. 10 and 11 are set as zero if ${\delta }_{ij,\mathrm{n}}<0$. Furthermore, $k_{\mathrm{n}}$ and $k_{\mathrm{t}}$ are elastic parameters, which are computed from the particle Young’s modulus, $Y$, and the Poisson’s ratio, $\nu$, throughwith $G=Y/\left(2+2\nu \right)$ denoting the shear modulus and $R_{\mathrm{e}\mathrm{ff}}=d/4$, while $A_{\mathrm{n}}$ and $A_{\mathrm{t}}$ are dissipative parameters, ${\mu }_{\mathrm{s}}$ is the (Coulomb) sliding friction coefficient, and the integral in Eq. 11 is performed over the relative displacement of the particles at the point of contact and over the entire contact duration (15, 41).

The dissipative parameter $A_{\mathrm{n}}$ further depends on material viscosities that are not directly available (11). Therefore, to determine $A_{\mathrm{n}}$, we use a relation between this dissipative parameter and the coefficient of restitution for the collision of two isolated particles (42, 43) and employ the Padé approximation from ref. 44. Furthermore, following ref. 45, we choose the tangential dissipative parameter, $A_{\mathrm{t}}$, such that the prefactors of ${\dot{\delta }}_{ij,\mathrm{n}}$ and $v_{ij,\mathrm{t}}$ in Eqs. 10 and 11, respectively, are of the same order of magnitude, thus yielding $A_{\mathrm{t}}=3A_{\mathrm{n}}{k}_{\mathrm{n}}/2$. Indeed, ref. 46 showed that this assumption leads to excellent agreement between simulation results and experimental values of particle velocity profiles in a gravity-driven shearing experiment.

Furthermore, the nonbonded attractive van der Waals interaction force between two particles in the system, labeled $i$ and $j$, is given by (15, 47, 48)where $A_{\mathrm{H}}$ is the Hamaker constant, which is computed from the surface energy density, $\gamma$, using the expression (49)

In Eq. 13, $D_{max}$ is the maximal (cutoff) distance of the van der Waals interaction, which is set as $1\hspace{0.17em}\mathrm{\mu }$m (15), while the constant $D_{min}$ is introduced to avoid the singularity of the Hamaker equation at ${\delta }_{ij,\mathrm{n}}=0$ and accounts for the fact that particles are not perfectly smooth. We take the value $D_{min}=1.65$ Å, which corresponds to the minimal distance between particles at contact (49, 50).

The total torque ${\mathbf{M}}_{ij}$ is given by ${\mathbf{M}}_{ij}={\mathbf{M}}_{ij}^{\mathrm{t}}+{\mathbf{M}}_{ij}^{\mathrm{r}}$, where ${\mathbf{M}}_{ij}^{\mathrm{t}}=\mathrm{\Theta }\left({\delta }_{ij,\mathrm{n}}\right)\left[d/2-{\delta }_{ij,\mathrm{n}}\right]{\mathbf{e}}_{ij,\mathrm{n}}\times {\mathbf{F}}_{ij,\mathrm{t}}^{\mathrm{c}}$, and ${\mathbf{M}}_{ij}^{\mathrm{r}}$ is the torque associated with the rolling resistance at the contact point between particles $i$ and $j$ and is computed through the model (14, 17)where $M_{ij}^{\mathrm{r}\mathrm{m}}={\mu }_{\mathrm{r}}{R}_{\mathrm{e}\mathrm{ff}}\left|{\mathbf{F}}_{ij,\mathrm{n}}^{\mathrm{c}}\right|$, with ${\mu }_{\mathrm{r}}$ standing for the coefficient of rolling resistance. Furthermore, ${\mathit{\omega }}_{ij}={\mathit{\omega }}_i-{\mathit{\omega }}_j$, $\mathrm{\Delta }t$ denotes the time step, and $k_{\mathrm{r}}=k_{\mathrm{t}}\sqrt{{\delta }_{ij,\mathrm{n}}}{R}_{\mathrm{e}\mathrm{ff}}^2$ (17, 51, 52). The magnitude of the spring torque ${\mathbf{M}}_{ij}^{\mathrm{r}}$ is, thus, upper bounded by the full mobilization torque ${\mu }_{\mathrm{r}}{R}_{\mathrm{e}\mathrm{ff}}\left|{\mathbf{F}}_{ij,\mathrm{n}}^{\mathrm{c}}\right|$, which means that rolling resistance behaves perfectly plastic beyond this upper bound.

The integration was performed using the open source Discrete Element Method particle simulation software LIGGGHTS (which stands for LAMMPS Improved for General Granular and Granular Heat Transfer Simulations) (52), along with the extensions introduced in ref. 15, to account for the computation of the dissipative term in the contact model (Eqs. 10 and 11), following ref. 44, and the attractive particle interaction force model (15, 47, 48) defined in Eq. 13. The integration time step $\mathrm{\Delta }t$ has been about $10\%$ of the collision time $T_{\mathrm{c}\mathrm{o}\mathrm{l}}$ between two particles, defined from Hertz’s elastic theory for the undamped, noncohesive situation; i.e., $T_{\mathrm{c}\mathrm{o}\mathrm{l}}\approx 3.21{\left(m_{\mathrm{e}\mathrm{ff}}/k_{\mathrm{n}}\right)}^{2/5}\cdot v_{\mathrm{i}\mathrm{m}\mathrm{p}}^{-1/5}$ (12), where $m_{\mathrm{e}\mathrm{ff}}=m/2$, and $v_{\mathrm{i}\mathrm{m}\mathrm{p}}$ denotes the relative normal velocity of the particles at impact.

The values of Young’s modulus, Poisson’s ratio, particle density, and surface energy density, listed in Table 1, are consistent with glass beads (15, 49), and we take the coefficients of sliding and rolling friction following calibration results from previous DEM simulations (6, 14–17, 46, 53, 54). The collisions of the particles with the frictional bottom wall are modeled using the same material parameters, while the wall is regarded as cohesionless and one of the collision partners has infinite mass and radius (15).

Furthermore, the normal damping coefficient $A_{\mathrm{n}}$ is computed based on the coefficient of restitution ${\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{f}}$ associated with a reference collisional velocity $v_{\mathrm{i}\mathrm{m}\mathrm{p}}$, assuming only contact forces (no cohesion) (44). Based on previous calibration of the coefficient of restitution for use in DEM studies of the angle of repose (53, 54), here we adopt ${\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{f}}$ in the range between 0.2 and 0.4, for $v_{\mathrm{i}\mathrm{m}\mathrm{p}}=1.0$ m$/$s. All results presented in this article refer to ${\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{f}}=0.2$, but we have found that using ${\epsilon }_{\mathrm{r}\mathrm{e}\mathrm{f}}=0.4$ still leads to a reasonable quantitative agreement between our numerical predictions for ${\theta }_{\mathrm{r}}$ and the experimental data in Fig. 2. However, we reinforce that the main contribution of our article consists of a mathematical expression that is robust with respect to details of interparticle interactions, encoded in the parameters ${\mu }_{\mathrm{e}\mathrm{ff},\infty }$ and $\beta$. The role of various ingredients that have been neglected in our simulations, such as torsion (13) and long-range electrostatic interactions, should be elucidated in future work.

Fig. 5 displays the results from numerical experiments to estimate the packing fraction $\varphi$ as a function of $d$ and $g$. Each of these experiments consists of releasing from rest a monodisperse packing of $20\%$ initial volume fraction composed of $\mathrm{11,039}$ (initially noncontacting) particles randomly distributed over the simulation domain—a vertical rectangular box of horizontal dimensions $L_x\times L_y$, with $L_x=L_y=17\hspace{0.17em}d$ and height $100\hspace{0.17em}d$. Boundary conditions are the same as in the experiments of the angle of repose, except that periodic boundary conditions are applied in the horizontal directions (15). Once the particles have settled due to the action of gravity, we compute the packing density of the granular column using $\varphi =n_{\varphi }\cdot \left(\pi /6\right)\hspace{0.17em}{d}^3/\left[L_x\cdot L_y\cdot \left(H_u-H_l\right)\right]$, where $n_{\varphi }$ is the number of particles with center-of-mass position between $H_l=0.3z_{max}$ and $H_u=0.7\hspace{0.17em}{z}_{max}$, with $z_{max}$ denoting the center-of-mass position of the highest particle in the granular column (15). We find that $\varphi$ differs negligibly over different realizations of the same experiment (each realization with different initial random positions of the particles). Therefore, each value of $\varphi$ presented in Fig. 5 has been obtained based on one realization for the corresponding value of $d$ and $g$, but by averaging over the four vertical subcolumns obtained by halving $L_x$ and $L_y$ (15), the associated SD corresponding to the error bars in Fig. 5 (i.e., the same protocol as in ref. 15 but including the effect of gravity in the investigation).