# Sink versus tilt penetration into shaken dry granular matter: The role of the foundation

L. Alonso-Llanes, G. Sánchez-Colina, A J Batista-Leyva, C. Clément, E. Altshuler, R. Toussaint

<figure><img src="image_1.png" /></figure>

## To cite this version:

L. Alonso-Llanes, G. Sánchez-Colina, A J Batista-Leyva, C. Clément, E. Altshuler, et al., Sink versus tilt penetration into shaken dry granular matter: The role of the foundation. Physical Review E, 2022, 105 (2), pp.024903. 10.1103/PhysRevE.105.024903. hal-03708677

HAL Id: hal-03708677

https://hal.science/hal-03708677

Submitted on 29 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Sink vs. tilt penetration into shaken dry granular matter:
the role of foundation

L. Alonso-Llanes,¹,²,* G. Sánchez-Colina,¹ A. J. Batista-Leyva,¹,³

C. Clément,² E. Altshuler,¹ and R. Toussaint²,⁴,†

¹Group of Complex Systems and Statistical Physics,
Physics Faculty, University of Havana, 10400 Havana, Cuba

²Université de Strasbourg, CNRS, Institut Terre et Environnement de Strasbourg,
UMR7063, 67000 Strasbourg, France

³Instituto Superior de Tecnologías y Ciencias Aplicadas
(InSTEC). University of Havana, 10400 Havana, Cuba

⁴SFF PoreLab, The Njord Centre, Department of Physics,
University of Oslo, P.O. Box 1074 Blindern, 0316 Oslo, Norway

(Dated: June 29, 2022)
Abstract

We study the behavior of cylindrical objects as they sink into a dry granular bed fluidized due to lateral oscillations. Somewhat unexpectedly, we have found that, within a large range of lateral shaking powers, cylinders with flat bottoms sink vertically, while those with a “foundation” consisting in a shallow ring attached to their bottom, tilt besides sinking. The latter scenario seems to dominate independently from the nature of the foundation when strong enough lateral vibrations are applied. We are able to explain the observed behavior by quasi-2D numerical simulations, which also demonstrate the influence of the intruder’s aspect ratio. The vertical sink dynamics is explained with the help of a Newtonian equation of motion for the intruder. Our findings may shed light on the behavior of buildings and other man-made constructions during earthquakes.

# I. INTRODUCTION

The Kocalei earthquake occurring on August 17, 1999 affected in various ways many constructions in the city of Adapazari, Turkey. Following observers, some buildings sank vertically into the soil, others tilted, and some even suffered lateral translation over the ground [1–3]. This case illustrates well the diversity of damage that earthquake fluidization of soils may cause to man-made structures [4].

Liquefaction in the ground may be triggered dynamically by waves emitted during earthquakes, generating cyclic shear stresses that lead to the gradual build-up of pore water pressure [5, 6]. The shaking produced by seismic events is a trigger for extensive liquefaction, as was observed recently in Belgium [7]. Ground fluidization [8, 9] has been investigated in different kinds of media like sand [9], dry granular soils [10] and sediments [11]. Of immediate interest for engineering and for the geosciences is to understand how man-made structures such as buildings, and massive rocks laying on granular soils respond to fluidization associated to seismic waves.

Granular matter itself displays a variety of puzzling phenomena [12–25], but during the last decade or so, our understanding of the dynamics of objects penetrating into granular

<footnote>* lalonso@fisica.uh.cu</footnote>

<footnote>† renaud.toussaint@unistra.fr</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. Experimental setup. At the upper right, we have illustrated the intruder consisting in a cylinder with ring.</figcaption></figure>

media has advanced quickly [26–43]. While laterally shaken granular beds have received a certain degree of attention [44, 45], the performance of objects initially laying on the surface of a granular bed submitted to lateral shaking has been rarely studied [46–50].

In this paper we perform systematic experiments associated to the latter scenario, which may help understanding the performance of human constructions and rocks laying on granular beds during earthquakes. In particular, using a cylinder as a simplified model, we study its settling dynamics on a granular bed submitted to lateral vibrations. Somewhat unexpectedly, we have found that, within a large range of lateral shaking powers, cylinders with flat bottoms sink vertically, while those with a “foundation” consisting in a shallow ring attached to their bottom, tilt besides sinking. The latter scenario seems to dominate independently from the nature of the foundation when strong enough lateral vibrations are applied. Quasi-2D simulations mimicking the experiments were also performed. The settling dynamics of the simulated intruders, with or without foundation, reproduces the corresponding experimental results. Our simulations also reveal how the difference in force-chain distributions between flat and non-flat bottom cylinders produces different torques justifying the two types of penetration. In addition, we present a simple phenomenological model that reproduces well the sinking dynamics and helps understanding how the tilting process influences the sinking one.
## II. EXPERIMENTAL

The penetration experiments were performed on a granular bed contained in a test cell of approximately 25 × 25 × 25 cm³ filled with Ugelstad spheres of non expanded polystyrene with a bulk density 1050 kg/m³, and diameter 140 µm (monodisperse within a 1 percent), type Dynoseeds ®, produced by Microbeads, Norway. The box was horizontally shaken at different amplitudes of motion (A), and a frequency (f) of 5.0 Hz (a value commonly found in seismic waves), using a TIRA ®TV51120-M shaker, see figure 1. By controlling the voltage of the shaker input signal we varied the amplitude of the oscillations up to a maximum value corresponding to a peak ground acceleration of A(2πf)² ≈ 12.2 m/s² [51]. This acceleration range covers most potentially damaging earthquakes, from weak to strong [52], though there has been reports of larger peak ground accelerations [53]. The time the shaker needs to reach the steady state depends on the dimensionless acceleration, being longer for larger accelerations. The time intervals can range from one (0.2 s) to three periods.

Two types of intruders were used in the experiments: (a) a hollow 3D printed cylinder of 44 mm diameter, 44 mm height $h_c$ (external dimensions), and 5 mm thick walls, and (b) the same cylinder with a ring of 5 mm height and 3 mm thickness glued to its bottom (illustrated in the upper right corner of Fig. 1). Intruders (a) and (b) will be called “No-ring” and “Ring”, respectively, from now on. Their masses were adjusted with ballast in such a way that their densities matched the average effective density of the granular medium, which was measured as 430 kg/m³. As far as the ballast used has a density near the effective density of the granular material, it was almost evenly distributed inside the cylinder. Note that, using a flat bottom cylinder and a ring-like bottom cylinder, we are modifying the “foundation” of our intruder.

A digital camera *Hero 2* made by GoPro ® was fixed to the electromagnetic shaker, in such a way that it could take a video of the sinking process from an oscillating reference frame locked to the test cell, as proposed in [48]. This method allowed a much better quality of the cylinder’s images, and made easier their digital processing. Videos were taken at a maximum rate of 120 frames per second, with a resolution of 1920 × 1080 pixels.

The images were processed as follows. We first converted the videos to image sequences in *.jpg format, and cropped each picture, excluding irrelevant space. Then, the images were binarized through an appropriate threshold. Using the tool *regionprops* from *MatlabR2014a*
<figure><img src="image_4.png" /><figcaption>FIG. 2. Sketch of sinking and tilting processes. The top row illustrates the sinking process of a No-ring cylinder in three moments during the experimental interval from $t = 0$ to a final time $t = \Delta t$. The bottom row shows the same temporal sequence for a Ring cylinder, which tilts in addition to sink.</figcaption></figure>

