Induced and endogenous acoustic oscillations in granular faults

TL;DR

This paper reviews how external perturbations induce acoustic oscillations in granular fault models, affecting their frictional properties and potentially explaining seismic phenomena like fault weakening and earthquake triggering.

Contribution

It provides a comprehensive review of recent findings on the response of granular systems to perturbations, highlighting non-monotonic weakening effects and re-entrant transitions.

Findings

External perturbations can induce acoustic oscillations in granular faults.

Frictional weakening may not increase monotonically with perturbation frequency.

Re-entrant transitions between fluidized and frictional states are observed.

Abstract

The frictional properties of disordered systems are affected by external perturbations. These perturbations usually weaken the system by reducing the macroscopic friction coefficient. This friction reduction is of particular interest in the case of disordered systems composed of granular particles confined between two plates, as this is a simple model of seismic fault. Indeed, in the geophysical context frictional weakening could explain the unexpected weakness of some faults, as well as earthquake remote triggering. In this manuscript we review recent results concerning the response of confined granular systems to external perturbations, considering the different mechanisms by which the perturbation could weaken a system, the relevance of the frictional reduction to earthquakes, as well as discussing the intriguing scenario whereby the weakening is not monotonic in the perturbation frequency, so that a re-entrant transition is observed, as the system first enters a fluidized state and then returns to a frictional state.

Equations

$${\mathbf{F}}_{ij}={\mathbf{F}}_{ij}^n+{\mathbf{F}}_{ij}^t=\left[k_n\delta {\mathbf{n}}_{ij}-m_{\mathrm{eff}}{\gamma }_n{\mathbf{v}}_n)\right]+\left[-(k_t\mathbf{t}+m_{\mathrm{eff}}{\gamma }_t{\mathbf{v}}_t)\right].$$

$${\chi }_{\alpha }(t)=\frac{1}{\alpha P_0}\underset{\tau \to \infty }{\lim }{\left[\frac{1}{N}\sum _i{\left({\mathbf{r}}_{\mathbf{i}}^{\alpha }(\mathbf{t}+\tau )-{\mathbf{r}}_{\mathbf{i}}^0(\mathbf{t}+\tau )\right)}^2\right]}^{1/2},$$

$$Z_p(x,t)=A\sin (2\pi ft)+A_x\sin (k_{x}x),$$

$$m\ddot{z}=k_n(Z_p-z)\Theta (Z_p-z)-{\gamma }_n(\dot{z}-\dot{Z_p})\Theta (Z_p-z)-mg,$$

$$m\ddot{x}=-k_d(x-V_{\text{d}}t)-{\gamma }_t\dot{x}-k_t\mathcal{F}(\dot{x},Z_p-z),$$

$$m{\ddot{z}}^{\prime }$$

$$m{\ddot{x}}^{\prime }$$

$$m{\ddot{z}}^{\prime }$$

$$m{\ddot{x}}^{\prime }$$

$$\varphi =\frac{⟨(\dot{x}-V_d{)}^2{⟩}_{t_a}}{V_d^2},$$

$$Z(t)=A\sin (2\pi ft),$$

$$f_1=\frac{\eta }{2}\left(\frac{M_p}{M}+\sqrt{\frac{M_p^2}{M^2}+4\frac{F_N}{MA{\eta }^2}}\right).$$

$$f_2=\sqrt{2\pi \frac{F_N}{MA}}.$$

$$2\pi f_1=\sqrt{g/A},$$

$${\omega }_{AF}=\frac{\pi }{W}\sqrt{\frac{k_n}{6\rho d}}.$$

$$P_{+}(t_p,t)=+\frac{\alpha P_0}{2}\left[1-\sin (\frac{\pi }{2}+\omega (t_p-t))\right],$$

Full text

I Introduction

Seismic occurrence is an intermittent phenomenon which is mainly controlled by macroscopic friction. The observation of earthquake triggering induced by relatively small changes in the stress hill93 ; kilb2000 ; gomberg2004 ; gomberg2005 suggests that friction weakens when faults are subject to external perturbations. This phenomenon can be attributed to the presence of crushed and ground-up rocks produced during past sliding events, known as fault gouge. This can be treated as a granular material which can act either as a lubricant or, via the formation of strong force chains, can inhibit the relative slip between the fault walls. The two different behaviours correspond to the double nature of granular materials which can be found either in a fluid like ’unjammed’ state or in an amorphous solid ’jammed’ state. An earthquake can therefore be interpreted as an unjamming transition from a jammed state, in which the gouge resists the existing stresses, to a flowing one. The transition from jammed to unjammed states plays a central role not only in earthquake triggering but also in the physics of avalanche dynamics herrmann as well as in the manufacturing process in material, alimentary and pharmaceutical industries coussot . In the latter case understanding the role of external perturbations takes on a great practical relevance, due to the possibility of controlling the mechanisms which reduce friction or enhance fluidization in different conditions persson ; krim ; rmpvanossi .

External perturbations are well known to be able to induce the transition from a jammed to a flowing state in confined granular systems experiencing a high normal stress and a small shear stress, and which would therefore be jammed in the absence of perturbations. This perturbation induced unjamming transition might actually explain the shaking induced fault weakening  UvH05 ; urb2004 ; JJ05 ; DB05 ; DB06 ; JSKGM08 ; BJM08 ; BJJ08 ; CVVZ09 ; krim ; GLPC12 ; CRBUF11 ; GDGJMC11 ; DWDDH11 ; VEBBJ12 ; CVVZ12 ; GFGDJMC13 ; XHM13 ; maksephyschemestry ; LELC15 ; ROLJJ15 ; LGV15 ; GSdALP15 ; corwin2015 ; LEC16 ; WDH16 ; dGW17 . Interestingly, external perturbations can fluidize granular systems through two distinct mechanisms. On the one side, external perturbations in given amplitude and frequency ranges can produce the detachment between the confining plates and the fault gouge. When this occurs the gouge does not oppose to the movement of the plates, which therefore flow. This behaviour can be understood, at least at the qualitative level, ignoring internal degrees of freedom of individual grains and treating the granular material as a single block, i.e. describing the fault as a block subject to a normal and a shear stress, resting on a oscillating plate. This simplified description permits analytical results which can be used to explain the non-monotonic behaviour of the friction as function of the vibration frequency CVVZ09 ; CVVZ12 ; GLPC12 ; prl2018 .

