Earthquake dynamic induced by the magma up flow with fractional power law and fractional-order friction

In this work, we examine the dynamical behaviour of the “single mass-springs” model for earthquake subjected to the strength due to the up flow of magma for the period of volcanism, considering the fractional viscous damping force, the fractional weakening friction and fractional power law of elastic force. The numerical simulation method used in this paper is that of Grünwald-Letnikov based on the generalization of the classical derivative, and the approximately analytical solution obtained by the harmonic balance method. The results have shown that the fractional-order derivative can affect the dynamical properties of fault rock, which is characterized by the equivalent damping coefficient and the equivalent stiffness coefficient. Moreover, the amplitude-frequency equation for the steady-state solution was established. It appears that the resonant amplitude and resonant frequency are strongly dependent on the fractional-order damping r, fractional-order friction q, the fractional deflection 𝛼 , the nonlinear stiffness coefficient and the fractional viscous coefficient. We have also shown that, the recurrence time of an event, the duration time of an event and the slip size of an earthquake can be controlled by the fractional-order derivative, the fractional-order deflection and the magnitude of the magma strength. The model allowed us to better interpret the earthquake as a stick-slip motion.


Introduction
An earthquake is a ground shock resulting from the sudden release of energy accumulated by the stresses on the rocks [Kanamori, 2001]. The breaking along a fault provokes this release of energy. Most earthquakes come from pre-existing faults or new breaks in the rocks that make up the earth's crust, along which rocks overlap each other.
As the rock is stretched, potential energy builds up. Eventually, it becomes so large that it is suddenly released in the form of earthquakes. The seismic waves move away from the source of the earthquake, causing shaking on the surface of the earth, also called ground acceleration. During a seism, the potential energy stored in the earth's crust is released and can produce heat, cracks, and deformations. The seismic energy is radiated in the form of seismic waves. This radiated energy is defined by [Choy and Boatwright, 1995]. In the mechanics of faulting, we still ignore the specific laws which control the released energy during an earthquake . Indeed the energy stored in the earth's crust can be controlled by dynamic friction [Bhattacharya and Rubin, 2014], pore fluid pressure [Scholz, 1990], thermal pressurization [Bizzarri, 2011] and brittle-ductile fracture rheology [Wang, 2013]. To understand more the earthquake dynamics and other phenomena that appeared during an earthquake, some authors used the fractional calculus method [El-Misiery, 2006;Pelap et al., 2018;Tanekou et al., 2018].
Fractional calculus has been extensively studied in recent years and many results have been given: Baleanu et al. in [2010] Caputo and Fabrizio in [2015], Atangana et al. in [2016], Podlubny in [2009], Torres and Malinowska in [2012], Kilbas et al. in [2006], Petras in [2011], Tanekou et al. in [2018], Pelap et al. in [2018], and others. In the last decade, the fractional calculus mainly received attention due to its essential role in modelling the anomalous dynamics of various processes related to complex systems in several fields of science and engineering. Fractional differential equations have found many applications in physics, control engineering, and sig[al processing [Miller andRoss, 1993, Oldham andSpanier, 1974;Laskin, 2002;Baleanu and Golmankhaneh, 2009;Chung, 2015;Eslami, 2016]. Fractional derivation has a particular advantage for the mechanical modelling of materials whose behaviour is called viscoelastic, such as rubbers, gums, and rocks [Bagley and Trovik, 1984;Koeller, 1985;Bagley and Calico, 1991;Bohdana, 2007;Ali, 2012]. The materials that preserve the memory of a past deformation and whose behaviour is said to be viscoelastic, the fractional derivatives are naturally introduced them [Mohammad and Yousef, 2017].
The seismic fault being made up of the viscoelastic rocks, its modelling will not be done without introducing the fractional derivative. Also, the earthquake occurs throughout the fault and not at a given point of the fault, the whole-order derivative having a local character will not be appropriate for earthquake modelling, hence the use of the fractional derivative. Since earthquake is only the result of past deformations (load phase of the fault), for its excellent, we will use an instrument that will take into account past events [Hong et al., 2018]. Given all this, the fractional derivation is introduced naturally because it has a memory effect. Several scientists had already studied this memory effect, so from 1979 to 1983, Dieterich and Ruina introduced a state variable θ in their model Dieterich [1979; and Ruina [1983], this variable θ was interpreted as a function of the history of sliding [Pomeau and Le Berre, 2011]. Well after this memory effect has been studied by Kostic et al. [2013] through the temporal delay τ introduced in the friction term.
Recent experimental observations have shown that inhomogeneous micro media manifest nonlinear elastic behaviour [Guyer and Johnson, 2000] which cannot be explained by the classical theory of nonlinear elasticity proposed by Landau and applicable to homogeneous and crystalline media. The nonlinear mesoscopic elasticity was introduced to describe the properties of such materials, which include most types of rocks, concrete, bones, some damaged metals and composites, etc. [Adams, 1998;Guyer and Johnson, 1999].  [Guyer and Johnson, 2000;Scalerandi et al., 2008]. Recently an energy harvesting system with fractional-order viscoelastic material has been studied by Cao et al. [2014] and Kwuimy et al. [2015]. They showed that the fractional order property of the material enhances high-energy chaotic motion as well as inter-well periodic oscillation.
Active fault near the potentially active volcanoes can be destabilized during the volcanic eruption. Complex interactions between fluids and rocks may generate long-period ground vibrations that precede and accompany eruptions. Volcanic eruptions require a flow of magma and/or aqueous fluids through rock, and there is potential for long-period seismic signals to provide valuable information on changes in the location, velocity, and types of fluids (e.g., gas, magma, bubbly magma) under volcanoes. However, such analysis requires understanding potential source mechanisms of the ground oscillations and the characteristics of the resulting signals [Julian, 1994]. There are numerous potential sources of the pressure disturbance required to trigger resonance including an earthquake, Christian Foka Fogang et al.
a new crack network connection, shock waves from "choked" flow [Morrissey and Chouet, 2001], or bubble coalescence leading to a rising gas slug [Cruz and Chouet, 1997]. These are all credible sources of LP-events, which decay after seconds or tens of seconds. However, tremor requires a disturbance that is sustained for minutes or even months and thus precludes transient resonance triggers such as earthquakes. Sustained resonance could be caused by the continued formation of shock waves or gas slugs however these seem to require exceptional circumstances and do not explain the near-ubiquity of tremor during volcanic eruptions of all styles and compositions [Julian, 2000]. If fluid flows through a crack induced oscillations in the conduit walls, this could be a source of sustained seismicity, lasting as long as the flow continued at a sufficient speed.
Julian [1994] explored this tremor mechanism using an approach similar to [Pedley, 1980] for blood flow through collapsible arteries. Julian set up a lumped-parameter model involving the two-dimensional flow of an incompressible Newtonian fluid through a slot-like constriction with two movable walls (Figure 1). At each end of the constriction, the conduit is wide in such a way that it can be considered to be a fluid reservoir with constant fluid pressure despite oscillations of the constriction. Each wall is modeled as a mass, whose motion is controlled by a spring representing the elasticity of the country rock, and dashpot for inelastic effects and radiation damping.
Separation of the walls changes as a function of time only, and thus the walls of the constriction are always parallel.
Viscous damping plays a significant role in the frictional strength and stability of tectonic faults. Some authors [Bizzarri, 2012a;Dragoni and Santini, 2015] consider that viscosity plays a role in causing seismic radiation to release strain energy during faulting. Bizzarri and Cocco [2006] revealed that the melting of rocks and fault gouge is likely to occur even with the inclusion of the thermal pressurization of pore fluids. When fluids are present in faults, thermal pressurization can play a significant role in earthquake rupture and also results in resistance on the fault plane [Wang, 2013]. When melting occurs, the rheological behavior of the fault zone no longer obeys the Coulomb-Amonton-Mohr formulation [Bizzarri, 2011] and during the co-seismic phase, all of the strain accumulated during the inter-seismic duration does not recover [Namiki et al., 2014]; a fraction of this strain remains as a result of viscous relaxation, in that a viscous rheology is needed to describe the rupture processes of an earthquake. The frictional strength increases with the depth as a function of the pressure, whereas the viscous strength decreases with the depth as a function of the temperature; consequently, rocks deform ductility beyond the depth known as the brittle-ductile transition [Burgmann and Dresen, 2008]. The rheology of the lower layer may not be a simple Newtonian viscosity [Freed et al., 2006] and is commonly modeled using Maxwell or Burgers bodies [Haberland, 1981], which respond to a sufficiently slow deformation as a fluid-like material, i.e., permanent displacement can remain without accumulation of elastic energy and shallow mantle viscosities can become as low. The viscosity description is mentioned in more detail by Wang [2017a, b].
Rocks of earth's crust present an anomalous elastic behaviour, and earthquakes are complex systems with longmemory, and some of their faults have fractal properties. It is fair to claim that fractional-order derivative and fractional power-law are very efficient for earthquake modelling. In this paper, we assume that the friction forces have fractional-order derivatives, while the elastic properties of the rock fault are modelled by using fractional order power law as discussed by Cvetivanin [2009], Lewis and Monasa [1981] and Lee [2002]. 4 The present paper is organized as follows. The next section presents the preliminaries on the fractional power law. Section 3 presents the physical system and the corresponding mathematical model. In Section 4, the resonance of the fractional-order system is analytically studied by the harmonic balance method. The numerical results are presented in Section 5. The last section is the conclusion.

Power law in the non-linear elastic material
Nonlinear elastic deformation of damaged rocks and other brittle materials is well-established by observations at various scales [Johnson et al., 1996;Lockner and Stanchits, 2002;Basaran and Nie, 2004]. Distributed rock damage in the form of cracks, joints, and other internal flaws develops as part of the rock formation and can increase significantly during tectonic loading. Several experimental studies under conditions of brittle deformation revealed that the elastic parameters strongly depend on the extent of the damage, leading to vanishing effective elastic modulus at significant stresses just before failure [Johnson et al., 1996;Hamiel et al., 2004]. For flexible material, Debnath [2003] reported that the stress-strain law of elasticity theory could be generalized using fractional order power law. In their study of the deflection of a cantilever beam made of a Ludwick type material under one concentrated vertical load at the free end, Lewis and Monasa [1981] experimentally obtained for the nonlinear part of a stress-strain relationship = 66.500 ⁰ . ⁴⁶³ ( being the stress and the strain] by measuring the deflection at the free end. In Refs. [Lewis & Monasa, 1981;Lee, 2002;Shatarat et al., 2009;Bank et al., 2010], various fractionalorder power laws were experimentally obtained for the stress-strain relation of different flexible materials.
The elastic force with a fractional-order power law in the flexible material is given by [Cvetivanin, 2009]: (1) where 1, 3 and are, respectively, the linear and nonlinear stiffness coefficient and the fractional-order. There are many materials for which the fractional-order is larger than 1 and not an integer. The negative nonlinear stiffness coefficient ( 3 ≻ 0) models a hardening nonlinearity, and ( 3 ≺ 0) models a softening nonlinearity [Miller, 2005].

Mathematical model
The model (see Figure 2) considered here consists of a slider with a mass m connected to a rigid pulling rod moving at a small constant velocity V by a linear spring with a constant and the nonlinear viscous damping . The block is connected to the wall conduit by a spring with a constant subjected to the magma thrust strength defined by [Pelap et al., 2016]: (2) Figure 2. Springs-block model of an earthquake subjected to the magma thrust strength [Pelap et al., 2016].
The friction force between the block and the surface is a function of the instantaneous velocity of the block concerning a characteristic velocity . The equation of motion with nonlinear viscous damping via the fractionalpower law of the velocity is given by the following expression: ( 3) where denotes the position of the block measured with respect to the position where the sum of the elastic forces on the block is zero, sgn(.) is the sign function. and are respectively the nonlinear stiffness coefficient, and the fractional-order, and 1 are respectively the viscous damping coefficient and the nonlinear damping coefficient, and are the fractional-order derivatives and denotes the displacement of the wall conduit with mass , defined as: (4) ( / / ) is the fractional-order weakening friction force define as: (5) where represents the friction force amplitude along the dry fault ( is the force that enables the block to resist the movement induced by tectonic plates under dry conditions) and ( ) designates a continuous function for ≥0 that vanishes for large values of and satisfies two conditions, namely, (0)=1 and '(0)=-1, wherein the prime denotes differentiation with respect to the argument and the second condition expresses the velocity-weakening effect of the friction [De Sousa Vieira, 1992]. By integrating Equation (4), we obtained: The fractional-order derivative is introduced to the damping force in Equation (3) due to the fact that many materials or devices exhibit obvious visco-elasticity, which could be precisely modelled by fractional-order derivative.
According to Stokes' law, the damping coefficient of a sphere of radius R in a fluid of viscosity is given by =6 [Kittel et al., 1968]. The viscosity, , of a fault zone depends upon the temperature, pressure, water content, SiO 2 content, rock types, etc. Spray [1993] showed that the calculated values of for friction melts are low and decrease with increasing temperature. In general, lies in the range 10 15 Pa.s and 10 20 Pa.s when the temperature varies from 500° C to 1200° C. Savage and Lachenbrush [2003] mentioned that for a seismogenic zone with a thickness of 14 km, the viscosity on the bottom is (0.31 -3.12)*10 20 Pa.s. Obviously, the values for lower crustal materials from Savage and Lachenbrush [2003] are higher than those for upper crustal materials from Spray [1993].
The parameter R is almost the dimension of the deformed volume around a rupturing fault. Based on stick-slip behaviour and elastic rebound theory of a fault, Turcotte and Schubert [1982]  N/m [Wang, 2012;Wang, 2017a]. Hence, the value of lies in the range 10 -8 to 10 -6 for of 10 2 to 10 15 Pa.s. The parameters adopted in this study are regrouped in the Table 1.

Parameters Units
Fault stiffness, and /
Before proceeding with the analysis, let us introduce the dimensionless variables [De Sousa Vieira, 1992] = / ₀ and = so that the equation of motion Equation (3) can be written in the following dimensionless form: where and = ₁ ₀⁻¹/( ) / ².
The value of must be much smaller than 1 because of ≈10⁻¹⁰ . ⁻¹. Since the value of mainly influences the recurrence time between two events, and can only makes a very small influence on the pattern of time variations in velocity and displacement of events. In order to study long-term earthquake recurrence, there must be numerous modelled events with clear and visualized time functions of displacement and velocity for an event in the computational time period. If = 10⁻¹⁰ is considered, the recurrence time of events is very long and thus the duration time of an event is much shorter than recurrence time. This makes the time function of an event displayed in the long-term temporal variation in slip looks like just a step function for the displacement and an impulse for the velocity [Wang, 2018]. Hence, in order to get fine visualization a larger value of is necessary. The value of is usually very small during an event and cannot influence the rupture because of very tiny value of . Numerical test shows that when > 10⁻², the value of is not small during an event and can influence the rupture [Wang, 2018].
Increased interest in stick-slip instability was renewed during the late 1970s through laboratory rock experiments to understand earthquake rupture processes. Dieterich [1978], Perrin et al. [1995], Rice [1983], Cocco [2005, 2003], and others used these experiments to formulate constitutive laws capable of describing the frictional stresses when rocks are sheared against each other or over a surface.
The fault behaviour strongly depends on the fault constitutive model used. Therefore, in the following text, we use the velocity-weakening frictional force given by Langer, 1989a, 1989b]:

Amplitude-frequency of the system
In this section, we used the method of harmonic balance to investigate the frequency response of the Equation (7).
Several analytical techniques are available for the analysis of dynamics of systems with fractional derivative and fractional power law, such as the residual harmonic balance method, multiple time scales method and the homotopy method [Leung et al., 2012;Xiao et al., 2013]. In this part, we will focus on the weak displacements of the block. (i.e) ( / ≺≺ ), the function will be developed in the linear approximation in terms of the Taylor series according to Langer and Tang [Langer and Tang, 1991]  Fractional earthquakes under magma flow (9) Owing to this model of friction, the dimensionless equation, Equation (7) becomes (10) The solution of the governing equation, Equation (10) is found in the form: and are amplitude and phase respectively. Considering, Equation (11) for each term of Equation (7), one obtains: Definition [Podlubny, 1999] of the fractional derivatives allows us to obtain the expression of (12) where (13) Although the integrals Equation (14) The approximation leads to (16) with Ψ( )=Ω -( ) Where (20) The nonlinear stiffness term is evaluated neglecting higher-order term as follows: The right-hand side of Equation (21) can be decomposed in a Fourier series according to Where and are defined as: Where ( , ) denotes the beta function.
One obtains finally (26) The nonlinear viscous term is given as follows: The right-hand side of Equation (36) can be rewritten as: With Θ = Ψ+ /2. Equation (28) has the same form as Equation (21). Therefore Equation (28) can be written in the form: wher.
Substituting thes e previous transformations into Equation (10), we obtain the following first-order differential equations: Christian Foka Fogang et al.

Fractional earthquakes under magma flow
After one development and simplification, Equation (31) becomes: Where: (34a) We define the parameters and as the equivalent linear damping coefficient and the equivalent linear stiffness coefficient respectively.
In the integer model ( = 1, = 1), the equivalent linear damping coefficient and the equivalent linear stiffness coefficient are given by Equation (35).
Based on Equation (34), it could be found that the fractional-order derivatives, and and the fractional damping coefficient are all significant in system dynamics. The fractional damping coefficient is linear with an equivalent linear damping and equivalent stiffness in this system. Moreover, it could be found that the equivalent linear damping and stiffness coefficients are also associated with the excitation frequency. The equivalent linear damping and linear stiffness coefficient versus the fractional-order derivatives and are shown in Figure 3. It appears in this figure that, the equivalent linear damping coefficient is little and the equivalent stiffness coefficient is more significant when the fractional-order derivative → 1. These results are remarkably different from the general consideration in most literature, where the fractional-order derivative was only considered as a special damping force.
Now we study the steady-state solution, which is more important and meaningful in vibration rock. For steadystate response, = 0 and = 0, Equation (33) becomes (36) The amplitude and the phase in the steady-state are given by the follow equations: The potential energy stored in the rocks surrounding fault is given by: The maximum energy stored is given as: Equation (37a) is an implicit equation for the amplitude of the response as a function of the frequency and the amplitude of the forcing term ₀: it is called the frequency-response equation. The Newton-Raphson algorithm is used to solve it in order to obtain the amplitude response curves (Ω) which are plotted in Figures 4-12 for three different values of the parameter. , , , , and respectively. In these figures, the response amplitude presents a maximum (infinite response) for Ω = Ω. We recognize such behavior as a resonance phenomenon and one can notice that the parameters , , , , and have a real effect on the amplitude of such oscillatory states.
The behaviour of amplitude-frequency is represented in Figure 4, where the comparison between the analytical and the numerical response frequency of the model is shown. We notice that both curves have approximately the same profiles; therefore; we can conclude that the harmonic balance method chose to solve our problem is appropriate.
To analyze the effect of the parameters on the response system, the amplitude-frequency and stored energyfrequency curves are shown in Figures 5 -7 , 9 -11. From the observation of amplitude-frequency and energy-frequency curves in Figure 5, 6, one find that the resonance amplitude (Ω) system decreases when the fractional order and the fractional damping coefficien increase. This is due to the increase of the equivalent linear damping coefficient. This shows that in the specific condition of the earth's crust, the stored energy in the rock can be controlled by the viscosity. From Figure 7, the resonance amplitude (Ω) changes with increasing the fractional -order r. This is very clear in Figure 8b which represents the resonance amplitude versus the fractional-order r. In this figure, the resonance frequency Ω of the fractional-order system decreases with increasing the fractional-order r. This is very clear in Figure 8a which illustrates the resonance frequency Ω. This is due to the decrease of the equivalent linear stiffness coefficient. Figure 9 and Figure 10 show the effect of the fractional power deflection and the fractional power viscous on the amplitude-frequency and the stored energy.
From the observation of the amplitude response and the stored energy-frequency in Figure 9, one finds that the resonance amplitude increases and the resonance frequency decreases with increasing the fractional power . This can be due to the increase of the equivalent nonlinear stiffness and the equivalent nonlinear viscous damping coefficients.
The increasing fractional power , increase the resonance amplitude of the system (Figure10). This may be due to the decrease in the equivalent nonlinear damping coefficient. The effects of the nonlinear stiffness coefficient on the amplitude-frequency are shown in Figure 11 and it appears on this figure that the resonance frequency increases and the resonance amplitude decreases with increasing the nonlinear stiffness coefficient. The amplitude decrease can be for the reason that the equivalent nonlinear damping coefficient increases with the nonlinear stiffness coefficient.

Christian Foka Fogang et al.
10     Figure 12 illustrates the effect of the resonance frequency Ω (Figure 12 a) and the resonance amplitude (Figure 12 b).
The amplitude response curves ( ) and ( ) are plotted in Figure 13. The resonance phenomenon is observed in this figure.
In the earthquake analogy; the resonance phenomenon can be explained by the fact that the rock of earth crust has stored much energy which created the news fissures that involve news small faults. In the sandblaster zone, the water can up flow through this fissure to provoke the liquefaction phenomena.

Christian Foka Fogang et al.
12 Figure 7. Effect of the fractional-order (damping term) on the amplitude response of the system (a) on the maximum energy (b) as a function of the excitation frequency. Other parameters are as given in Figure 5. It appears that increasing the fractional-order (damping term) results a variation in the resonance amplitude of vibration and decrease in the resonance frequency.

Numerical results
There are many numerical methods for fractional-order differential systems, including computation in time or frequency domains. The numerical formula for the fractional derivative used here is [Petras, 2011;Podlubny, 1999;Tanekou et al., 2018] (40) In this is the memory length, ℎ is the time integration step, = ℎ and ⁽⁾ are the binomial coefficients satisfying the following relations: (41) Figure 9. Amplitude response of the system (a) and stored energy (b) as a function of the excitation frequency. Other parameters are as given in Figure 5. It appears that increasing the fractional power deflexion results a increase in the amplitude of vibration and an decrease in the resonance frequency.    Therefore, the numerical solution of the fractional differential equation 0 ( ) = ( , ) using Equation (40) is finally given by the following equation: (42) Based on the above investigations, the equation of motion Equation (7) can be decomposed into a set of lower degree equations: For simulation purposes, we will use a numerical solution of Equation (43) obtained by using Equation (42), which leads to the equations in the following form: In which, = 1,2,3,... = /ℎ, is a simulation time and the initial conditions are: (0), (0). By combining the two first equations of the above system Equation (44), we obtain: Equation (43) has been solved numerically with initial conditions ( (0), (0)) = (0,0)and the step integration ℎ = Δ = 0.001. 14, 15, and 16 show the numerical simulations of the temporal variations in slip displacement and slip velocity of the slider in which ⁻¹ = 0.39, = 0.5, = 3.0, 0 0.6, Ω = 0.9, = 3.0. It appears on these figures that system exhibited stick-slip motion (slider moves with "adhesion" to the rough surface).

Figures
In Figure 14 we see that the amplitude of slip displacement and slip velocity in the integer model are smaller than those in the fractional model, this can be because decrease the fractional-order derivative can increase the block acceleration, thus increasing the stored energy in the rock as mentioned in Figure 5. We also remark that decreasing the fractional-order derivative increases the recurrence of the earthquake and reduce the recurrence time of the earthquake. This is because frictional force can decrease with decreasing fractional-order derivative. The displacement along a fault is controlled by the fault rheology [Bizzarri, 2009;Rathbun and Marone, 2010].
Fractional derivative can have some interesting consequences, including one or more repeats slip events or "aftershocks" (Figure 14), these "aftershocks" are caused by re-rupturing the same surface instead of stress transfer to an adjoining part of the fault.
In our model, the introduced fractional-order in the weakening friction should play the same role as the parameter of the fault rheology. As specified above, fractional-order derivative replaces the state variable incorporated in the Dieterich-Ruina friction law [Ruina 1983] which is usually interpreted as a function of the sliding history and the fault's age [Pomeau and Le Berre, 2011]. The decrease of recurrence time with decreasing fractionalorder could have significant implications. In the earthquake analogy, the increase of recurrence slip and decrease of recurrence time might indicate that in the specific conditions of the earth's crust, it is more difficult to produce large earthquakes along a fault when it becomes more mature. This implicates that seismic hazard is higher for a young fault than a mature one.
Increasing the fractional power deflection diminishes the values of the slip size of displacement and the velocity of the slider (Figure 15). This can be for the reason that increase the fractional power decreases the potential energy stored in the rock's fault. This result is consistent with this presented on the (Figure 8). Figure 16 that the values of the displacement and the velocity with high viscosity are smaller than those with low viscosity because viscosity can decrease the amplitude of vibration. The viscous effect on slip and the velocity of the model is consistent with that mentioned by the system's response (Figure 7). Also, the presence of viscosity results in smaller velocities of sliders than the absence of viscosity. We can also see that the loading time of the earthquake decreases with increasing viscosity and the duration time of an event increases. This suggests that viscosity can cause slow earthquakes or creep of faults. The viscous effect on the model is consistent with that mentioned by Wang [2007] who applied the model to study the difference in ground motions between the northern and southern segments of the Chelungpu fault generated by the 1999 Chi-Chi, Taiwan, Ms 7.6 earthquake.  From the observation of Figure 19, one could remark the fractional-order affects the slip size of fault ( Figure   19a) and the recurrence of an earthquake (Figure 19b) considerably.

It appears in
We remark that increasing the viscosity decreases the stored energy of the rocks. This show that increasing viscosity can cause oscillation death, therefore, tends to suppress the seismic activity, high viscosity can produce non-volcanic tremors, silent slip or steady slip. Viscosity is responsible for the reduction of seismic radiation during seismic ruptures [Daub et al., 2011;Wang, 2017b], this behavior is most often observed in the lower crust "plastosphere" where rocks are deformed by flow and also in the fault segments that incorporate the clayey rocks.
In these zones, slip is controlled by plastic, ductile or viscous yielding. The fractional viscous model can be used to model the earthquakes dynamics within the lower crust or within the upper mantle where rocks are deformed permanently without energy accumulation and the fractional power model can be used to model the rigid rocks of the upper crust. Our result can find its application in civil engineering: indeed, to build at the level of the faults and volcanoes active, one must have an idea of the fault behaviour in order to construct by taking into account the earthquake-resistant conditions. We know that when the oscillation frequency of the ground is near that of a building, there is resonance.

Conclusion
The dynamics of an improved one-dimension nonlinear spring-block model for earthquake with fractional viscous damping subjected to the strengths due to the motion of the tectonic plates and the up flow of magma is studied. The approximately analytical solution is obtained by the harmonic balance method. In the literature, the resonance amplitude is only affected by the linear and nonlinear damping coefficients, and the resonance frequency is only affected by the linear and nonlinear stiffness, contrarily, in this work the resonance amplitude and the resonance frequency are also affected by fractional-order derivative and the fractional power. It is because the fractional-order parameters can affect the linear and nonlinear equivalent damping coefficient and linear and nonlinear equivalent stiffness coefficient simultaneously. When ≠ 0, ≠ 0 and ≠ 0, fractional-order derivative and will possess both the function of stiffness and the function of damping simultaneous.
Our results show that resonant amplitude of rock vibration decreases with the increase of the order of the time derivative, the fractional viscous coefficient, fractional power deflection and the nonlinear stiffness coefficient and increases with increasing fractional power of viscous damping. Note that the resonant frequency decreases with the increase of the order of the time derivative and increases with the increasing the fractional power deflection and the nonlinear stiffness coefficient. This resonance is responsible for the transition from creep to stick-slip-like behaviour [Perfettini and Schmittbuhl, 2001]. In our model, the resonance case (Ω=Ω) is more interesting since stick-slip is observed, showing that the fractional-order parameters variations can destabilize the creeping fault.
The model gives us an excellent interpretation of the earthquake as a stick-slip motion. We have shown that the recurrence time of the earthquake and the slip size of the fault can be controlled by the fractional-order, the fractional viscous coefficient, fractional power deflection and the magma strengths. These results may be hugely crucial in a dynamic system with the fractional-order derivative and power-law. Our model can be very useful for the system design and control of the dynamical response of a structure subjected to the ground vibrations by choosing suitable fractional-order parameters.