(c), we identified and assigned coordinates to several bright marks we had glued to certain points of the cylindrical intruder. The coordinates of the marks were used to calculate the position of the intruder's geometrical center and inclination relative to the vertical in each picture. In some experiments where the sinking was particularly big, it was difficult to obtain the tilt angle, since part of the marks sank below the level of the sand surface, and they were impossible to follow. In such cases the upper border of the cylinder was identified using the Matlab's tools *find* and *bwtraceboundary*, and then fitted to a polynomial using the function *polyfit*. The fit was used to find the inclination. In the case of experiments ending in a very inclined position, the reference to calculate the inclination was the cylinder's corner above the sand surface, that was identified as the intersection of the two polynomial fits of the upper and one lateral borders of the cylinder.

As the cylinder oscillates due to the vibration of the box, it is difficult to determine the final position, particularly when there is a big tilting. Then, in order to determine the sinking depth and tilting, we observe in the videos the onset of a cyclic movement of a reference point in the cylinder. Then, the final position could be measured in the frames filmed after the shaker was stopped.

The experimental protocol can be described as follows: (I) preparing the granular medium by stirring it evenly with a long rod, (II) gently depositing the cylinder in the upright position on the free surface of the granular bed, (III) turning ON the camera, (IV) switching ON the shaker after setting the desired frequency and amplitude (V) turning OFF the shaker and
<figure><img src="image_5.png" /><figcaption>FIG. 3. (color online) Snapshots of the initial (a) and final (b) positions of a No-ring intruder in a typical quasi-2D simulation using a shaking frequency of 5 Hz.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 4. (color online) Snapshots of the initial (a) and final (b) positions of a Ring intruder in a typical quasi-2D simulation using a shaking frequency of 5 Hz.</figcaption></figure>

the camera after the penetration process had visibly finished.

In Fig. 2 we define the main parameters describing the sinking process of a No-ring
cylinder (upper row), and the tilting and sinking of a Ring cylinder (bottom row), during
the experimental lapse, defined as $\Delta t$. As the figure indicates, in the following we will call
$h$ the penetration of the geometrical center at a time $t$ and $\Delta h$ the final penetration at
time $t = \Delta t$. Note that both magnitudes are defined as the vertical displacement of the
geometrical center of the cylinder (without taking the ring into account). In the same way
we will call $\theta$ the inclination of the intruder at time $t$ and $\Delta\theta$ the final inclination at $t = \Delta t$.

We also explored the phenomenology through numerical simulations. They were based
on a discrete element method code (DEM) for the computation of granular systems [20,
22, 24, 49, 50, 54, 55]. We modeled a quasi-2D granular medium, made of finite-sized hard spheres with radii between 1.0 and 1.5 mm, to avoid the crystallization effect. The medium contains 4000 particles and is prepared by placing the latter randomly in a space 30 cm wide and 25 cm high and then allowing them to settle under the action of gravity $g = 9.81 \text{ m/s}^2$. Once the medium reaches equilibrium, it occupies a virtual space 30 cm wide and about 8 cm high, laterally delimited by flat walls that define the Hele-Shaw cell. The simulation box used was created narrow in order to have a single plane of particles in the direction perpendicular to the images shown in Figures 3 and 4. The components of the velocities and forces along this direction are set to zero at each time step. To mimic the experimental conditions, we simulate particles of density 1050 kg/m$^3$.

The two intruders are made of cohesive particles. One is a square of 40 mm side, made of $N = 1681$ particles with diameter 1 mm, placed in a quasi-2D square arrangement, which simulates the No-ring intruder of the experiments. The second one is also a square of 40 mm side to which two "small feet" are attached. Each foot has a size of $4.5 \times 2.7 \text{ mm}^2$, so the simulated Ring intruder contains a total of $N = 1705$ particles. The density of the spheres $\rho_p$ which form the intruders is adjusted so that the bulk density of the rigid body matches the effective density of the granular medium. The latter is calculated once the medium has settled down and is stable, and was always found to be around $\rho_m = 566 \text{ kg/m}^3$. Then, the density of the particles forming the intruder is obtained as $\rho_p = \rho_m V_i / N V_p$, where $V_i$ is the volume of the intruder and $V_p$ the volume of a sphere.

Once our granular medium is created, we place the intruder 1 mm above the medium, with its bottom parallel to the horizontal direction. We release it, under the action of the force of gravity, and wait until the whole system becomes motionless (i.e. its total kinetic energy reaches a value under $10^{-7} \text{ J}$). Then, we apply horizontal oscillations of different amplitudes and a frequency of 5 Hz ($\Delta t \approx 8 \text{ s}$) to the walls of the cell and compute the time evolution of the position and tilting angle of the intruder. The amplitudes were chosen in such a way that the dimensionless acceleration $\Gamma = A(2\pi f)^2/g$ had the values of 0.16, 0.25, 0.5, 1.0, 1.25 and 1.5.

The contact between spheres was modeled as a linear spring-dashpot $F_{ij} = (k_n \delta \mathbf{n}_{ij} - m_{eff} \gamma_n \mathbf{v}_n) - (k_t \Delta \mathbf{S}_t + m_{eff} \gamma_t \mathbf{v}_t)$ [54, 56], where $k$ and $\gamma$ are the elastic and viscoelastic damping constants, $\delta \mathbf{n}_{ij}$ is the overlap distance along the line connecting the centers of the two spheres, and $\mathbf{v}$ their relative velocity. $\Delta \mathbf{S}_t$ is the tangential displacement vector
between two spheres, which is truncated to satisfy a frictional yield criterion, and $m_{eff} = m_i m_j / (m_i + m_j)$ is the effective mass of two spheres of mass $m_i$ and $m_j$. We considered normal ($n$) and tangential ($t$) forces components between the particles and, in order to model hard spheres that interact on contact (i.e. spheres whose deformation during collisions is less than a small fraction of their radii), we used the following parameter values: $k_n = 1.2 \times 10^7$ N/m, $k_t = 2/7k_n$, $\gamma_n = 12$ s$^{-1}$ and $\gamma_t = 0.1\gamma_n$. The interaction force between the walls and the particles touching them is the same as the corresponding for two particles but considering the wall of infinite radius and mass (flat wall). The microscopic friction coefficient between spheres, and between spheres and boundaries, was taken as $\mu = 0.3$. The time step $dt$ was chosen to guarantee that there are at least 50 steps during one characteristic time of a collision $dt = t_c/50$, where $t_c = \pi / \sqrt{(k_n/m_{eff}) - \gamma_n^2}$.

Figures 3 and 4 show the initial and final positions of both types of intruders in two typical runs. Fig. 3 indicates that the No-ring cylinders are slightly inclined, while in Fig. 4 the large inclination of the Ring one becomes obvious.

We also performed an additional set of simulations aimed at elucidating the influence of the intruder's aspect ratio in the sink - tilt behaviour. Two intruders, one No-ring with aspect ratio 1.125, and a Ring with aspect ratio 1 (see Fig. 5), were submitted to the same range of dimensionless accelerations. Note that the new intruders (in dark gray in Fig. 5) have the same dimensions of the former ones (in light gray) but the geometry of the bottom is interchanged.