On the other side, external perturbations can also induce fault weakening in the absence of detachment. In this case the single block description is no longer useful and it is necessary to take into account the reorganization of the contact force network among grains induced by sound propagation. This reorganization can explain the experimentally observed modulus softening JJ05 and hysteretic behaviour identified for sufficiently large acoustic perturbations BJM08 ; BJJ08 . Among several mechanisms proposed to explain fault weakening induced by transient waves, in this Chapter we focus on the hypothesis of Acoustic Fluidization (AF). The AF concept was introduced by Melosh in 1979 to explain the transition from simple to complex craters on the Moon Mel79 and also the low coefficients of friction observed in large-volume rock avalanches Mel86 . According to Melosh Mel79 ; Mel96 dynamic fault weakening can be attributed to the activation of short-wavelength vibrations in the fault core. These could generate stationary oscillations at a characteristic frequency $\omega_{AF}$ counteracting the applied stress and eventually “lubricating” the system. The AF mechanism is therefore expected to strongly depend on the value of the confining stress as well as on the thickness of the granular medium which affects $\omega_{AF}$. This scenario is consistent with laboratory experiments on confined grains which have documented XHM13 a clear transition, induced by an external shaking, from elastic like to fluid like behaviour. These experiments, indeed, find that the unjamming transition is dramatically controlled by the thickness of the sample and the striker impact velocity, which modifies the confining pressure, in agreement with the AF scenario. The AF mechanisms has been also proposed to explain the change in density of flowing materials VEBBJ12 as well as the velocity-weakening of macroscopic friction dGW17 .

In this chapter we consider the different fludization mechanisms in numerical models of a granular fault. At variance with experiments, numerical simulations allow us to follow the trajectory of single grains obtaining insights not available in experiments. For example, numerical studies FGGJMC13 have shown that vibrations affect both weak and strong contacts among grains. However when the perturbation is switched off the weak contacts essentially returns to their initial state while the strong ones remain weakened. The weakening of strong contacts can be responsible of the clock advance of large slip events induced by the shaking. Furthermore, in numerical simulations it is also possible to explore the endogenous activation of oscillation at the frequency $\omega_{AF}$. This is the mechanism proposed by Melosh to justify the unexpectedly small ratio between shear and normal stress in real fault systems. The idea is that the elastic strain accumulated inside the fault is able to excite short wave-lengths which in turn can promote the activation of additional vibrations inside the fault core leading to self-fluidization.

In the following, we first present the details of the considered model for the granular faults, and briefly discuss some of its static features. Then in Sec. III we focus on the conditions under which the external perturbation induces the detachment between grains and the confining plate. In particular we use the results of a single block dynamics to explain the outcomes of experiments measuring the angular velocity $\omega$ of a vane coupled to an external motor and embedded inside a confined granular material. In Sec. IV we focus on the AF mechanism, considering the role of of external perturbations not able to detach the confining plates, and on the endogenous activation of acoustic vibrations at the onset of stick-slip instabilities. Conclusions are drawn in the last Section.

II A model for granular faults

Geophysical phenomena as earthquakes occur over very large length scales and involve high values of energy, implying that it is clearly impossible to reproduce them at the laboratory scale. Nevertheless, earthquakes are characterised by scale-invariant laws, as the ubiquitous Gutenberg-Richter law for the magnitude distribution, which indicates the absence of characteristic scales. This implies that the experimental and numerical investigation of model seismic faults might give access to the relevant physical properties of real faults. Many experiments and simulations of model seismic faults are therefore present in recent literature.

II.1 Numerical model

The numerical model of seismic fault presented here focuses on the fault core dynamics, where most of the shear displacement occurs. Results are obtained under the assumption that the region outside the core is a rigid body, so that the fault of width $W$ can be modelled as two parallel rigid plates of area $L_{x}L_{y}$ confined by a normal stress $P_{0}$ and subject to a shear stress $\sigma$ (fig.1). See Giacco2014b ; 3ddegriffa2014 for results considering elastic confining boundaries. The fault gouge, that in real faults consists of rocks produced in past wearing events, is modelled as a collection of frictional granular particles over a width $W$. Particles are mono-dispersed spheres of mass $m$ and diameter $d$. This choice does not lead to crystallization because of the rough confining plates, however poly-dispersed particles would represent a more realistic modelization of real gouges, numerically more demanding. In real faults the shear stress $\sigma$ slowly increases as a consequence of the convective motion of the upper mantle, and quickly decreases when a fault slips. In simulations, an analogous scenario is realised by applying the stress through a spring mechanism. One end of the spring is attached to the plate, and thus moves with the plate velocity $v_{p}(t)$, while the other end is driven with a constant velocity $V$. Assuming the spring to have zero equilibrium length and zero length at time $t=0$, the force exerted by the spring at time $t$ is $F_{t}=k_{d}(Vt-x(t))$, where $x(t)=\int_{0}^{t}v_{p}(t^{\prime})dt^{\prime}$ and $k_{d}$ is the spring stiffness. Then the shear stress results to be $\sigma_{xz}=\sigma=F_{t}/L_{x}L_{y}$. The stick-slip dynamics is recovered when the time over which the stress increases, which depends on the shear rate $k_{d}V/L_{x}L_{y}$, is much longer than the duration of a slip event, which depends on the pressure $P_{0}$ as well as on the dissipative properties of the system. In numerical models, the confining plates are rough rigid objects, generally realised as a dense disordered assembly of granular particles. The plate is rigid since all relative distances between the particles in the plate are kept constant and particles in the plate are not allowed to rotate.

The interaction force between two granular particles has a normal and a tangential component. To describe the interaction, two particles $i$ and $j$ are considered, with radii $R_{i}$ and $R_{j}$, in position $r_{i}$ and $r_{j}$, linear velocities $v_{i}$ and $v_{j}$ and angular velocities $\omega_{i}$ and $\omega_{j}$. The particles interact only when in physical contact, i.e. when $|r_{i}-r_{j}|<d_{ij}=R_{i}+R_{j}$, i.e. when the overlap $\delta=d_{ij}-|r_{i}-r_{j}|>0$. We use the Harmonic spring-dashpot model silbert ,

[TABLE]

The first term of Eq. (1) corresponds to the normal interaction. Here $k_{n}$ is the stiffness of the particles, whereas $\gamma_{n}$ controls the dissipation and $m_{\rm eff}$ is the effective mass of spheres with mass $m_{i}$, $m_{j}$. In this model, the normal collision between two grains dissipates a constant fraction $(1-e^{2})$ of the overall kinetic energy, where $e(k_{n},\gamma_{n},m_{eff})$ is the restitution coefficient. The second term of Eq. (1) corresponds to the tangential interaction force, where $k_{t}$ is a tangential stiffness, ${\bf t}$ the tangential shear displacement, and $\gamma_{t}$ a damping parameter. The tangential shear displacement ${\bf t}$ is the integral of the relative tangential velocity ${\bf v}_{t_{ij}}$ at the point of contact, where ${\bf v}_{t_{ij}}={\bf v}_{ij}-{\bf v}_{n_{ij}}-1/2({\bf\omega}_{i}+{\bf\omega}_{j})\times{\bf r}_{ij}$ depends on the angular velocities of the particles. Here the relative velocity is ${\bf v}_{ij}={\bf v}_{i}-{\bf v}_{j}$ and ${\bf r}_{ij}={\bf r}_{i}-{\bf r}_{j}$. The tangential force acting on a particle contributes to the overall torque acting on it, which induces its rotational motion. The presence of a finite Coulomb’s friction coefficient $\mu$ is implemented by introducing an upper bound for the tangential force, which is capped at $\mu|F_{n}|$. This is done by appropriately rescaling ${\bf t}$, thus mimicking the slipping of the contact. We measure the mass in units of $m$, the lengths in units of $d$ and time in units of $\sqrt{m/k_{d}}$. Unless specified, results are obtained for a restitution coefficient $e=0.8$, a coefficient of static friction $\mu=0.2$, the confining pressure is $P_{0}={k_{d}/d}$, $k_{n}=2\cdot 10^{3}\,k_{d}$, $V=0.01\,d\sqrt{m/k_{d}},\gamma_{t}=0$, the temporal integration step of the equations of motion is $5\cdot 10^{-3}\,\sqrt{m/k_{d}}$ and $W\simeq 10\,d$ as in Ref. prl2010 . The model reproduces the stick-slip dynamics with slipping event sizes distributed according to the Gutenberg-Richter law characterising real seismic occurrence, independently of model parameters. In particular, the dynamics consists of almost periodic large events, called slips, and of power-law distributed smaller events, called microslips (Fig. 7) epl2011 .

II.2 Micromechanics of failure

Given that the numerical model gives access to all possible quantities of interest, it is useful to investigate how does the system fail. There are two possibilities, both of them based on the consideration that, as the shear stress increases, the inter-particle forces slightly change, since the system is in mechanical equilibrium at all times. In the local scenario, failure occurs as one of the inter-particle contacts reaches the Coulomb threshold and starts slipping. This local failure then propagates leading to the macroscopic failure of the sample. In the global scenario, conversely, the shear stress increase leads to a deformation of the energy landscape of the system. Thus the energy minimum, in which the system is trapped when in a jammed configuration, might gradually evolve into a saddle. When the minimum becomes a saddle, the system fails. This second scenario is certainly at work in the absence of frictional forces, that are known not to influence the statistical properties of the model epl2011 . Formally, in the absence of friction this second scenario occurs when an eigenvalue of the dynamical matrix of the system vanishes. It has been shown prl2010 that also in the presence of frictional forces failure results from a global instability of the system. Indeed, the system failure time is the first time at which, if the shear stress is kept constant, at least one contact reaches its Coulomb threshold.

In frictional systems one cannot explicitly show that failure is a global process investigating the dynamical matrix of the system, due to the absence of a Hamiltonian. As an alternative, the global instability scenario suggests an analogy with a second order phase transition, where a minimum evolves into a saddle point. On approaching the failure point one therefore expects the system to become increasingly more sensitive to external perturbations. It is therefore interesting to investigate the evolution of the response of the system to external perturbations. A perturbation is considered as a decrease of the confining pressure $P_{0}$ by $\alpha P_{0}$ for a short time $\delta t_{pert}=0.1$, with $\alpha\ll 1$ or small enough to probe the linear response regime. In order to guarantee the separation of time scales between the external drive and the mechanical relaxation, the external drive of the spring is kept constant by fixing $V=0$ when performing this study. As a consequence of the applied perturbation, grains move. The response of the system is monitored by studying the susceptibility

[TABLE]

where $\bf{r}_{i}^{\alpha}$ and $\bf{r}_{i}^{0}$ are the asymptotic particle positions in the perturbed and unperturbed system, respectively. Here the $t$ dependence indicates solely the time at which the perturbation is applied since the susceptibility is evaluated asymptotically and is therefore a static quantity. In the unjammed phase the susceptibility is divergent since global position rearrangements occur, whereas in the jammed phase only limited regions give contribution to $\chi_{\alpha}$ which therefore provides a measure of the correlation length.

By monitoring the response of the system for different $\alpha$ close to slips and microslips a very different behaviour is detected prl2010 . Indeed, close to microslips the response of the system is linear in $\alpha$, namely $\chi_{\alpha}$ increases proportionally to $\alpha$ as the microslip time is approached and drops to zero at the subsequent jamming time. Interestingly, the $\chi_{\alpha}$ value at the onset of the microslip is proportional to the microslip size, indicating that the information on the slip size is contained in the global state of particle positions. Conversely, as the large slip occurrence time is approached, the susceptibility suddenly diverges, the sooner the larger $\alpha$. Therefore, in the limit of vanishing perturbation $\chi_{\alpha}$ diverges at the slip time. This divergence, reminiscent of the response of a system close to a critical point, indicates that global rearrangements of the particle positions occur leading to a correlation length of the order of the system size. In order to better simulate the mechanical conditions of seismic faults, the system should undergo periodic perturbations representing the passage of seismic waves triggered by earthquakes just occurred even far in space. The frequency dependence of the response will be discussed in the next Sections.

III Oscillations inducing the detachment between grains and the confined plates.

In this section we review some recent results on the response of granular packings to mechanical perturbations focusing on the role played by the vibration frequency on the frictional properties of the medium. In particular, we first discuss a single block model under vibration, to introduce the issue of friction reduction in a simple context. Then we consider a granular system similar to the gouge fault model introduced in Sec. II analysing the central role played by the vibration frequency. Finally we report some recent experimental results on shaken granular systems in a different setup, suitable for rheological studies, where a vane coupled to an external motor, immersed in the medium, is used to probe its viscosity properties.

III.1 A case study: The single block

Friction between sliding objects is largely affected by mechanical vibrations. To introduce the problem on simple grounds, a numerical study GLPC12 on a spring-block model in the presence of vertical vibrations is considered. This allows to bring to the fore some fundamental ingredients that influence the frictional properties of sliding solids, such as the external drive, the geometry of the surfaces over which the block moves and the kind of confining force. In particular, one observes non-trivial behaviour as a function of the vibration frequency, such as friction reduction, that will be analysed in more complex systems in Subsection III.2. The solid-on-solid model consists in a block of mass $m$ that is pulled along the $x$ direction by a spring of elastic constant $k_{d}$ driven at constant velocity $V_{d}$ and moves along a surface which is vibrated in the $z$ vertical direction according to

[TABLE]

where $Z_{p}$ is the vertical coordinate of the plate (see Fig. 2). One can identify three main variants: In model A, the surface of the substrate is flat ($A_{x}=0$) and the confining force is vertical; in model B, the surface has a sinusoidal shape ($A_{x}>0$) and the confining force is still vertical; in model C, the surface is sinusoidal ($A_{x}>0$) but the confining force is always normal to it. In particular, the equations of motion along the two directions $(x,z)$ with respect to a fixed reference system in model A are given by:

[TABLE]

[TABLE]