The new intruders are also rigid bodies made of 1 mm diameter cohesive particles and their bulk densities also correspond to that of the granular medium. The new Ring intruder is composed by $N = 1498$ particles forming a 40 mm wide and 35.5 mm high rectangle to which two $4.5 \times 2.7 \text{ mm}^2$ feet were added. Note that the size of these feet and those of the Ring intruder are the same. The new No-ring is a rectangle of 40 mm wide and 44.5 mm high composed of $N = 1845$ particles.
<figure><img src="image_7.png" /><figcaption>FIG. 5. Comparison of the geometry of the intruders used to test the influence of the aspect ratio. In light gray the intruders previously described and in dark gray the new ones.</figcaption></figure>

### III. RESULTS AND DISCUSSION

#### A. Sink vs. tilt penetration in experiments

Figure 6(a) shows the time variation of the sinking depth for selected values of the dimensionless acceleration $\Gamma = A(2\pi f)^2/g$ (where $g = 9.81 \text{ m/s}^2$ is the gravitational acceleration and $A(2\pi f)^2$ is the horizontal peak acceleration of the sand box) for No-ring cylinders. It is easy to see that the penetration of the No-ring cylinders follows a common pattern for all the accelerations. A first process of fast sinking is followed by a slow creep. Only the penetration depth increases with $\Gamma$. In this figure we do not show the total creep process, due to its long duration. As the height of the cylinder is $h_c = 44$ mm, it is possible to check from Fig. 6(a) that, for a dimensionless accelerations of 1.24, the cylinder sinks completely. An important characteristic of the sinking process in this type of cylinder is that the intruder penetrates the granular medium with almost no tilting, and a final inclination smaller than 5°.

Figure 6(b) is similar to the previous one, but measurements were performed with Ring cylinders. The general features of both graphics are similar, but there is a difference, that will be better observed in the following figures: the dimensionless acceleration at which the cylinder sinks completely in the medium is bigger for the Ring cylinders than for the No-ring ones.

Figure 6(c) presents the time evolution of the tilting angle for a Ring cylinder, a process that occurs simultaneously with the sinking. The sinking and tilting dynamics of Ring cylinders is more irregular than that of the No-ring ones. This is illustrated in Fig. 6(b) and (c), even after being submitted to a moving average process, to get a smoother graph.

No-ring cylinders tend to sink vertically as the granular soil is fluidized by horizontal
<figure><img src="image_8.png" /><figcaption>FIG. 6. (color online) Experiments. Time evolution of penetration depths and tilt angles. Time dependence of the penetration depth of a No-ring cylinder (a), the penetration depth of a Ring cylinder (b) and the tilting angle of a Ring cylinder (c), for different dimensionless accelerations. The long-time creep process is not completely shown. The tilting angle of No-ring cylinders is not displayed, due to the fact that it oscillates around angles not larger than 5° relative to the vertical direction.</figcaption></figure>

shaking, while cylinders with rings tend to tilt. Figure 7 quantifies the differences between the initial and final stages of the process, for almost all the range of accelerations our experimental setup was able to reach.

Figure 7(a) shows sink data for No-ring cylinders. As can be seen, for dimensionless accelerations up to Γ = 0.27, there was no significant penetration of the intruder into the
<figure><img src="image_9.png" /><figcaption>FIG. 7. Experiments. Sinking and tilting: heights and angles for different dimensionless accelerations $\Gamma$. Final sink heights for No-ring cylinders (a) and for Ring cylinders (b). Final tilt angles for Ring cylinders (c). Tilt angles of No-ring cylinders are not shown for the same reasons of the previous figure.</figcaption></figure>

granular bed. Vertical penetrations started to increase significantly above this acceleration, reaching a plateau around $\Gamma \approx 0.7$. At the plateau, the cylinder has sunk completely, but stays “floating” into the fluidized granular medium, as expected for an object isodense relative to it, so there is no further sinking.

In Fig. 7(b) the sinking process of the Ring cylinders is summarized. Though the low acceleration part is similar to Fig. 7(a), now the plateau is not observed for the range of accelerations recorded. Notice that from the depth reached at $\Gamma \approx 1.2$, approximately the
height of the cylinder $h_c$, it would not sink any further, and that from this value of $\Gamma$ onwards a plateau would appear.

Figure 7(c) shows the tilt data for Ring cylinders. No significant tilting is observed for $\Gamma$ smaller than approximately 0.25. With the increase of the dimensionless accelerations, the cylinder significantly tilts, increasing abruptly the tilting angle with $\Gamma$, until it slows down at $\Gamma \approx 0.7$. We do not show the tilting angle of No-ring cylinders, because it is always smaller than 5°, with a random distribution of values around the vertical direction.

Figures 7(b) and (c) are closely related, because they are two descriptions of the same process: the motion of Ring cylinders into the granular medium, that includes both sinking and tilting. The fact that at the accelerations shown in this figure the plateau in the sinking depth is not reached while for the tilting angle at higher values of $\Gamma$ the inclination almost saturates, could be explained by the increase of the friction of the intruder with the granular medium when the tilting angle increases. At $\Gamma \approx 0.7$ the intruder has reached a large inclination, but is not completely submerged into the medium. An increase in the acceleration does not increase significantly the angle, because the resulting torque has diminished due to the influence of both sinking and tilting, but the increase in fluidization helps further sinking, until most of the cylinder is submerged.

The sinking process can be understood taking the experimental results in Ref. [46] into account. When the system is submitted to lateral shaking, a fluidized layer appears in the upper part of the granular cell. This layer reaches a depth $h_f$ that depends on the dimensionless acceleration $\Gamma$. Below this layer exists a “solid” layer. For accelerations in the range spanned in our experiments, $h_f$ varies almost linearly with $\Gamma$ (see Fig. 3(a) in Ref. [46]), so we can write

$$h_f(\Gamma) = \alpha(\Gamma - \Gamma^*); \Gamma > \Gamma^* \qquad (1)$$

where $\Gamma^*$ is the onset of fluidization and $\alpha$ is the slope of the linear dependence. If $\Gamma \le \Gamma^*$ the depth of the fluidized layer is zero.

Then, at low values of $\Gamma$ the granular medium is not fluidized, and the cylinder almost does not sink (merely 5 mm at $\Gamma = 0.27$; see Fig. 7(a)). For accelerations above the fluidization threshold, the cylinder sinks until it gets in contact with the solid layer. The larger is the acceleration, the deeper is that layer, so the bigger is $\Delta h$. But as soon as the solid layer appears at a depth larger than the cylinder’s height, it does not sink further:
instead, it “floats” due to isodensity with the sand, so a plateau is reached.

According to reference [57], $\Gamma^*$ can be taken as proportional to the friction coefficient $\mu$
between the cylinder and the granular medium. In these experiments we can approximate
$\mu \approx 0.3$, which is the value we use in the simulations. The authors also conclude that the
final depth of intrusion depends on isostasy, and on the severity of shaking. It can be entirely
determined by isostasy, when the shaking completely unjam the medium and suppresses the
average friction around the intruder [50].

To better understand the differences in the dynamics of both types of intruders, we
performed numerical simulations and their results are described below.

## B. Sink vs tilt in quasi-2D numerical simulations

Figure 8 shows the time dependence of the penetration depth (a) and tilting angle (b) for both types of intruders at the dimensionless acceleration $\Gamma = 1.0$. In both figures the thick curves represent the average value of six repetitions varying the initial conditions and the surrounding zone represents $\pm 1$ standard deviation.