where $\Theta$ is the Heaviside step function, $k_{n}$ is the elastic constant, $\gamma_{n,t}$ are the viscoelastic constants, normal and tangential respectively, while the quantity $k_{t}\mathcal{F}$ is a frictional term that is present when block and surface are in contact. It is assumed to be proportional to the shear displacement over the contact time interval $(t,t_{0})$ between the block and the plate cundall and is given by $\mathcal{F}(\dot{x}(t),Z_{p}-z)=\int_{t_{0}}^{t}\dot{x}(t^{\prime})dt^{\prime}$. The Coulomb friction is taken into account through the condition $|k_{t}\mathcal{F}|<\mu_{s}N$, where $N$ is the normal load; when this condition is violated $\mathcal{F}$ is set to zero.

In model B and C we write the equations of motion in a frame of reference which moves along with the plate, with the horizontal axis tangential to the plate, see Fig. 2. In this frame the equations of motion for model B are

[TABLE]

where $\alpha=\arctan\left(\frac{\partial{Z_{p}}}{\partial{x}}\right)$ at fixed time, while model C is described by

[TABLE]

with $P_{l}$ a confining force which is always perpendicular to the surface.

Spring block models without vibrations present a sliding (fluid) phase and a stick-slip (solid) phase vasco and the presence of vibrations affects the transition between these two phases. The three models introduced above show different behaviour: The main properties of model B only depend on the driving velocity and are almost independent of the presence of vibrations, whereas models A and C show a transition from the stick-slip to the sliding phase for increasing frequency (or amplitude) of the oscillating substrate, occurring when its maximum acceleration overcomes the gravity acceleration $g$. Remarkably, in model C one observes that a further increase of the frequency leads to a second transition whereby the system re-enters the stick-slip phase. More specifically, in order to better identify the transition one can introduce an order parameter, defined as

[TABLE]

where the brackets $\langle\cdots\rangle_{t_{a}}$ indicate temporal averages over a period $t_{a}$. During the flowing phase, the block moves with the external drive velocity ($\dot{x}=V_{d}$) and therefore $\phi=0$. Conversely, in the stick-slip phase , $\phi$ takes a finite value that depends on the period $t_{a}$ (or the number of occurred slips during $t_{a}$). The different behaviours of models A and C are summarized by the phase diagrams in the parameter space $(f,A)$, obtained from numerical simulations in GLPC12 and reported in Fig. 3. Model A is characterized by the line $A=g(2\pi f)^{-2}$, separating a stick-slip region and a sliding region that originates from the detachment condition. The phase diagram of model C is more complex, featuring a second transition from the stick-slip to the flowing phase at high frequencies. This friction recovery transition is related to the balance condition between inertial and dissipative forces GLPC12 . The results obtained for these single block models indicate that friction recovery can occur even in simple systems and can be related to the modulation of the surface over which the block slides. In the following, a case is presented where a similar phenomenology takes place in systems of granular particles.

III.2 Complex rheological behaviour of a granular system under external perturbation

Rheological properties are typically investigated by monitoring the dynamics of a driven probe immersed in the medium. A measure of the effective viscosity can be obtained by applying a constant force (or torque) to the probe and measuring its stationary angular velocity, as in the setup shown in Fig. 4. The dynamics of the probe is then studied as a function of the relevant parameters such as vibration frequency, amplitude or velocity, mechanical properties of the materials, or density and pressure. One of the main aims is to pinpoint the key quantities responsible for the different observed phenomena, such as friction reduction and jammed-unjammed transitions in the medium.

As discussed in the previous Sections, the relevance of the study of this kind of systems is twofold: On the one hand, one is interested in understanding how the rheological properties of the granular medium are modified in the presence of different kinds of perturbation; on the other hand, the response to external stimuli is relevant to clarify the microscopic mechanisms responsible for sliding frictional properties or for gouge failure mechanisms (see Sec. II). Concerning the first issue, the effect of flow rate and vibration amplitude on the rheological properties of glass beads has been recently investigated in several works, showing non-monotonic flow curves DWDDH11 ; gnoli2 and critical behaviour WDH16 . The effect of mechanical fluctuations on the probe has been studied in PADCC15 , where a generic rheological model is proposed. Other recent works addressed the interesting issue of non-local rheology (see e.g. BITCCA15 ), discussing different choices for the fluidity parameter, or the sound waves propagation (see e.g. santibanez2016 ).

Regarding the study of the effect of mechanical vibrations on frictional sliding, a common model system is represented by a particle layer confined between two substrates in relative motion persson , the granular medium acting as an intermediate lubricant layer in this case. The phenomenon of friction reduction in this kind of systems was discovered by Capozza et al. CVVZ09 ; CVVZ12 in numerical simulations of repulsive particles confined between a top horizontally driven plate and a bottom vertically vibrated substrate. In particular, focusing on the response of the driven plate upon varying the vertical vibration frequency, they observed suppression of friction in a well-defined range of frequencies. In the studied model, the vertical coordinate $Z$ of the flat bottom substrate follows the law

[TABLE]

while the top plate is driven through a spring of elastic constant $K$ moving at constant velocity $V$, similarly to the model described in Sec. II. The friction coefficient is defined as $\mu=F_{L}/F_{N}$, where $F_{L}=K[X(t)-Vt]$, with $X(t)$ the horizontal coordinate of the top plate and $F_{N}$ a normal external force acting on it. The average friction coefficient $\langle\mu\rangle$ shows a non-monotonic behaviour as a function of $f$, with a marked decrease in the interval $f\in[f_{1},f_{2}]$. Friction reduction in a well-defined range of frequencies shares strong similarities with the phenomenon observed in the spring block model C discussed above. The first transition frequency $f_{1}$ is related to the detachment condition and can be predicted imposing that the inertial force $F_{in}=M\ddot{Z}$, where $M=M_{p}+M_{top}$ with $M_{p}$ the mass of the particle layer and $M_{top}$ the mass of the top plate, overcomes the sum of the normal load $F_{N}$ and the damping force $F_{damp}=M_{p}\eta\dot{Z}$, where $\eta$ is the damping coefficient accounting for viscous dissipation. This condition gives CVVZ09

[TABLE]

The viscosity recovery frequency $f_{2}$ is obtained requiring that the detachment time from the bottom plate equals the period of the external oscillation and is related to the condition of maximum momentum transfer from the vibrating plate to the confined particles. This leads to

[TABLE]

These theoretical predictions are in very good agreement with numerical simulations and indicate that friction suppression is related to the reduction of effective interface contacts in the system due to the external vibrations. Finally, let us mention that recent studies LGV15 ; Lasta2 have highlighted the role played by the velocity of the imposed mechanical vibrations on the frictional properties of sheared granular media.

The friction reduction phenomenon occurring in a range of vibration frequencies in vertically shaken granular systems has been recently observed in experiments prl2018 . The setup is shown in Fig. 4: The granular particles (steel, glass or delrin spheres) are confined by an aluminium plate, with packing fraction $\sim 49\div 52\%$. The system is vertically vibrated by an electrodynamic shaker following Eq. (11), where $Z$ is the coordinate of the shaker plate. The probe is represented by a Plexiglas vane suspended in the medium and subject to an external torque. Further details on the experimental setup can be found in gnoli ; gnoli2 . For constant applied torque, the average angular velocity $\omega$ of the vane is proportional to the inverse of the macroscopic viscosity of the system and therefore provides information on its frictional properties. The measured values of $\omega$ are reported in Fig. 5 (left), as a function of $f$, for three values of the vibration amplitude $A$.

In agreement with the scenario described in Subsections III.1 and III.2, one observes two regimes, at low and high $f$, respectively. The first is characterized by high viscosity of the medium (corresponding to small values of $\omega$), whereas in the second fluidized regime, for intermediate values of $f$, $\omega$ reaches its maximum value $\omega_{max}$, corresponding to viscosity reduction, see Fig. 5 (left). The frequency threshold $f_{1}$ of the first transition from the solid (high frictional) state to the fluid (low frictional) state is very well estimated by the theory presented in Section III.2, leading to Eq. (12). In the experimental setup, the largest force provided by the shaker is $F_{in}=M\ddot{Z}_{max}$, with $\ddot{Z}_{max}=A(2\pi f)^{2}$, and $F_{N}=Mg$, with $M$ the total mass of the system (granular particles and aluminium plate). Thus, from Eq. (12), one has

[TABLE]

that is in very good agreement with the experimental results, see dashed lines in Fig. 5 (left). The generality of the underlying mechanism, related to the detachment condition, is demonstrated by experiments with different granular materials, reported in Fig. 5 (right). It is interesting to note that this fluidization phenomenon is different from the acoustic fluidization discussed in Section IV.

Increasing the vibration frequency one observes viscosity recovery for, say, $f\gtrsim f_{2}$, as shown in Fig. 5 (left). This phenomenon is analogous to the one observed in the numerical simulations described in Subsec. III.2. However, the main underlying physical mechanism is different. Indeed, the viscosity recovery frequency shows a strong dependence on the material properties, see Fig. 5 (right). The physical mechanism responsible for the phenomenon relies on a balance between dissipative and inertial forces and on the dissipation rate in the system, analogously to the case of spring-block model discussed in Section III.1. The dissipation rate is affected by the dissipative forces characterizing both the grain-grain and the grain-interface interactions. As shown in Fig. 5 (right), the recovery frequency is significantly reduced in the case where the bottom plate is covered with a thick layer of rubber tape (compare blue squares to green squares), namely $f_{2}$ decreases upon increasing the dissipation in the system.

The role of the dissipation rate has been investigated in more detail in numerical simulations prl2018 , considering a geometry where the granular medium is enclosed between two plates, the bottom one oscillating according to Eq. (11), and confined by the gravitational force. As a probe, a rigid cross-shaped subset of grains, subject to a constant force $F$ along the horizontal direction, is used. The velocity $v$ of the probe in the force direction for different $f$ and $A$ shows the fluidization transition at $f_{1}$ and the viscosity recovery at higher frequencies (see Fig. 6). The viscoelastic properties of the system are explored by changing the restitution coefficient of each grain $e$ silbert . In particular, two different restitution coefficients can be introduced to model the experimental situation: $e_{g}$ for grain-grain collisions, and $e_{b}$ for collisions between grain and bottom plate. The behaviour of $v/v_{max}$ as a function of $f/f_{1}$ for different values of $e_{g}$ and $e_{b}$ is shown in Fig. 6: $f_{1}$ is independent of $e_{b}$ and $e_{g}$, in agreement with Eq. (14), whereas the recovery frequency $f_{2}$ depends on the dissipation properties.

Insights in the physical mechanisms at the origin of viscosity recovery in this system can be obtained by the following reasoning. Viscosity recovery occurs when the rate of energy dissipation becomes larger than the rate of energy provided to the granular medium by the bottom plate oscillations. Indicating with $E_{in}$ the average energy input to the grains, and considering for simplicity a unique value of $e$, the rate of energy dissipation is given by $\dot{E}_{dis}\sim E_{in}(1-e^{2})f_{c}$, where $f_{c}$ is the collision frequency. It is reasonable to assume that $f_{c}$ is proportional to the vibration frequency $f$, yielding $\dot{E}_{dis}\sim(1-e^{2})fE_{in}$. On the other hand, at sufficiently large frequency, the rate of energy input $\dot{E}_{in}$ is controlled by the oscillation amplitude of the bottom plate and becomes independent of the vibration frequency $f$. As a consequence, also $E_{in}$ and $\dot{E}_{in}$ are $f$-independent and the condition $\dot{E}_{dis}\sim\dot{E}_{in}$ gives a recovery frequency $f_{2}\sim(1-e^{2})^{-1}$. This dependence of $f_{2}$ on $e$ is confirmed by numerical simulations in the case $e=e_{g}=e_{b}$, see inset of Fig. 6, showing the behaviour $f_{2}/f_{1}\sim(1-e^{2})^{-1}$. This analysis confirms that the second transition, leading to viscous friction recovery, relies on the dissipation mechanisms in the medium and between the medium and the bottom plate, explaining the observed dependence on the materials.

IV Acoustic Fluidization and Dynamical weakening

In this Section, the AF mechanism is considered as a possible explanation for the weakness of seismic faults. Indeed, many seismic faults exhibit a resistance to shear stress which is astonishingly lower than the one measured in laboratory experiments on both intact or ruptured rock specimens. AF represents an explanation for the observed weakness and can be considered an alternative to other interpretations, mostly based on the presence of water inside the fault zone. These interpretations, however, are in disagreement with measurements of porosity changes in sheared rocks Mel96 . On the other hand, the AF scenario does not invoke the presence of water or fluids but assumes the existence of an in-cohesive region inside the fault. If strong elastic (acoustic) waves perpendicular to the fault plane are activated, the normal stress is reduced and the in-cohesive region can slide under much smaller shear stress than required in the absence of vibrations. This process is based on the heterogeneous and noisy nature of the fracture propagation where the small wavelength component of the seismic radiation is scattered by small-scale heterogeneities inside the fault and eventually generates perpendicular oscillations. Melosh proposes a feedback mechanism: Vibrations promote the failure of a limited region of the fault which slips releasing its internally stored elastic energy. However, a fraction ${\cal E}$ of this energy becomes available to activate additional vibrations causing further failures. The self-sustainability of this process is based on the energetic balance between ${\cal E}$ and the fraction ${\cal E}_{diss}$ of energy scattered outside the fault or converted in heat. Using experimental values for ${\cal E}_{diss}$ and ${\cal E}$, Melosh has shown that acoustic standing waves of amplitude comparable to the confining pressure can exist in real seismic faults. Nevertheless, some initial acoustic energy must exists to trigger this cascading process. Melosh invokes stress-drop events in a limited area of the fault, generating a vibration over a sufficiently large volume, to allow acoustic energy to initially grow and self-sustain. The necessity of sufficiently large volumes is also proposed by Melosh as an explanation for the absence of AF in laboratory experiments but deeper insights on the initiation process are still missing. Furthermore, not only laboratory experiments but also instrumental measurements of radiated patterns of real earthquakes cannot provide a clear proof of the AF scenario. Indeed short wavelengths responsible for AF cannot propagate far from the fault zone and therefore can be observed only by instruments very close (few meters) to the rupture surface. Numerical simulations probably represent the only approach where the AF hypothesis can be concretely tested. The remaining part of this Section is devoted to the study of the AF scenario in numerical simulations of the fault model introduced in Sec. II.