Regarding the vertical sinking in Fig. 8(a), we do not observe major changes between Ring
and No-ring intruders; both types of intruders sink less than in the experiments. This may
be related with the lower dimensionality of the simulations relative to the real experiment.
Quasi-2D granular media allow less choices of readjustment than in 3D: they are easily
jammed, which makes it more difficult for an object to sink. Moreover, the size ratio of
the intruder over the particles is 8 times smaller in the simulations than in the experiments
(experiments: 44 mm/0.140 mm $\approx$ 300; simulations: 40 mm/1 mm= 40), which means
that if one particle is stuck under the intruder during the simulations, it will slow down the
intruder significantly more than if the particle were 8 times smaller.

Figure 8(b) indicates that the presence of a foundation at the bottom of the intruder
causes a large tilting. Indeed, for the shaking with no ring, the intruder tilting angle is
around 10°, but during the shaking with ring, the intruder tilts fast, reaching an angle
around 50°. This resembles what happens in the experiments (see Fig. 4): the intruder
almost ends up lying on one of its sides. Of course, the tilting is also limited by the
diminished dimensionality in the quasi 2D simulations.

Figure 9 compares the penetration depth reached for both types of cylinders at different
<figure><img src="image_10.png" /><figcaption>FIG. 8. (color online) Simulations. Time dependence of sinking depth (a) and tilting angle (b) for both types of intruders with $\Gamma = 1.0$. The central lines are the average of six simulations, while the colored bands cover $\pm 1\sigma$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 9. (color online) Simulations. Final depth reached for No-ring (squares) and Ring (circles) cylinders as a function of the dimensionless acceleration. Symbols represent the average of six experiments and the error bars $\pm\sigma$.</figcaption></figure>
<figure><img src="image_12.png" /></figure>

FIG. 10. (color online) Simulations. Maximum tilting angle reached for No-ring and Ring cylinders as a function of the dimensionless acceleration. Symbols represent the average of six simulations and the error bars $\pm\sigma$.

values of $\Gamma$. The conclusions obtained from Fig. 8 are valid for all the dimensionless accelerations used in the simulations: there are no significant differences in the final sinking depth between both types of intruders.

In Fig. 10, on the contrary, the difference in tilting angles between the two types of intruders is clearly seen. For all the values of $\Gamma$ the simulated Ring intruder tilts more than the No-ring one. For $\Gamma=1.50$, the No-ring tilts up to an angle that is closer to the Ring's one, corresponding to preliminary observations found in experiments with frequencies above 5 Hz for dimensionless accelerations $\Gamma >1.25$.

According to our simulations, the difference in tilting between intruders lies in that one type of intruder, the No-ring one, is somehow more capable of rectifying its rotation during sinking, while the other, the Ring intruder, is not. This rectification can be understood as the process of returning to, or recovering, the initial rotation angle once one oscillation of the cell has concluded and, as can be seen in the temporal evolution of $\theta$ (Fig. 11(a)), the difference in the rotation angles between the No-ring and Ring intruders is produced by a non-rectifying cumulative process taken by the latter.

To understand why the tilting dynamics is affected by the presence of the legs, which makes the No-ring intruder able to further rectify its rotation – at least for the values of gamma between 0.25 and 1.25–, we calculate from the simulations the torque about the center of mass and the angular velocity.

Fig. 11 shows the time evolution of $\theta$, the angular velocity, the torque about the center
<figure><img src="image_13.png" /></figure>

FIG. 11. (color online) Simulations. Temporal evolution of the calculated magnitudes for the intruders during the first second of a simulation using $\Gamma = 1.0$. From top to bottom: the evolution of the rotation $\theta$ (a), the angular velocity $\Omega$ (b), the torque about the center of mass $\tau$ (c) and finally, the position of the cell (d) showing the oscillation described by it. The vertical line indicates the instant $t = 1.84$ T. $\tau_{max}$ is maximum torque on a non-tilting but horizontally accelerated No-ring intruder.

of mass and the oscillation of the cell for the two types of intruders in a simulation with $\Gamma = 1.0$. Focusing on the $\theta$ curve, the difference in tilting can be noticed during the first second of the simulation (as in Fig. 8 (b) for the averaged values) as well as the aforementioned Ring intruder non-rectification process. Unexpectedly, the values of the torques shown in Fig. 11 (c) are very similar for the two intruders, contradicting the intuitive idea that the sole presence of the legs would produce higher torques about the Ring intruder's center of mass making it rotate more. However, the slight differences may affect, to a greater or lesser extent, the rotational movement of the latter with respect to the No-ring intruder (see Fig. 11 (b)).

In general, during the first half of a cell oscillation, both intruders rotate in the same direction –counterclockwise due to the fact that the granular bed moves from left to right. But, during the second half, some forces appear in the bottom and/or the inside of the
<figure><img src="image_14.png" /></figure>

FIG. 12. (color online) Simulations. Pressure field in the granular medium represented as a jet color map for No-ring (a) and Ring (c) intruders. The color of the intruders represents their angular velocity using a seismic color map. In (b) and (d) the outer segments of the intruders are represented as a whole with a color corresponding to the resulting torque about the center of mass that is generated on all grains belonging to them. Notice that for the Ring intruder each leg is divided into 3 segments. The grains in these segments are magnified for better viewing. The resulting torque about the center of mass is represented as the interior color of the intruder. For both angular velocity and torque, blue color represents clockwise.

Ring intruder's right leg that do not allow it to rotate in the same way as the No-ring does. In some cases, these forces completely prevent it from rotating clockwise, as in both first oscillations shown in Fig. 11. This process is prone to occur in each of the oscillations during the simulation and its repetition causes the differences in inclination observed for the two intruders after 8 simulated seconds (see Fig. 10).

Fig. 12 illustrates in more detail what is described above. In (a) and (c) it shows the pressure field in the granular medium at $t = 1.84$ T (time indicated in Fig. 11 by the vertical line), where the force chains are represented using a sequential color map. In them, the color of the intruders represents their angular velocity, which in turn is displayed as a diverging color map where blue indicates clockwise rotations. Figures 12 (b) and (d) show the contribution to the torque about the center of mass of each of the intruder segments as a result of the forces acting on them. In this case, the color of the intruders is associated with the resulting torque about the center of mass. Note that all the grains of each outer edge in
<figure><img src="image_15.png" /><figcaption>FIG. 13. (color online) Simulations. Grain tracking. Initial (left) and final (right) positions of the No-ring (top) and Ring (bottom) intruders as well as the grains of the granular medium next to them. The grains right below the intruders are represented in blue while the yellow and cyan colors represent the grains initially at the left and right side of the intruders, respectively. Note that for the Ring intruder, the blue grains initially located between the legs remain in the same region until the end of the simulation</figcaption></figure>

(b) and (d) are represented with one color corresponding to the resulting torque about the center of mass obtained from the torques of all the particles on this same edge, though only the parts of the edge in contact with the granular medium are effectively interacting. These figures help identify what type of torque, clockwise or counterclockwise, is generated in each part of the intruders (including the legs), how it is generated, and how representative it is for the resulting torque about the center of mass. In this particular case, it is observed for the Ring intruder that the force associated with the lower part of the right leg is responsible for almost canceling the torque about the center of mass. Therefore, the Ring intruder remains inclined while the No-ring one rotates back to its original position.