IV.1 Laboratory observation of AF

A laboratory investigation of the AF hypothesis has been conducted in ref. XHM13 . In their experimental settings, acoustic excitation via compressive stress pulses, generated by high speed impact, are applied to a granular sample used to simulate fault gouge. Xia et al. have studied the stress-strain curve of the sample for different thickness of the granular sample as well as for different choices of the striker impact velocity $V_{im}$. For a fixed thickness the experiment shows a dramatic change in the rheological properties at a characteristic impact velocity $V_{c}$ such that the granular sample exhibits a solid-like behaviour for $V_{im}>V_{c}$ and a fluid-like behaviour for $V_{im}<V_{c}$. This result is consistent with the AF scenario since larger values of $V_{im}$ cause an increase of the confinement stress. As a consequence, for $V_{im}>V_{c}$ the internal pressure produced by acoustic oscillation of grains is no longer able to counterbalance the confining pressure and a solid-like behaviour is observed. This interpretation is supported by the decrease of $V_{c}$ as function of the granular thickness, consistently with the AF scenario, as well as its decrease when pre-stress granular samples are considered. An indirect evidence of AF in laboratory investigation is the clear reduction of the granular layer thickness produced by acoustic vibration and documented in ref.VEBBJ12 . This reduction, indeed, has been attributed to vibrations generated inside the sample which lead to the auto-acoustic compaction according to a mechanism similar to what predicted by the AF hypothesis.

IV.2 The AF frequency in the numerical model of seismic faults

In this Section it is explored the possibility that vibrations perpendicular to the fault plane can activate inside the fault and are able to self-sustain. The most natural hypothesis is that these vibrations correspond to standing waves that travel vertically between the two confining plates. If these waves propagate with velocity $v_{a}$ their frequency is given by $\omega_{AF}=\pi v_{a}/W$, where $W$ is the width of the layer. In order to provide an expression for $\omega_{AF}$ in terms of the model parameters, we consider that the propagation velocity is $v_{a}=\sqrt{M/\rho}$, where $M$ is the $P$-wave modulus and $\rho$ the system density. The evaluation of $M$ in confined granular media is very complicated and, indeed, experimental and numerical studies Makse04 ; ROLJJ15 indicate that it increases for increasing confining pressure. However, since the shear modulus of a granular packing is negligible, compared to the bulk modulus, in first approximation $M$ coincides with the bulk modulus of a single grain. In the mechanical model each grain under normal compression can be considered to deform as a cube. More precisely, a compressional stress $\sigma_{ii}$ applied in the $i$-th direction, on the two faces perpendicular to the $i$-th direction, produces a deformation $2\delta x_{i}$, along the $i$-th direction, with $\delta x_{i}$ given by $k_{n}\delta x_{i}=\sigma_{ii}d^{2}.$ For $\delta x_{i}\ll d$, one has $\frac{\delta V}{V}=-2\sum_{i=1}^{3}\frac{\delta x_{i}}{d}=-\frac{2d}{k_{n}}\sum_{i=1}^{3}\sigma_{ii}=-\frac{6d}{k_{n}}P$, where $P$ is the applied pressure and $V$ the volume. The $P$-wave modulus then is equal to $M=k_{n}/(6d)$, leading to the AF frequency

[TABLE]

IV.3 Identification of AF in the numerical model

IV.3.1 The response to an external perturbation

In this Section the response of the system to an external perturbation is discussed as function of the perturbation frequency GSdALP15 . A granular fault of width $W=10d$, confined by $P_{0}$, with normal spring constants $k_{n}=210^{3}P_{0}d$ is subject to an external perturbation at frequency $\omega^{*}=1.4\pi$ which, according to Eq.(15) corresponds to $\omega_{AF}$in our model. More precisely, at a certain time $t$ the external drive is stopped, i.e. $V$ is set to zero, and for a temporal period $\tau$ the confining pressure is changed of a quantity $P_{-}(t_{p},t)=-\frac{\alpha P_{0}}{2}\left[1-\sin(\frac{\pi}{2}+\omega(t_{p}-t))\right]$, with $t_{p}\in[t,t+\tau]$ and $\omega=\omega^{*}$. The duration of each perturbation is fixed to $\tau=10$, which for varying frequencies leads to $n=\tau\omega/2\pi\in[1,10^{3}]$ pulses. The reduction of the confining pressure is expected to promote a displacement of the top plate $\Delta x(t)=x_{\alpha}(t+\tau)-x_{0}(t+\tau)$, where $x_{\alpha}(t+\tau)$ is the position of the top plate after the perturbation has been applied, and $x_{0}(t+\tau)\simeq x_{0}(t)$ is the unperturbed top plate position. The response to the external perturbation is then quantified by the function $\Pi_{-}(t,\omega^{*})=\Delta x(t)/\alpha P_{0}$ which, for sufficiently small values of $\alpha$ ($\alpha<0.05$), becomes $\alpha$-independent. The evolution of $\Pi_{-}(t)$ is monitored in the temporal window $[t_{s_{0}}:t_{s_{3}}]$ (Fig. 7) which presents three slips of size $\Delta x_{\rm slip}(t_{s_{i}})<0.1$, at times $t_{s_{0}}=0,t_{s_{1}}$ and $t_{s_{2}}$, followed by a large slip $\Delta x_{\rm slip}(t_{s_{3}})=20.3$, at time $t_{s_{3}}$.

The quantity $\Pi_{-}(t,\omega^{*})$ (left panel Fig. 8) increases as $t$ approaches the slip occurrence time. More precisely the external perturbation induces a displacement of the top-plate $\Delta x(t)$ which is comparable to the slip in the unperturbed evolution if $t$ is sufficiently close to the slip occurrence time $t_{s_{i}}$. As a consequence, the main effect of the external perturbation is to anticipate the slip occurrence and one can define a time advance $\Delta t_{a}$, by the condition $\Delta x(t_{s}-\Delta t_{a})=0.2\Delta x_{\rm slip}(t_{s})$. This result shows that external perturbations can reduce the confining pressure and can weaken the fault in such a way that slip occurs for a smaller value of the applied shear stress similar to what observed in the numerical study FGGJMC13 .

In the following evidence is provided that the observed weakening is promoted by acoustic oscillations at the frequency $\omega_{AF}$, as in the AF hypothesis, by showing that

• •

(a) A behaviour qualitatively similar to $\Pi_{-}(t,\omega^{*})$ is observed if the external perturbation increases the confining pressure or reduces the shear stress;

• •

(b) The response $\Pi_{-}(t,\omega)$ strongly depends on the perturbation frequency $\omega$ and $\Pi_{-}(t,\omega)\simeq 0$ if $|\omega-\omega_{AF}|\gg 0$;

• •

(c) The response only weakly depends on the duration $\tau$ of the applied perturbation.

To address item (a) the same analysis is performed under an external perturbation which increases the confining pressure

[TABLE]

with $t_{p}\in[t,t+\tau]$ or with a perturbation which forces the shear stress to vary by $\sigma_{\pm}(t_{p},t)=\pm\frac{\alpha\sigma(t)}{2}\left[1-\sin(\frac{\pi}{2}+\omega(t_{p}-t))\right]$. As in the definition of $\Pi_{-}$, the response of the system in the linear regime is quantified by $\Pi_{\pm}(t)=\Delta x(t)/\alpha P_{0}$ and $\Sigma_{\pm}(t)=\Delta x(t)/\alpha\sigma(t)$. Results (right panel Fig. (8)) show that the system response to the external perturbation does not depend on the orientation of the perturbation and even perturbations which should stabilize the system as $P_{+}$ or $\sigma_{-}$ are able to anticipate the slip occurrence.

Concerning item (b) in Fig. (8) $\Pi_{\pm}(t,\omega)$ and $\Sigma_{\pm}(t,\omega)$ are plotted for different values of $\omega$. Results clearly show that the larger the difference between $\omega$ and $\omega^{*}$ the weaker is the response. This effect can be quantified by the dependence of $\omega$ on the advance time $\Delta t_{a}$ (left panel Fig.9). For all kinds of perturbation $\Delta t_{a}$ presents a peak at $\omega=\omega^{*}$. In particular for $\Pi_{+}$ and $\Sigma_{-}$, $\Delta t_{a}$ approaches a Dirac-delta function indicating that the system is weakened only by perturbation at the frequency $\omega^{*}=\omega_{AF}$. Finally, to support item $(c)$ $\Pi_{\pm}(t,\omega)$ and $\Sigma_{\pm}(t,\omega)$ are evaluated under perturbations of different durations $\tau$. Results plotted in Fig.(9) confirm that $\Delta t_{a}$ only weakly depends on $\tau$.

In order to verify that fault weakness is promoted by the acoustic oscillations described in Subsec.IV.2, the same analysis is performed considering systems of different width $W$ and with different grain elastic properties (different $k_{n}$) leading to different values of $\omega_{AF}$ according to Eq.(15). The frequency $\omega_{max}$ which produces the maximum value of the advance time $\Delta t_{a}$ (right panel Fig.9) is $\omega_{max}\simeq\omega_{AF}$ for all values of $k_{n}$ and $W$ clearly supporting the AF scenario. Finally, it is important to emphasize that item (a) indicates that AF can be triggered by perturbations applied along any direction. This provides a possible explanation of triggering caused by transient seismic waves regardless the fault orientation.

IV.3.2 Endogenous activation of acoustic oscillations

In this Section, following ref.GSdALP15 , it is explored whether oscillations at the frequency $\omega_{AF}$ can be self-activated promoting slip instabilities at shear stress levels smaller than the one required in the absence of oscillations.

For this study, at each time $t$, a replica of the system is created under constant external drive ($V=0$) and its spontaneous relaxation is followed in the subsequent time interval. The autocorrelation function of the particle velocities $\bf{v}_{i}$ is evaluated as $C(t,t^{\prime})=\frac{\sum_{i=1}^{N}\bf{v}_{i}(t)\cdot\bf{v}_{i}(t^{\prime})}{\sum_{i=1}^{N}\bf{v}_{i}(t)\cdot\bf{v}_{i}(t)}$. Given a slip with occurrence time $t_{s}$, in Fig.10 $C(t,t^{\prime})$ is plotted as function of $t^{\prime}$ focusing on nine different values of $t$: Three times (blue circles, panels (a-c)) are located much before the slip time $t<t_{s}$, three times (red circles, panels (d-f)) immediately before the slip instability $t\lesssim t_{s}$ and three times (green circles, panels (g-i)) after the slip time $t>t_{s}$. The panels (d-f), corresponding to $t\lesssim t_{s}$, show that $C(t,t^{\prime})$ presents a clear oscillating pattern with a typical frequency $\omega\simeq\omega_{AF}$. Oscillations at the same frequency are also observed after the slip (panels (g-i)), but with an amplitude decreasing with increasing $t$. Conversely, before the slip (panels (a-c)) $C(t,t^{\prime})$ exhibits an irregular pattern. The presence/absence of oscillations at $\omega=\omega_{AF}$ can be detected by the Fourier transform $\hat{C}(t,\omega_{AF})$ of $C(t,t^{\prime})$ with respect to the variable $t^{\prime}$. Results (upper left panel) show that $\hat{C}(t,\omega_{AF})$ is very small in the short temporal period anticipating the slip. However, as soon as the slip time is approaching $t_{s}$, $\hat{C}(t,\omega_{AF})$ rapidly increases, reaches its maximum value in correspondence to the slip occurrence time and then remains substantially constant. The evolution of $\hat{C}(t,\omega_{AF})$ on a longer temporal period is plotted in the lower right panel of Fig. 10. This figure confirms that $\hat{C}(t,\omega_{AF})$ presents a very fast increase as the slip instability is approaching. After the slip, conversely, it decays reaching small values at times $t$ distant from slip instabilities. In the same figure $\hat{C}(t,\omega)$ is plotted for other values of $\omega$: If $\omega\neq\omega_{AF}$, $\hat{C}(t,\omega)$ is a noisy quantity fluctuating around a very small value. This figure clearly shows that characteristic oscillations appear only at a frequency $\omega\simeq\omega_{AF}$ at the onset of the slip instability. These oscillations then tend to disappear as the dynamics goes on. This pattern is confirmed by the contour map plot in the upper right panel of Fig. 10.

IV.4 Mechanisms producing AF in the numerical model of seismic faults