Those forces acting on the Ring intruder's legs could be associated with the grain jamming between them. Fig 13 shows that almost all the blue grains that were initially in the region between the legs remain there throughout the simulation. Furthermore, it can be seen that some of the grains (cyan) around the Ring intruder accumulate between its legs: this is because during the first part of the oscillations these grains tend to move towards the legs, however, they cannot leave them during the second part. In contrast, grains in the region
below No-ring intruders have more freedom to exit. The previous process suggests that the Ring intruder along with the grains between the legs could be treated as a No-ring intruder with increased friction at the bottom. This increase in friction will make the grains underneath more likely to get stuck, preventing them from coming out and thus creating force chains capable of stopping the clockwise rotation (restoring towards vertical position) of the intruder.

A final observation from the simulations for dimensionless accelerations of $\Gamma = 1.5$ is that the No-ring intruder rotates almost 90 degrees in a 100% of the cases where it reached 45 degrees, doing so in an abrupt way. The Ring intruders, however, show a constant growth up and, in some cases (after reaching 90 degrees), increase the slope reaching values of up to 180 degrees. It is worth noting that once Ring intruders turn 90 degrees they begin to resemble No-ring ones as the presence of the legs loses importance in the penetration dynamics. Behaviors such as those described before were not observed experimentally since the $\Gamma$ values used in the experiments did not exceed 1.24 due to technical limitations of the shaker used.

The experimental findings are explained not only by the numerical simulations, but also by a Newtonian model developed in the Appendix 1. This model is based in the force balance on a cylinder sinking in a granular medium. The forces considered in the 1D model are gravity, a frictional force proportional to velocity and a pressure like force, proportional to the depth $h$, as expressed in Eq. (A.13). Though this model does not include the degree of freedom associated with tilting, the consideration of its influence in the lineal and surface dimensions of the intruder is enough to explain why a tilted intruder reaches a final depth smaller than that reached by a non tilted intruder for a given $\Gamma$.

Now, we examine the influence of the intruder's aspect ratio on the penetration dynamics. The results of the simulations performed with the intruders of different aspect ratios are summarized in Fig. 14. Fig. 14(a) compares the dependence of the final tilt of two intruders with the same aspect ratio (1), one with legs and the other with a flat bottom. With the exception of the smaller values of $\Gamma$ where no noticeable differences are apparent, the intruder with legs always tilts more. This behavior is repeated in the results of Fig. 14(b) for a larger aspect ratio (1.1125): again the intruder with legs has a final tilt larger than the one with a flat bottom. Therefore, the above suggests that the presence of a ring at the bottom of the intruder causes a higher final tilt, although it might become less important with
<figure><img src="image_16.png" /><figcaption>FIG. 14. Simulations. Final tilt angles at different $\Gamma$ for two intruders with (a) aspect ratio 1 and different foundations; (b) aspect ratio 1.1125 and different foundations; (c) a flat bottom and different aspect ratios and (d) a ring at the bottom and different aspect ratios.</figcaption></figure>

increasing aspect ratio. Figures 14(c) and (d) compare intruders with the same foundation and different aspect ratios. Both figures lead to the same conclusion: the larger the aspect ratio, the larger the tilt angle. Summarizing the results of Fig. 14, the intruder with higher aspect ratio and ring placed on the bottom has the largest tilt angle for all dimensionless accelerations, and the intruder with flat bottom and lower aspect ratio has the smallest tilt angle. Interestingly, the other two intruders show approximately equal final tilt angle values for equal values of $\Gamma$.
## IV. CONCLUSIONS

We have studied the behavior of cylindrical objects as they sink into a dry granular bed fluidized by horizontal oscillations, as a model system to understand the effects of earthquake-related fluidization of soils on human constructions and other objects like rocks.

We have found that, within a relatively large range of lateral shaking amplitudes at a frequency of 5 Hz, cylinders with flat bottoms sink vertically, while those with a “foundation” consisting in a shallow Ring attached to their bottom, tilt laterally besides their vertical sinking.

We have been able to mimic the above described behaviors by quasi-2D numerical simulations. With their help we found that these differences are not necessarily due to the sole existence of the legs that generate bigger torques about the center of mass. Instead, they can be associated with the jamming of the particles in the region between the legs, which can increase the friction at the bottom of the Ring intruder generating force chains capable of preventing the total recovery of its initial angle of rotation. Numerical experiments also helped to clarify the influence of the intruder aspect ratio on the tilt dynamics: of two intruders with the same foundation, the one with higher aspect ratio will have a larger tilt angle. So, the aspect ratio and the foundation type cooperate to establish the penetration dynamics of the intruder.

We have also reproduced the vertical sink dynamics of cylinders with a flat base using a Newtonian equation of motion for an object penetrating a fluidized layer of granular matter, where the granular effective density increases with depth, eventually reaching a solid phase. The same model allows to understand the sinking even in the present of tilting (Appendix 1).

Finally, it is worth noting that preliminary experimental data and quasi-2D numerical simulations suggest that, when strong enough lateral shaking is applied, the tilting scenario tends to dominate regardless the nature of the intruder's foundation.

## ACKNOWLEDGMENTS

L. A. thanks the Photovoltaic Research Lab of the University of Havana for allowing the use of its facilities. E. A. drew inspiration from the late M. Álvarez-Ponte. We acknowledge support from Project 29942WL (Fonds de Solidarité Prioritaire France-Cuba), from the EU
ITN FlowTrans, and from the INSU. We thank the SCAC of the French Embassy in Havana and the University of Strasbourg for their support. RT acknowledges the support of the Research Council of Norway through its Centres of Excellence funding scheme, Project No. 262644. We thank Mustapha Meghraoui, Einat Aharonov, Knut Jørgen Måløy, Eirik G. Flekkøy, for fruitful discussions. This research was made in the frame of the University of Havana's institutional project "Medios grandes: creando herramientas para prevenir catástrofes".

# Appendix

## 1. Sink dynamics: a phenomenological Newtonian model

The model to be formulated should account for two related processes, the sinking in the vertical direction and the oscillations perpendicular to it. But as was shown above, the No- ring intruders have only small oscillations that end fast, being the overall sinking process almost vertical. Regarding the Ring intruders, though they strongly oscillate, the tilting process ends first, so we will consider only the equation controlling the vertical sinking, figuring out how the tilting angle affects the sinking dynamics.

In order to formulate a model to describe analytically the sinking process, let us consider the forces acting on the cylinder. As soon as the shaking starts, if the dimensionless acceleration is above threshold, the upper part of the granular bed is fluidized, and the intruder sinks.

Let us assume that the cylinder just sinks vertically, and let us name the vertical downward axis as $z$. The force balance on the intruder can be written as

$$m\vec{a} = m\vec{g} + \int (-P)\hat{n}dS + \int \sigma_s \cdot \hat{n}dS \quad (A.1)$$

where $P$ is the pressure, $\sigma_s$ the shear stress tensor, $\hat{n}$ is the vector normal to the intruder's surface, and the integrals run over the boundary of the intruder that is inside the granular material. Assuming a hydrostatic pressure profile, we can write:

$$P = \int_{0}^{h} \rho(z')g dz' \quad (A.2)$$
where $h$, as previously, is the depth reached by the cylinder below the surface of the granular medium. In Eq. (A.2) we have made explicit that the density of the material varies with depth. Let us assume that it varies as a power law between zero and the density of the solid layer, $\rho_{sl}$, that is reached at a depth $h_f$:

$$ \rho(z') = \rho_{sl} \left(\frac{z'}{h_f}\right)^p \qquad (A.3) $$

where $p \in [0, 1]$. The selection of the value of $p$ is discussed below (see also Appendix 2).

By combining (A.3) and (A.2) and integrating, we find the hydrostatic buoyancy force acting on the cylinder with a length $h$ under the (average) level of the granular bed, as:

$$ \int (-P) \hat{n} dS = - \frac{\rho_{sl} S g}{(p+1)h_f^p} h^{p+1} \hat{h} \qquad (A.4) $$

where $S$ is the characteristic area of the intruder cross section, and $\hat{h}$ is a unit vector pointing downwards. It is easy to see that the buoyancy force depends on the volume submerged into the granular medium.

Neglecting the inertial forces, which according to our simulations is typically two orders of magnitude smaller than the contact forces, the shear stress component goes as

$$ \int \sigma_s \cdot \hat{n} dS = -D\gamma v \hat{h} \qquad (A.5) $$

where $\gamma$ has the dimensions of a viscosity, $D$ is the characteristic size of the cross section of the intruder and $v$ is its sinking speed [58, 59]. By substituting Eq. (A.4) and Eq. (A.5) into Eq. (A.1), and only recovering the modular values, we get:

$$ m \frac{d^2 h}{dt^2} + D\gamma \frac{dh}{dt} + \frac{\rho_{cs} S g}{(p+1)h_f^p} h^{p+1} = mg \qquad (A.6) $$

Before solving Eq. (A.6) we will assume that the sink velocity is constant, which follows quite well the behavior during the fast sink regime, as seen in Fig. 15 (i.e., we neglect the inertial term). So,

$$ \frac{dh}{dt} + \frac{\rho_{sl} S g}{D\gamma(p+1)h_f^p} h^{p+1} = \frac{mg}{D\gamma} \qquad (A.7) $$

which can be written as
<figure><img src="image_17.png" /><figcaption>FIG. 15. (color online) Time dependence of sinking depth for the No-ring cylinder from experiment, compared with that determined from Eqs. (A.9, A.10). The inset shows the solutions of Eq. (A.6) for different values of p (see text).</figcaption></figure>

$$ \frac{dh}{dt} + ah^{p+1} = b \qquad (A.8) $$

The definitions of a and b are easily deduced by comparing Eqs. (A.7) and (A.8).

Equation (A.8) has analytical solutions if $p = 0$ or $p = 1$, which correspond to the extreme cases of constant density and a linear density profile with depth, respectively. The solutions are

$$ h(t) = \frac{b}{a}(1 - e^{-at}) \qquad (A.9) $$

if $p = 0$, and

$$ h(t) = \sqrt{\frac{b}{a}} \tanh(\sqrt{ab}t) \qquad (A.10) $$

if $p = 1$.

It is easy to see that both expressions correspond to an exponential growth that saturates.

Figure 15 shows the experimental results (continuous line) obtained for a dimensionless acceleration $\Gamma \simeq 0.3$. It is possible to see in more detail the initial fast sinking process, followed by the slow creep. Fig. 15 also shows the fitting of equations (A.9, A.10) to experimental data. Both solutions reproduces well the main features of the sinking process.

It is almost impossible to determine experimentally the exact density profile. But we do not need to know it in order to validate our model, if we use the following rationale.
Firstly, we fit Eqs. (A.9, A.10) to the experimental data and obtain the values of a, b that correspond to $p = 0$ ($a(p = 0), b(p = 0)$) and $p = 1$ ($a(p = 1), b(p = 1)$). Let us assume that a and b vary linearly with p between the extremes values which were obtained from the fitting process. For an intermediate value of p (say, $p_1$) we can calculate the corresponding values of $a(p_1)$ and $b(p_1)$. With them, we can in turn determine the constants of Eq. (A.6). Then, we solve this equation numerically. This procedure is repeated for values of p between 0 and 1, with a step of 0.1.

The inset in Fig. 15 shows some of the numerical solutions for the values of p in the legend. The main conclusion is that the density profile has small influence on the first (and most important) part of the sinking process. Of course, the final depth is influenced by the value of p, but due to experimental uncertainties, it is almost impossible to choose any particular value.

Let us now study the influence of the values of p in the quality of the fit of the solution of Eq. (A.6) to the experimental data. For doing this we notice that the values of a and b in Eq. (A.9) can be easily obtained from the experiments. Considering Eq. (A.8) in the first moments of motion, as h is small, $h'(t) \simeq b$, so b can be evaluated as the initial slope. As at large times $h(t) \sim h_{eq}$ ($h_{eq} = \Delta h$ if $h(0) = 0$) then $a = b/h_{eq}^{p+1}$. Then solving Eq. (A.6) for a given value of b, p, a(p) and naming the result $h_{mod}$, the best value of p arises from the minimization:

$$p_{opt} = \arg\min \sum_{i=1}^{N} [h_{mod}(t_i, p) - h_{exp}(t_i)]^2 \quad (A.11)$$

where $h_{exp}(t)$ are the experimental values of h.

The result for $\Gamma \le 1.0$ is indifferent to p: the fit is equally good no matter which is the value of $p \in [0, 1]$. For $\Gamma = 1.24$ there are differences in the quality of the fits for various values of p, but Eq. (A.11) gives a minimum for $p = 0$, so, we will use $p = 0$ in the following (Appendix 2 supports the selection of p from the simulations). Then, Eq. (A.6) becomes:

$$m \frac{d^2 h}{dt^2} + D \gamma \frac{dh}{dt} + \frac{\rho_{sl} S g}{h_f} h = mg \quad (A.12)$$

that can be taken as the simplest equation of motion describing the vertical sink dynamics of our cylinders. It is worth noticing that Eq. (A.12) reproduces quite closely the results reported in Fig. 15, and can be used to qualitatively describe the vertical sinking of Ring-
cylinders while tilting, as we will see below. It is possible to demonstrate that Eq. (A.12),
developed for a granular bed fluidized by shaking, is closely related with that proposed in [30]
to describe the penetration of an intruder into ultra-light granular material that eventually
behaves like a fluid medium even in the absence of shaking.

In order to understand the influence of the tilting dynamics in the sinking process, it is
useful to note that, when applying Eq. (A.12) to a tilted cylinder, the values of both *D* and *S* change. The reason is that when we calculate the surface integrals, the result will be
proportional to the cylinder's immersed volume. As the cylinder tilts, both the immersed
surface and linear dimensions increase more than in the case of sinking without tilting, so the
drag force is bigger in the former case. Considering, for instance, the situation represented
in the lower row of Fig. 2, when the cylinder sinks a distance *Δh*, the surface and linear
dimensions increase as the inverse of cos *θ* (of course, other intruder geometries may follow
different laws).

To test it, let us assume a simplified model: the increase factor of *S* and *D* is proportional
to the characteristic size of the cross section of the cylinder projected on the horizontal plane,
i.e., it is proportional to the inverse of cos θ. Then, instead of *D* and *S*, we will solve Eq.
(A.12) using *D*/cos θ(t) and *S*/cos θ(t), where θ(t) is a function that grows from zero to the
maximum angle θmax reached by the cylinder, mimicking Fig. 6(c), i. e. with the same time
constant.

The consequences can be seen in Fig. 16. While in the beginning the sinking process
in all situations occurs with the same dynamics, as the cylinder approaches the final angle,
the behavior changes, being the final depth larger for the situations corresponding to low
tilting.