In order to identify the mechanisms responsible for the activation of these oscillations, the trajectories of each grain inside the stick phase can be monitored in numerical simulationsGdAPL18 . Because of the high granular density, the large majority of particles is always in contact with their neighbours forming an almost rigid structure, i.e. the backbone. Conversely, a small fraction (less than $10\%$) of particles, the rattlers, are located inside the cages formed by the particles in the backbone TS10 and most of the time do not interact with other particles, even if they contribute to stabilize the backbone structure GdAPL17 . In Fig. 11 (first panel) the $x$-position of four neighbouring particles is plotted during a short time window far from slip: Three particles $(x_{1},x_{2},x_{3})$ exhibit regular oscillations along the $x$-direction. Differently, particle $4$ is a rattler and moves along a straight line up to an abrupt change in direction caused by a collision. As shown in the left panel of Fig.11, backbone particles exhibit an oscillating behaviour along the $x$-direction with a characteristic frequency $\omega\simeq\omega_{AF}$. A similar oscillating behaviour is also observed for the $y$ and $z$ component of the particle velocity, with oscillations along the $z$-axis much smaller in amplitude. Superimposing the centers of each trajectory in a common point, as in Fig.11 (second panel), trajectories appear to be roughly confined in a plane and exhibit an elliptic-like shape.

These elliptic trajectories cannot be stable during the whole stick-slip dynamics. Fig.10 indeed has shown that the correlation function $C(t,t^{\prime})$ presents irregular behaviour with $\hat{C}(t,\omega)\simeq 0$ in a temporal window anticipating the slip instability $t_{s}$. This indicates that grain velocities decorrelate during the system evolution. The origin of the decorrelation is in the change of orientation of the elliptic trajectories. In order to prove this point, the angle $\theta_{i}$ formed by the $i$-th particle velocities with the $z$ axis is evaluated, together with the distribution $P(\theta)$, i.e. the number of particles whose velocity has orientation $\theta_{i}\in[\theta,\theta+\Delta\theta)$. Fig. 12 (left panel) shows $P(\theta)$ evaluated for different slips and at different times before and after $t_{s}$. For all slips, at large temporal distances from $t_{s}$, $P(\theta)$ is sharply peaked at $\theta\simeq 90^{o}$ (open symbols), corresponding to an oscillatory motion in the $x-y$ plane. This distribution does not change significantly during the evolution and only in proximity of the slip time it spreads towards smaller values of $\theta$ (filled symbols). Therefore, as $t$ approaches $t_{s}$ oscillations become present also in the direction parallel to the $z$-axis ($\theta\simeq 0$). This is confirmed by the behaviour of the $z$-coordinate as function of time and as function of the $x$-coordinate (Fig.11). Far from the slip (first panel of Fig.11), the displacement in the $z$-direction presents oscillation at the frequency $\omega_{AF}$. As already observed $z-$displacements are small compared to the $x-$displacements and the trajectory is mostly confined in the $x$-$y$ plane ($\theta=90^{o}$) (second panel of Fig.11). At the onset of slip instability (third and fourth panel of Fig.11) the angle $\theta$ is no longer stable and $z-$ displacements of size comparable to $x-$displacements are indeed observed. The above findings support the hypothesis of weakening by AF. Indeed when oscillations are confined in the $x$-$y$ plane ($\theta\simeq 90^{o}$) they do not affect the confining pressure. Conversely, when $\theta\simeq 0^{o}$, oscillations can reduce the confining pressure promoting failure.

The overall picture is confirmed by the behaviour of $P(\theta)$ (right panel of Fig.12) when the perturbation $P_{+}(t,\omega)$ (Eq.16) is applied at the resonant frequency $\omega_{AF}$, with $\alpha=0.05$. Differently from the unperturbed evolution when $P(\theta)$ dramatically changes as the slip instability is approaching (left panel of Fig.12), in the presence of the perturbation $P(\theta)$ is substantially time independent (right panel of Fig.12). In particular the presence of vertical oscillation $\theta\simeq 0$ are observed at all times. Combining this observation with the behaviour of $\Pi_{+}$ (Subsec. IV.3.1), it is possible to conclude that external perturbations at the frequency $\omega\simeq\omega_{AF}$ induce the alignment of the orientation of ellipses along the vertical direction reducing the confining pressure and promoting the time-advance in the slip occurrence. Interestingly, close to the slip instability, only small differences are found in the angle distribution $P(\theta)$ with or without the perturbation. This indicates similarities in the mechanisms leading to the endogenous activation of vertical oscillations in the unperturbed evolution and the activation triggered by the applied perturbation.

V Conclusions

Occasionally, a major earthquake can trigger additional quakes up to 1000 km from the epicentre of the first event. How that happens is not clear because the strength of the seismic waves decreases the farther they travel. A numerical approach proposes an explanation for this observation in terms of seismic waves from the original earthquake that can create a lathering effect in the grains of the gouge between tectonic plates. The fluidization of the granular medium reduces the friction between the plates, which increases the likelihood of slipping. The model reveals that the frequency of the seismic waves is the only important variable. Even if the amplitude of the seismic waves is very small, regardless their propagation direction, they are able to trigger an earthquake if their frequency matches the resonance frequency of the fault gouge $\omega_{AF}$. Results of Sec. IV indicate that vibrational modes at the characteristic frequency $\omega_{AF}$ do not form at the onset of slips but are already present inside the system at all times. Reasonably, the energy responsible for these oscillations originates from the energy stored, during the stick phase, through the spring which couples the system to the external drive. Most of this energy is released very rapidly during the slip but a significant fraction contributes to the activation of harmonic oscillations. Because of the vertical confinement, only the mode at the frequency $\omega_{AF}$ (Eq.(15)) can be a standing wave inside the fault. Even if these modes have been explained in terms of compressional waves propagating along the $z-$direction, because of the heterogeneous structure of the granular packing, these waves induce also displacements along the $x-$ and $y-$ directions. Far from the slip, the confinement along the $z-$direction and periodic boundary conditions along $x$ and $y$, lead to $x-$ and $y-$ displacements larger than $z$-displacements. Conversely, as the system evolves there exists a finite probability that ellipses rotate activating oscillations along the $z-$direction. A possible mechanism responsible for the ellipse rotation can be identified in the collisions of rattlers with backbone particles. Indeed, these collisions can be sufficiently energetic to destabilize oscillations originally confined in the $x-y$ plane. Within this scenario, it is impossible to forecast in advance the occurrence time $t_{s}$ of the next slip. Indeed, rattlers follow chaotic trajectories and the occurrence time of sufficiently energetic collisions appear to be non-predictable.

Beside the relevance for seismic occurrence, the fluidization of a granular medium under periodic perturbations is of extreme relevance in a variety of fields, for many phenomena, from avalanche dynamics to the manufacturing process in material, food, and pharmaceutical industries. Both experimentally and numerically, different frictional regimes have been observed in the system, from very large viscosity at low vibration frequencies to fluidized states (corresponding to viscosity reduction) at intermediate $f$, with a viscosity recovery at higher values of $f$. The first transition to the fluidized state is well characterized by the detachment condition of the granular medium from the confining plate and is independent of the material properties. The second transition, leading to viscosity recovery, is related to dissipation mechanisms in the medium and between the medium and the bottom plate and therefore shows a strong dependence on materials. These observations suggest the possibility to control the viscous properties of confined granular media by tuning the shaking frequency in the system, with important practical application in several fields, from tribology to geophysics and the material industry.