The upper curve, calculated for θ = 0 coincides with the upper curve in the inset of Fig.
15 (calculated for p = 0). Subsequent curves are calculated for values of θmax varying in
steps of π/15, the lowermost curve corresponds to θmax = π/3. As the inclination of the
cylinder increases, both the buoyancy and the viscous drag do. The effect of these factors
on the sinking process of Ring cylinders was already noted in Fig. 6(b): an immediate
consequence is the decrease of the sinking depth (for a given Γ), compared with that of the
No-ring ones, which can be easily observed in the experiments. From the inset it is possible
to deduce that, for the larger angles, a small decrease in the depth is observed.

In spite of the simplifications assumed, it is worth noting that one of the basic differences
<figure><img src="image_18.png" /></figure>

FIG. 16. (color online) Time dependence of sinking depth as calculated solving numerically Eq. (A.12) considering the variation of *S* and *D* provoked by tilting (see text). Upper curve is for *θ*<sub>max</sub> = 0 while the lower one is for *θ*<sub>max</sub> = π/3. Between them, *θ*<sub>max</sub> varies in steps of π/15. The inset shows the last three seconds.

between Fig. 7 (a) and (b) –for a given Γ the No-ring cylinders sink deeper than the Ring ones– could be qualitatively described by our model.

Finally, there is another element that was neither considered by us: as the container
shakes horizontally, it produces a horizontal drag that changes periodically its direction.
According to [60], it creates an additional lift force, and also a dependence of the drag force
with depth, which, of course, must influence the detailed penetration dynamics of the Ring
cylinders. The results of Li et al [40] also support these ideas.

2. Density profile of the granular medium

At each time step the simulation box is subdivided into a fixed number of rectangles for which the density is calculated as $\rho = m/V$, where $m$ is the sum of the masses of all the particles of the granular medium within each rectangle and $V$ is the volume of the rectangle calculated as its area multiplied by the average diameter of the particles of the medium inside it.

The density profile for each dimensionless acceleration $\Gamma$ is obtained by calculating the average of the density profile at each time step. The latter is obtained by averaging the profiles of a set of rectangle columns near the intruder. These columns are chosen avoiding the presence of the intruder in them in order not to affect the density profile since only the particles that compose the granular medium are taken into account in its calculation.

Figure 17 shows the density profiles as a function of depth for the different values of
<figure><img src="image_19.png" /></figure>

FIG. 17. (color online) Density profiles for the different values of dimensionless acceleration $\Gamma$ obtained from the simulations. The background shows the subdivision of the simulation box where the color of each rectangle goes sequentially from blue to red and represents its density. The intruder is represented in light gray.

$\Gamma$. As can be seen in the region of the graph enclosed by dashed lines, the density of the medium saturates rapidly with increasing depth for all dimensionless accelerations. This fact supports from the simulations the use of $p = 0$ in expression A.3 of Appendix 1.

[1] J. D. Bray, J. P. Stewart, M. B. Baturay, T. Durgunoglu, A. Önalp, R. B. Sancio, J. P. Stewart, D. Ural, A. Ansal, J. B. Bardet, and et. al., Damage patterns and foundation performance in adapazari, Earthquake Spectra **16**, 163 (2000).

[2] R. B. Sancio, J. D. Bray, J. P. Stewart, T. L. Youd, H. T. Durgunoglu, A. Önalp, R. B. Seed, C. Christensen, M. B. Baturay, and T. Karadayılar, Correlation between ground failure and soil conditions in adapazari, turkey, Soil Dynamics and Earthquake Engineering **22**, 1093 (2002).

[3] R. B. Sancio, J. D. Bray, T. Durgunoglu, and A. Önalp, Performance of buildings over liquefiable ground in adapazari, turkey, in Proc., 13th World Conf. on Earthquake Engineering (Canadian Association for Earthquake Engineering Vancouver, Canada, 2004).

[4] N. Ambraseys, Engineering seismology: Part i, Earthquake Engineering & Structural Dynamics **17**, 1 (1988).

[5] S. Ben-Zeev, L. Goren, S. Parez, R. Toussaint, C. Clément, and E. Aharonov, The combined effect of buoyancy and excess pore pressure in facilitating soil liquefaction, in Poromechanics
VI (2017) pp. 107–116.

[6] S. Ben-Zeev, E. Aharonov, R. Toussaint, S. Parez, and L. Goren, Compaction front and pore fluid pressurization in horizontally shaken drained granular layers, Physical Review Fluids 5, 054301 (2020).

[7] K. Vanneste, M. Meghraoui, and T. Camelbeeck, Late quaternary earthquake-related soft-sediment deformation along the belgian portion of the feldbiss fault, lower rhine graben system, Tectonophysics 309, 57 (1999).

[8] N. R. Council, *Liquefaction of Soils During Earthquakes* (The National Academies Press, Washington, DC, 1985).

[9] J. Berrill and R. Davis, Energy dissipation and seismic liquefaction of sands: revised model, Soils and foundations 25, 106 (1985).

[10] E. Clement and J. Rajchenbach, Fluidization of a bidimensional powder, EPL (Europhysics Letters) 16, 133 (1991).

[11] C. Wang and M. Manga, *Earthquakes and Water* (Springer-Verlag Berlin Heidelberg, 2009).

[12] H. M. Jaeger, S. R. Nagel, and R. P. Behringer, Granular solids, liquids, and gases, Reviews of modern physics 68, 1259 (1996).

[13] T. Shinbrot and F. J. Muzzio, Reverse buoyancy in shaken granular beds, Phys. Rev. Lett. 81, 4365 (1998).

[14] E. Altshuler, O. Ramos, E. Martínez, A. J. Batista-Leyva, A. Rivera, and K. E. Bassler, Sandpile formation by revolving rivers, Phys. Rev. Lett. 91, 014501 (2003).

[15] I. S. Aranson and L. S. Tsimring, Patterns and collective behavior in granular media: Theoretical concepts, Reviews of modern physics 78, 641 (2006).

[16] B. Andreotti, Y. Forterre, and O. Pouliquen, *Granular media: between fluid and solid* (Cambridge University Press, 2013).

[17] E. Altshuler, R. Toussaint, E. Martínez, O. Sotolongo-Costa, J. Schmittbuhl, and K. J. Måløy, Revolving rivers in sandpiles: From continuous to intermittent flows, Phys. Rev. E 77, 031305 (2008).

[18] J. Boudet, Y. Amarouchene, B. Bonnier, and H. Kellay, The granular jump, Journal of Fluid Mechanics 572, 413 (2007).

[19] J. L. Vinningland, Ø. Johnsen, E. G. Flekkøy, R. Toussaint, and K. J. Måløy, Granular rayleigh-taylor instability: Experiments and simulations, Phys. Rev. Lett. 99, 048001 (2007).
[20] Ø. Johnsen, R. Toussaint, K. J. Måløy, and E. G. Flekkøy, Pattern formation during air injection into granular materials confined in a circular hele-shaw cell, Phys. Rev. E **74**, 011301 (2006).

[21] Ø. Johnsen, R. Toussaint, K. J. Måløy, E. G. Flekkøy, and J. Schmittbuhl, Phys. Rev. E **77**, 011301 (2008).

[22] M. J. Niebling, E. G. Flekkøy, K. J. Måløy, and R. Toussaint, Sedimentation instabilities: impact of the fluid compressibility and viscosity, Phys. Rev. E **82**, 051302 (2010).

[23] M. J. Niebling, E. G. Flekkøy, K. J. Måløy, and R. Toussaint, Mixing of a granular layer falling through a fluid, Phys. Rev. E **82**, 011301 (2010).

[24] M. J. Niebling, R. Toussaint, E. G. Flekkøy, and K. J. Måløy, Dynamic aerofracture of dense granular packings, Phys. Rev. E **86**, 061315 (2012).

[25] S. Turkaya, R. Toussaint, F. K. Eriksen, M. Zecevic, G. Daniel, E. G. Flekkøy, and K. J. Måløy, Bridging aero-fracture evolution with the characteristics of the acoustic emissions in a porous medium, Frontiers in Physics **3**, 70 (2015).

[26] J. Uehara, M. Ambroso, R. Ojha, and D. J. Durian, Low-speed impact craters in loose granular media, Phys. Rev. Lett. **90**, 194301 (2003).

[27] J. Boudet, Y. Amarouchene, and H. Kellay, Dynamics of impact cratering in shallow sand layers, Phys. Rev. Lett. **96**, 158001 (2006).

[28] H. Katsuragi and D. J. Durian, Unified force law for granular impact cratering, Nature physics **3**, 420 (2007).

[29] D. I. Goldman and P. Umbanhowar, Scaling and dynamics of sphere and disk impact into granular media, Phys. Rev. E **77**, 021308 (2008).

[30] F. Pacheco-Vázquez, G. Caballero-Robledo, J. Solano-Altamirano, E. Altshuler, A. Batista-Leyva, and J. Ruiz-Suárez, Infinite penetration of a projectile into a granular medium, Phys. Rev. Lett. **106**, 218001 (2011).

[31] L. Kondic, X. Fang, W. Losert, C. O'Hern, and R. Behringer, Microstructure evolution during impact on granular matter, Phys. Rev. E **85**, 011305 (2012).

[32] H. Torres, A. González, G. Sánchez-Colina, J. Drake, and E. Altshuler, Impact dynamics in “hard” and “soft” granular matter, Revista Cubana de Física **29**, 1 (2012).

[33] A. H. Clark and R. P. Behringer, Granular impact model as an energy-depth relation, EPL (Europhysics Letters) **101**, 64001 (2013).
[34] J. Ruiz-Suárez, Penetration of projectiles into granular targets, Reports on Progress in Physics **76**, 066601 (2013).

[35] T. A. Brzinski, P. Mayor, and D. J. Durian, Depth-dependent resistance of granular media to vertical penetration, Phys. Rev. Lett. **111**, 168002 (2013).

[36] E. Altshuler, H. Torres, A. González-Pita, G. Sánchez-Colina, C. Pérez-Penichet, S. Wait-ukaitis, and R. Hidalgo, Settling into dry granular media in different gravities, Geophysical Research Letters **41**, 3032 (2014).

[37] S. Joubaud, T. Homan, Y. Gasteuil, D. Lohse, and D. van der Meer, Forces encountered by a sphere during impact into sand, Phys. Rev. E **90**, 060201 (2014).

[38] R. Harich, T. Darnige, E. Kolb, and E. Clément, Intruder mobility in a vibrated granular packing, EPL (Europhysics Letters) **96**, 54003 (2011).

[39] É. Kolb, J. Cviklinski, J. Lanuza, P. Claudin, and É. Clément, Reorganization of a dense granular assembly: The unjamming response function, Phys. Rev. E **69**, 031306 (2004).

[40] C. Li, T. Zhang, and D. I. Goldman, A terradynamics of legged locomotion on granular media, science **339**, 1408 (2013).

[41] J. Aguilar and D. I. Goldman, Robophysical study of jumping dynamics on granular media, Nature Physics **12**, 278 (2016).

[42] W. Kang, Y. Feng, C. Liu, and R. Blumenfeld, Archimedes' law explains penetration of solids into granular media, Nature communications **9**, 1 (2018).

[43] V. L. Díaz-Melián, A. Serrano-Muñoz, M. Espinosa, L. Alonso-Llanes, G. Viera-López, and E. Altshuler, Rolling away from the wall into granular matter, Phys. Rev. Lett. **125**, 078002 (2020).

[44] M. Medved, H. M. Jaeger, and S. R. Nagel, Modes of response in horizontally vibrated granular matter, Europhysics Letters **52**, 66 (2000).

[45] C. A. Kruelle, Physics of granular matter: pattern formation and applications, Rev. Adv. Mater. Sci **20**, 113 (2009).

[46] S. Tennakoon, L. Kondic, and R. Behringer, Onset of flow in a horizontally vibrated granular bed: Convection by horizontal shearing, EPL (Europhysics Letters) **45**, 470 (1999).

[47] L. Liu and R. Dobry, Seismic response of shallow foundation on liquefiable sand, Journal of geotechnical and geoenvironmental engineering **123**, 557 (1997).

[48] G. Sánchez-Colina, L. Alonso-Llanes, E. Martinez, A. Batista-Leyva, C. Clement, C. Fliedner,
R. Toussaint, and E. Altshuler, Note: “lock-in accelerometry” to follow sink dynamics in shaken granular matter, Review of scientific instruments **85**, 126101 (2014).

[49] C. Clément, R. Toussaint, and E. Aharanov, Shake and sink: liquefaction without pressurization, arXiv preprint arXiv:1802.04391 (2018).

[50] C. Clément, R. Toussaint, M. Stojanova, and E. Aharonov, Sinking during earthquakes: Critical acceleration criteria control drained soil liquefaction, Phys. Rev. E **97**, 022905 (2018).

[51] The maximum oscillation amplitude that the shaker can generate is 13 mm. For a frequency of 5.0 Hz, the maximum possible acceleration is 12.83 m/s².

[52] U.S. Geological Survey, Shakemap scientific background, https://earthquake.usgs.gov/data/shakemap/background.php, accessed: 2021-11-5.

[53] H. Goto, Y. Kaneko, J. Young, H. Avery, and L. Damiano, Extreme accelerations during earthquakes caused by elastic flapping effect, Scientific Reports **9**, 1117 (2019).

[54] P. A. Cundall and O. D. Strack, A discrete numerical model for granular assemblies, Geotechnique **29**, 47 (1979).

[55] S. Parez, E. Aharonov, and R. Toussaint, Unsteady granular flows down an inclined plane, Phys. Rev. E **93**, 042902 (2016).

[56] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, Lammmps - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Computer Physics Communications **271**, 108171 (2022).

[57] R. Toussaint, C. Clément, C. Fliedner, M. Stojanova, E. Aharonov, G. Sanchez, E. Altshuler, A. Batista, and E. Grude Flekkoy, Soil liquefaction: how the granular medium evolves, macroscale and microscale study, in *EGU General Assembly Conference Abstracts* (2014) p. 12539.

[58] J. R. de Bruyn and A. M. Walsh, Penetration of spheres into loose granular media, Canadian Journal of Physics **82**, 439 (2004).

[59] M. Hou, Z. Peng, R. Liu, Y. Wu, Y. Tian, K. Lu, and C. Chan, Projectile impact and penetration in loose granular bed, Science and Technology of Advanced Materials **6**, 855 (2005).

[60] X. Zhang, D. Sheng, G. P. Kouretzis, K. Krabbenhoft, and S. W. Sloan, Numerical investiga-
tion of the cylinder movement in granular matter, Phys. Rev. E **91**, 022204 (2015).