Control instabilities and incite slow-slip in generalized Burridge-Knopoff models
CControl instabilities and incite slow-slip in generalizedBurridge-Knopoff models
Ioannis Stefanou Ecole Centrale de Nantes, Universit´e de Nantes, CNRSGeM (Institut de Recherche en G´enie Civil et M´ecanique), Nantes, France
Corresponding author: Ioannis Stefanou, [email protected] –1– a r X i v : . [ phy s i c s . g e o - ph ] A ug bstract Generalized Burridge-Knopoff (GBK) models display rich dynamics, characterized byinstabilities and multiple bifurcations. GBK models consist of interconnected massesthat can slide on a rough surface under friction. All masses are connected to a plate,which slowly provides energy to the system. The system displays long periods of qui-escence, interrupted by fast, dynamic events (avalanches) of energy relaxation. Duringthese events, clusters of blocks slide abruptly, simulating seismic slip and earthquakerupture.Here we propose a theory for preventing GBK avalanches, control its dynamicsand incite slow-slip. We exploit the dependence of friction on pressure and use it as abackdoor for altering the dynamics of the system. We use the mathematical Theoryof Control and, for the first time, we succeed in (a) stabilizing and restricting chaosin GBK models, (b) guaranteeing slow frictional dissipation and (c) tuning the GBKsystem toward desirable global asymptotic equilibria of lower energy. Our controlapproach is robust and does not require exact knowledge of the frictional behavior ofthe system. Finally, GBK models are known to present Self-Organized Critical (SOC)behavior. Therefore, the presented methodology shows an additional example of SOCControl (SOCC).Given that the dynamics of GBK models show many analogies with earthquakes,we expect to inspire earthquake mitigation strategies regarding anthropogenic and/ornatural seismicity. In a wider perspective, our control approach could be used forimproving understanding of cascade failures in complex systems in geophysics, ac-cess hidden characteristics and improve their predictability by controlling their spatio-temporal behavior in real-time.
Generalized Burridge-Knopoff models (GBK) (Burridge & Knopoff, 1967) mod-els consist of several interconnected masses that can slide on a rough surface underfriction. The sliding blocks are connected to a driver plate that slowly transfers loadand energy to the system. This system is characterized by slow and fast dynamics andshows spatio-temporal correlations due to the high interconnectivity of the masses.The continuous slow movement of the driver plate leads to instabilities expressed asabrupt sliding of clusters of blocks, which can be avoided following the proposed the-ory. GBK models are also frequently used as simplified analog models for earthquakes,as they combine the basic mechanisms of Reid’s elastic rebound theory (Reid, 1910)and statistical similarities.The main objective of this work is to show that the chaotic, rich dynamics ofgeneralized Burridge-Knopoff models can be altered and controlled. Then, slow-slipand smooth energy relaxation can be achieved in a controlled and designed manner.Our theory is inspired by the recent experience that humans do cause earthquakesby injecting fluids under pressure in the earth’s crust leading to fault reactivation(Raleigh et al., 1976; McGarr et al., 2002; Keranen et al., 2013; Rubinstein & Mahani,2015; Cornet, 2016; Guglielmi et al., 2015; Cappa et al., 2019). As a result of fluidpressure change, the spatio-temporal distribution of earthquakes can be perturbed (e.g.Petersen et al., 2015). However, until now, this is done in an uncontrollable manner,which results in significant criticism of active projects involving fluid injections in theearth’s crust (see conventional and unconventional energy production in oil and gasindustry, renewable energies like deep geothermal projects, CO sequestration etc.).Here, the problem of fluid injections is seen from another perspective.We exploit the frictional dependence on pressure and the possibility of alteringthe local equilibrium conditions by fluid injections. Consequently, fluid pressure isconsidered as a backdoor for altering the dynamics of the GBK system. The main –2– ngredient of our approach is the modern mathematical Theory of Control (Vardulakis,1991, 2012; Khalil, 2015), which is applied in order to:1. Stabilize the GBK system and2. drive it to lower energy levels and stable equilibria.This is achieved without precise knowledge of the mechanical parameters of the systemand of its heterogeneous and uncertain frictional properties.GBK models are considered as qualitative analogs of a single earthquake fault oras models of distributed seismicity (see Turcotte, 1999, among others). However, it isworth emphasizing that the dynamic behavior of real faults can be much richer thanthat of the Burridge-Knopoff model and of its generalization (cf. Barbot, 2019; Barbotet al., 2012). The main limitations of models of the Burridge-Knopoff type are exten-sively discussed and shown by Rice (1993) and later publications. However, the valueof using this analog model resorts to its simplicity and the fact that its rich dynamicbehavior is well documented and thoroughly discussed in the literature. Moreover,its mathematical structure allows the development of a general control approach thatcould be applied in more realistic situations of earthquake rupture and instability, pro-vided that a consistent discretization approach like the one proposed by Rice (1993);Chinnery (1963); Erickson et al. (2011) is followed (see also Appendix A). Therefore,from the mathematical point of view, more realistic cases could be tackled using thetheoretical developments presented herein, but this extends the scope of the presentwork. However, it is worth mentioning that the present approach has several limi-tations. These are related, among others, to the actual techno-economical means forfluid injections in the earth’s crust, the sampling-rate and frequency of observations,the exclusion of poroelastic effects and the in-situ hydrogeological and geomechanicalconditions. These limitations and their solution is explored in the frame of the ongoingERC project “Controlling earthQuakes - CoQuake” ( http://coquake.eu ).The paper is organized as follows. In Section 2 we present the main ingredientsof generalized Burridge-Knopoff models, we discuss their frictional properties and weshow how the dependence of friction on pressure can be exploited for achieving ro-bust control of the dynamics of the system. In Section 3 we confirm the well-known Self-organized critical (SOC) behavior for GBK models and we give an example ofsystem stabilization using our theoretical developments. Furthermore, we show howthe input pressure can be adjusted in real-time for assuring slow-slip and driving thesystem to lower (potential) energy states (robust tracking). The mathematical proofsfor stabilization and tracking under uncertainties are given in the Appendices. Theresults of our approach are discussed extensively in the last Section of this work, whereperspectives and implications of the current framework are given for SOC control andgeophysics, in general, and for man-made and natural earthquakes, in particular.
We consider an ensemble of n blocks of mass m each one, connected throughsprings and dampers. Each block can slide independently on a rough, horizontal planeas shown in Fig. 1. The blocks are connected together and with a plate through springsand dampers. The plate, which is called here driver plate , is translated under constantvelocity and provides energy to the system, which is progressively stored in the elasticsprings. During this phase the system is in stable equilibrium and its total potentialenergy increases. Due to slip and slip-rate frictional weakening (Fig. 1d,e), at a certainpoint this equilibrium becomes unstable and some blocks slide abruptly. During thisphase, a part of the stored energy is dissipated abruptly due to friction and damping.The dynamics of this system is described by the following set of non-dimensional –3– c δ ∞ ,v ∞ k c k l Driver plate F r1 F r2 F rn δ ∞ , v ∞ k c k l k l k l Driver plate d e p t h S u rf a c e Fault segment
Slip μμ s μ k Slip-rate μμ s (a) (b)(c) (d) (e)slip weakeningfriction slip-rate weakeningfriction Figure 1.
Schematic representation of (a) the Burridge-Knopoff model (b) its two-dimensionalgeneralization and (c) of a strike-slip fault discretized in N x × N z segments ( N x = 4, N z = 6in this Figure). Evolution of friction coefficient with slip (d) and slip-rate (e). The dashpotsconsidered in the mathematical model of the current GBK model are not drawn for the sake ofsimplicity. equations (see Appendix A): ¯ x (cid:48) = ¯¯ G ¯ x + ¯ H (¯ x ) , (1)where ¯ x is the state vector, whose first n components represent the dimensionlessdisplacements and the rest n components the dimensionless velocities of the masses, ( . ) (cid:48) denotes the dimensionless time derivative, the vector ¯ H (¯ x ) represents the forces appliedto the blocks due to the initial deformation of the springs, the displacement of the driverplate and friction. The matrix ¯¯ G contains information on how the blocks are connectedtogether and can describe different physical situations as shown in Fig. 1. Thesephysical situations span from the classical Burridge-Knopoff model, (see Fig. 1a andBurridge & Knopoff, 1967; Dieterich, 1972; Carlson & Langer, 1989, among others), its2D generalization (see Fig. 1b and Ito & Matsuzaki, 1990; Brown et al., 1991; Huanget al., 1992, among others) and a strike-slip fault discretized in N x × N z segments(Fig. 1c). The system is said to be in equilibrium when ¯ x (cid:48) = 0. All the physicalquantities of the system were scaled as described in Appendix A, where more detailsare also given.It is worth emphasizing that several variants of the original Burridge-Knopoffmodel exist in the literature (Ben-Zion, 2008; de Arcangelis et al., 2016, among manyothers). For this purpose, the above formulation is kept general and allows the in-corporation of most of these variants by adapting the connectivity matrix ¯¯ G and theforcing vector ¯ H . This general formulation allows to derive general proofs about thecontrol of the system as shown in Appendix C and Appendix D. These proofs accountfor robustness, meaning that the system can be controlled even in the absence of de-tailed information regarding the connectivity of the blocks (including elasticity andviscosity) and friction. Here we considered only uncertainties related to friction, butthe extension of the theory to cover uncertainties related to elasticity, viscosity andconnectivity is straightforward. –4– .2 Friction and instabilities We assume Coulomb friction, F ri = µ i N ( e ) i , where N ( e ) i is the effective normalforce applied on block i and µ i is the friction coefficient that may depend on slip,slip-rate, time and other internal/state variables.In the frame of fault mechanics, several rheological models have been proposedfor describing the apparent friction of faults and fault gouges (e.g. Byerlee, 1978; Di-eterich, 2013; Scholz, 2019). In tribology as well. Despite the existing wide experimen-tal and theoretical knowledge about friction and given the presence of heterogeneities(cf. Dieterich, 1979; Marone, 1998; Di Toro et al., 2011; Rattez, Stefanou, & Sulem,2018; Rattez, Stefanou, Sulem, Veveakis, & Poulet, 2018; Barras et al., 2019; Rattez& Veveakis, 2020; Kenigsberg et al., 2020; Collins-Craft et al., 2020), friction is nota well constrained quantity and it is hard to quantify without large uncertainties. Tothis extent we keep here a simple, but general rheology for friction and we adopt astatic/dynamic friction law. In particular, the coefficient of friction evolves from aninitial value µ s (static friction coefficient), to a residual one µ k (kinetic friction coeffi-cient). Fig. 1d-c shows schematically the transition between static and kinetic frictionfor a slip and slip-rate weakening law. This transition is made in a characteristic dis-tance D c (Kanamori & Brodsky, 2004). This characteristic distance, together with thefriction coefficient drop and the applied effective normal force, determine the apparentfrictional weakening during sliding and render the system unstable leading to suddenslip (stick-slip motion Ruina, 1983; Stefanou, 2019). Alternatively, the classical rate-and-state friction law could be used, but this model is mathematically singular forvery low (and very high) sliding velocities, as the ones developed herein when the sys-tem is under control, and therefore it is not appropriate without proper regularization(Ben-Zion, 2008).From the physical point of view, a dynamic instability takes place when theelastic unloading of the springs cannot be counterbalanced by friction. The sameholds when the friction shows velocity weakening that cannot be counterbalanced bythe viscosity of the dashpots. Stability analysis using Lyapunov methods can showthe exact conditions for which this system becomes unstable (cf. Stefanou, 2019).When several blocks are considered, the system presents a chaotic, SOC behavior,which makes its dynamical behavior extremely rich and challenging to control (see forinstance Carlson & Langer, 1989; Huang et al., 1992; Schmittbuhl et al., 1996; Becker,2000; Erickson et al., 2011). Key element for the control and arrest of the above mentioned instability is theeffective normal force N ( e ) i , which can give us a valuable access to the dynamics ofthe system (backdoor). By changing N ( e ) i , the apparent friction of the blocks canbe modified in a desired way. For instance, by reducing N ( e ) i , the friction is reducedresulting in a lubrication-like effect and enhancing slip. By increasing N ( e ) i , the ap-parent frictional force is increased and hampers sliding. These mechanisms are centralfor stabilizing the system and driving it smoothly to a desired equilibrium of lowerenergy, using the mathematical theory of control as shown below.For instance, the effective normal force and, therefore, the apparent friction, canbe altered by injecting fluids into the frictional interfaces where sliding takes place. Theinjection of fluids can change the fluid pressure and, therefore, the apparent friction.Adopting Terzaghi’s principle of effective stress (Terzaghi, 1925), N ( e ) i = N i − P i ,where P i is the interstitial fluid pressure change at the interface of block i and N i isa constant, reference normal force (e.g. the weight of the block). As a result, fluidpressure can play the role of input into the dynamical system, described by Eq. (1). –5– n the case of real faults, N i is a fraction of the effective overburden earth load,averaged over the area of the faut’s segment i . Its value depends on the tectonic settingand the type of the fault. P i corresponds to a change in the interstitial fluid pressure,provoked, for example, by fluid injections at the vicinity of the segment i . It is worthmentioning that injecting fluids in the earth’s crust and altering the local equilibriumby fluid pressure changes is nowadays a common practice in many industrial projects.Some examples are deep geothermal projects, CO sequestration and the oil industry(Rubinstein & Mahani, 2015). However, recent experience shows (e.g. Keranen et al.,2013; Cornet, 2016; Grigoli et al., 2018; Kwiatek et al., 2019; Cornet, 2019; Parisio etal., 2019; Hofmann et al., 2019; Cauchie et al., 2020) that, in many cases, earthquakesare triggered due to the reactivation of tectonic faults in the earth’s crust, which is in astate of marginal stability. Consequelntly, the results of this work could be used for thediscovery of systematic strategies for controlling unwanted seismicity and controllingthe complex dynamics that faults and earthquakes display (Turcotte, 1999). It shouldbe mentioned, though, that for the sake of simplicity poroelastic effects are neglectedin this work, despite their important role in fault reactivation and induced seismicity(Segall & Lu, 2015). Their control will be explored in more practical applications ofthe proposed framework in the future. Generalized Burridge-Knopoff models have rich dynamics and show chaotic be-havior. The are also frequently used in statistical physics and geophysics as paradigmsof criticality and self-organized criticality. In this Section we first illustrate the richdynamics of the system by showing its SOC behavior. Then we show how its dynamicscan be altered by using the general mathematical developments presented in AppendixC and Appendix D, which are based on the mathematical theory of control. The nu-merical examples show how the system can be stabilized and how it can be driven tostable equilibria of lower (potential) energy in a controlled way avoiding instabilitiesand cascade phenomena.
Many systems in nature show universality in their behavior. A large class ofthem is believed to exhibit self-organized criticallity . These systems are continuouslyin or close to a state of marginal stability , show chaotic behavior and obey to similarspatio-temporal correlations and statistical laws. Some examples of SOC are believedto be earthquakes, climate fluctuations, forest and wild fires, snow avalanches, rice-and sand-piles, traffic flows, power electric grids, living organisms (see
Game of Life (Gardner, 1969)), population dynamics and evolution, brain neural activity and sparks,stock markets, wars and pandemics (see Jensen, 1998; Turcotte, 1999; Watkins et al.,2016, and references therein).The term
Self-Organized Criticality (SOC) was coined in 1988, in the seminalpaper of Bak et al. (1988) (see also Bak & Chen, 1989), in order to describe theemergence of a critical state in dissipative, dynamical systems, which have no intrinsictime or length scale. The analysis of Bak et al. is based on a cellular automaton(Wolfram, 1983), in which a particle is added to a randomly selected cell in a squaregrid of cells. When a cell in the grid accumulates four particles, the particles areredistributed to their neighboring cells or they are lost (deleted), if they exceed thegrid. This conceptually simple model leads to a behavior characterized by long peri-ods of stasis (quiescence) interrupted by intermittent bursts of activity involving theavalanche of few or many particles. These instabilities follow a frequency-area power –6– fractal) distribution: N ∝ N − af , (2)where N is the number of avalanches, A the area, i.e. the number of particles involvedin the avalanche, and a ≈ u ∞ ≈ u ∞ . Eachstep in this plot corresponds to a dynamic event of abrupt sliding of either a single,a cluster or of all the blocks of the system. The magnitude of the observed jumps indisplacement depends on the number of the involved blocks in each event. Before anevent, a period of quiescence is observed. During this period, the energy is progressivelystored in the springs and no slip takes place (see plateaus in Fig. 2a). Then the systembecomes again unstable and sudden slip occurs, as previously explained. In this sensethe system is always in a state of marginal stability.Figures 2b-d depict slip events involving a single block and clusters of blocksof various sizes. Often, large events start with the sliding of only one or a couple ofblocks, which push their neighbors to sliding in a similar way to a chain reaction (seeFigures 2c-d). The reported slip velocities are high, compared to the slow time scaleof the movement of the driver plate. It is worth mentioning that similar behavioris observed in systems with more blocks, in different configurations (cf. Fig. 1) andfor real earthquakes. In the case of earthquakes, large events correspond to the mainseismic event, while smaller ones to foreshocks and aftershocks (Scholz, 2019).As shown in Fig. 3, the frequency - number of blocks involved in a slip event( N f ) is found to satisfy the power law distribution with a ≈ .
5. Similar exponentvalues were found for different frictional laws and for larger systems of blocks withhigher interconnectivity (see Carlson & Langer, 1989; Turcotte, 1999; Huang et al.,1992; Brown et al., 1991; Narkounskaia et al., 1992; Rundle & Brown, 1991, amongothers). The above power-law distribution can be also transformed to a frequency-event magnitude distribution, which is more common in seismology see Huang et al.(1992). Notice, that the observed divergence from the power-law for events of highersizes is due to the finite size of the system (here 24 blocks) and it is commonly observedin this kind of simulations in the relevant literature (e.g. Huang et al., 1992). In thiscase, the necessary conditions for SOC behavior, which are detailed in the Discussion(Section 4.1), do not hold and the system is not representative of SOC anymore.Nevertheless, our control strategy is independent of SOC manifestation and is alwaysworking independently of the size of the events and the statistics of the uncontrolledsystem. –7– uiescence cascadeof events single event involving 12 blockssingle event involving 4 blocks quiescence quiescence (a) (b) (c)
12 blocks (d)
Figure 2. (a) Evolution of average dimensionless accumulated slip, (cid:104) ˆ u i (cid:105) = n (cid:80) ni =1 ˆ u i , infunction of the dimensionless driver plate displacement. The jumps correspond to fast, dynamicevents of slip involving one or several blocks (avalanches). Cascade or single events are precededby large periods of quiescence (plateaus), where energy is accumulated into the system due tothe slow movement of the driver plate. (b-d) examples of the sudden, unstable sliding of a singleblock, of four blocks and of twelve blocks (out of twenty four). Slip, ˆ u i , and slip-rates, ˆ v i , arereported. –8– α Figure 3.
Power-law distribution of the frequency ( NN ) - number of blocks involved in slipevents ( N f ) / avalanches. The power law exponent is equal to a ≈ . Based on the above statistics and the discussion in Section 4, the presentednumerical example exhibits chaotic, SOC behavior and can be used for illustrating theefficiency of our control approach. The challenge is therefore to extend this periodof quiescence as long as possible and avoid the abrupt energy releases due to suddensliding. Exploiting the analogy of our system with earthquake faults (see for instanceRice, 1993; Bak & Tang, 1989; Turcotte, 1999, Appendix A and Figure 1 and FigureA1), this would mean that the earthquake instability could be prevented (in theory).
The dynamics of the generalized multiblock system, presented herein, can becontrolled using the mathematical theory of control (Vardulakis, 1991, 2012; Khalil,2015). The target is to update the input, which here is the pressure of fluids injected(added) or pumped (removed) at the frictional interfaces of the blocks, in order tostabilize it, i.e. to avoid abrupt slip and sudden energy release. The term stability isused here in the Lyapunov sense (i.e., the system remains close to its equilibrium stateunder small perturbations from it; for a rigorous mathematical definition of Lyapunovstability we refer to Lyapunov (1892) and Stefanou and Alevizos (2016)). + -- -
Figure 4.
Negative feedback control system Σ( P (cid:96) , C ). Σ( P (cid:96) ) is the GBK system (plant) to becontrolled with the controller Σ( C ). We assume a general negative feedback control system as depicted in Fig. 4. Σ( P (cid:96) )is the multivariable system (plant) to be controlled, i.e. the generalized multiblock –9– ystem in our case, and Σ( C ) is the stabilizing controller we need to design. y ( t ) isthe output of the closed-loop, controlled system Σ( P (cid:96) , C ), which here coincides with¯ x ( t ), ¯ y c ( t ) = ¯ P ( t ) the output of the controller and ¯ x c ( t ) = ¯ x ( t ) + ¯˜ x ( t ) the input of thecontroller Σ( C ), with ¯˜ x ( t ) being a possible perturbation. ¯ x d ( t ) is a desired state of thesystem, such that lim t →∞ ¯ y = ¯ x d . First, we seek the controller Σ( C ) that can immobilize(or stabilize at the origin in terms of Lyapunov stability) the generalized spring-slider(¯ x d ( t ) = 0). Then we will consider specific forms for ¯ x d ( t ) (e.g. constant velocity) inorder to drive the system smoothly to a desired stable equilibrium point and dissipatethe energy in a controlled manner. In the frame of the mathematical theory of control,this process is called tracking .The problem is challenging due to friction and the consequent nonlinearitiesit introduces. Moreover, the exact values of the frictional parameters are usuallyunknown. Our stabilizing controller takes into account this uncertainty and is effectiveeven in the absence of complete knowledge of the system’s parameters (robustness,(Khalil, 2015)). As a result, it manages to unravel in real time the unknown, due touncertainties, dynamics of the system and stabilize it by increasing or decreasing thefluid pressure in the required rate for assuring stability.This is achieved using the mathematical developments detailed in Appendix Cand Appendix D. Our controller guarantees robustness, provided that we can havean estimate a) of the minimum friction coefficient system and b) of the maximumslip and slip-rate softening that can take place. Under this condition, the controllercan stabilize and freeze the system at a desired state and consequently prevent suddensliding. Pressure is adjusted in real time in order to prohibit dynamic events as follows:¯ P = − s ¯¯ B T ¯¯Θ¯ x c , (3)where ¯¯ B and ¯¯Θ are real matrices of constant coefficients and ¯¯ s a diagonal matrixcontaining the sign of the blocks’ velocities. These matrices are defined and determinedin Appendix C. The presence of the diagonal matrix ¯¯ s in the above equation is justifiedby the fact that the friction is always opposed to the velocity and as such it wasincluded to the GBK model for completeness. Consequently, ¯¯ s appears in the controllaw, even though in practice it is expected to be equal to the diagonal matrix in mostcases due to the constant motion of the driver’s plate and the subsequent loading ofthe blocks. Moreover, only tracking under monotonous displacements is of practicalinterest herein, which implies again that the matrix ¯¯ s will coincide with the identitymatrix.The above equation specifies the control law Σ( C ) and assures global asymptoticstability of the closed-loop system Σ( P (cid:96) , C ), through negative feedback. Notice, thatthe above control law is a continuous time regulator and requires continuous monitoringof the state of the system, i.e. of the slip and of the slip-rate. However, in practicalsituations direct access to the exact state may not be possible. This difficulty couldbe bypassed by the derivation of observers (e.g. Khalil, 2015; Vardulakis, 1991),given the observability of our system. Regarding the sampling rate, control can bepossible if the system is transformed and studied as discrete-time one. Depending onthe technological limitations, the sampling rate could be high or not and could be aparameter to optimize. Quantifying the predictability horizon (cf. Gualandi et al.,2020) of the dynamics of the system could be a useful indicator for optimal control indiscrete time.Fig. 5 shows an example of stabilization of the controlled system. In this examplethe controller was activated after the initiation of unstable sliding of the two over thefour blocks presented in Fig. 2c. The controller automatically reduces the fluid pressureat the interfaces of the sliding blocks (Fig. 5b) and immobilizes them as illustrated inFig. 5a. It is worth mentioning that if the controller is activated before the sliding –10– blocks (a)(b) Figure 5.
Stabilization of the unstable movement of blocks as presented in Fig. 2b. Thecontroller was activated at ˆ t = 2 . P i corresponds to no fluid pressure change, positive to pumping of more fluidand negative to fluid withdrawal. event, then the instability is completely avoided by tiny decreases in the pressuresat the four frictional interfaces of the blocks having the tendency to slide (unstableblocks) and their neighbors. The pressure changes in this case are that small thatis not worth of plotting. It is worth emphasizing that this stabilization is achievedwithout knowing exactly the rheology and frictional properties.As far as the pressure can be maintained, the system will be stable and noinstability (sudden sliding and energy relaxation) will take place. However, if thecontroller is deactivated, the system will slide abruptly and its unstable character willbe restored. Therefore, there is a need to drive the system towards a new, stableequilibrium in a stable, smooth, quasi-static way. Once a stabilizing, robust controller has been determined, it can be extendedin order to make the blocks move to a desired new position with a desired velocity.In the mathematical theory of control, this is called tracking (Khalil, 2015). Themathematical extension of the controller is presented in Appendix D, it has the sameform with Eq. (3) and is robust, meaning that it can drive the system in a controlledmanner and in the absence of complete knowledge of its parameters.In order to illustrate how the controller moves the system to a new position,we focus on the avalanche involving twelve unstable blocks, as presented in Fig. 2d.During this event, the total potential energy change was ∆ E ˆ U ≈ − –11– n potential energy happens extremely fast. The maximum velocities reported duringthe movement of the unstable blocks was of the order of ˆ v ≈
1. Two control strategies will be investigated for guarantying at least the same drop in potential energy, butsmoothly, without any unstable, uncontrolled movement of any block.In the first control strategy, the controller will adjust automatically the inputpressure in order to assure translation of all the blocks under the same constant ve-locity, which is chosen to be equal to ˆ v i = 2 × − . This target velocity is threeorders of magnitude lower than the maximum velocity developed during the unstablemovement, but several orders of magnitude higher than the far-field driving velocity,ˆ v ∞ . We allow the system to evolve for a total time equal ˆ t d = 2000. In Fig. 6, wepresent the evolution with time of the displacements and velocities of all the blocks, theinput pressure determined by the controller, the potential energy drop and the energydissipation due to friction. We observe that the controller succeeds in regulating thevelocity of all the blocks to the desired value. A small overshoot in velocities at thebeginning of the controlled sliding is related to the parameters of the controller chosen(see Eq. (D5) in Appendix). This overshoot can be increased or decreased dependingon the desired rate of input pressure change.At the beginning the applied pressure change is zero, but the controller adjusts itautomatically, in order to allow for the blocks to attain the desired velocity, as shownin Fig. 5b. Two groups of blocks are distinguished, i.e. those that are on the vergeof unstable movement and those that are in a stable equilibrium. In the case of theformer, we observe that the input pressure is decreased (negative pressure change; fluidwithdrawal). As a result the regulator increases friction, decelerates the movement ofthe blocks and stabilizes the system. On the contrary, in the case of the latter, thecontroller increases the input pressure in order to accelerate their motion and achievethe desired target velocity. The pressure range (ordinate in Fig. 6b), determinedautomatically by the controller, depends on how far from equilibrium each block is,while the evolution of pressure with time (abscissa in Fig. 6b) on the target velocity,ˆ v i . It is worth emphasizing that the regulator is not based on any “if-then” state-ments - this would be impossible in general situations. On the contrary, based on themathematical developments presented in the appendix, it automatically regulates thefluid pressure of the blocks and stabilizes the system by unraveling its dynamics inreal time. This is accomplished by monitoring the motion of the system and adjustingthe input pressures changes in real-time. Moreover, it has to be mentioned that theregulator is agnostic to the exact frictional parameters and frictional rheology.Regarding the blocks that started from a stable equilibrium, at a certain pointthey also enter to a critical state due to slip accumulation. At that point, the controllerdecreases their input pressure, in order to alter the dynamics of the system, guaranteestability and achieve motion under the desired velocity. This is depicted in Fig. 6b bythe negative pressure change at all the blocks (pressure reduction).The stored energy in the system is dissipated almost linearly with time as shownin Fig. 6c. The same holds for the potential energy decrease in the system. Notice thatthanks to the desired very low target velocities, the kinetic energy of the system andthe viscous dissipation are negligible compared to the drop in potential energy. Withour approach, we finally manage to dissipate more than the released energy duringthe unstable, abrupt movement of the system, but in a slow, smooth and controlledmanner. Consequently, the system is not anymore in a state of a marginal stabilityand it cannot present a self-organized critical behavior (or richer, see Ben-Zion, 2008).Moreover, as its motion is actively controlled, no chaotic behavior can be observed. Inother words, we managed to completely alter the dynamics of the system in a desiredway. –12– igure 6. First control strategy: time evolution of (a) the displacements, ˆ u i , and velocities,ˆ v i , of all the blocks, (b) the input pressures, ˆ P i , determined by the controller and (c) the totalpotential energy drop, ∆ E ˆ U , and the energy dissipation due to friction, ∆ E ˆ F r . We observe thatthe controller succeeds in regulating the velocity of all the blocks to the desired value, ˆ v d , anddissipate the total potential energy in a controlled and smooth manner, avoiding any instabilities.–13– igure 7. Second control strategy: time evolution of (a) the displacements, ˆ u i , and velocities,ˆ v i , of all the blocks, (b) the input pressures, ˆ P i , determined by the controller and (c) the totalpotential energy drop, ∆ E ˆ U , and the energy dissipation due to friction, ∆ E ˆ F r . We observe thatthe controller succeeds in regulating the average velocity of the blocks to the desired value, ˆ v d .The blocks self-organize to a desired stable equilibrium state. Then, after ˆ t = 1250, the controlleris progressively deactivated. Any instabilities are prevented and the total potential energy isdissipated efficiently, as desired. –14– owever, the results show that the controller has to assure a negative inputpressure in order to prevent the unstable movement of some blocks (see Fig. 6b).These lower pressures than the initial ones have to be maintained for assuring stability.Therefore, if we decided to deactivate the controller, the system would become unstableagain. Nonetheless, it is possible to bring it in a state of stable equilibrium, in whichthe controller is no more necessary and can be deactivated. For this purpose, we followa different scenario and instead of setting the velocity of each one of the blocks to thedesired velocity, we set the average velocity equal to (cid:104) ˆ v i (cid:105) = n (cid:80) ni =1 ˆ v i = 2 × − .In this manner we incite a faster, but controlled movement for the unstable blocksfor which the input pressure is negative. As a result some blocks move faster thanothers, noting higher displacements (Fig. 7a). After some time we observe that theblocks self-tune and slide with identical slip-rates equal to the desired one. Once more,the system’s dynamics were altered and no self-organized critical or chaotic behavioroccurred. On the contrary, we incited the system to self-organize towards a stableequilibrium (see ˆ t > t = 1250 and 1750, as depicted inFig. 7. Notice, that the system remains now in a state of stable equilibrium of a lower(potential) energy level. Of course, if the controller remains inactive, the continuousslow movement of the driver plate will render again the system unstable after a (large)time interval. Therefore, it might be interesting after this point to set the controller’starget velocity equal to the driver plate’s velocity. In this case the regulator willautomatically adjust the fluid pressure and the blocks will follow the movement of thedriver plate sliding continuously in an aseismic way and by small pressure changes.The above numerical examples show that the system is finally stabilized by reduc-ing the pressure at the blocks’ interfaces and not by increasing it. This might be at firstcounter-intuitive, if stabilization is thought as a simple process of adjusting pressure inorder to satisfy the frictional stability condition of the system. For the single spring-slider model, frictional stability is guaranteed when the stiffness of the loading systemis higher than a critical softening stiffness, which depends on the exact frictional pa-rameters and rheology (see Ruina, 1983; Scholz, 2019; Stefanou, 2019). For the GBKmodel, the stability conditions are qualitatively similar. Nevertheless, attempting tostabilize the system by satisfying the frictional stability conditions (provided that thefault is reactivated) would require very high positive fluid pressure changes. For typicalvalues of stiffness and frictional weakening parameters of faults the pore fluid increaseshould be close to the insitu effective stress. Moreover, stabilizing the system in thatway wouldn’t mean that it would be immobilized. On the contrary it would slideaseismically, but uncontrollably. In the best case it would follow the far-field tectonicvelocity. Moreover, very high velocities and instabilities could be developed during thepore pressure increase after fault reactivation.However, with our approach we achieve stabilization with minimal pressurechanges. Additionally, we guarantee controlled sliding, with a desired target velocityprofile, decided by the operator. This can be quite important for industrial projectsinvolving fluid injections in depth. Some preliminary, practical numerical exampleson the base of the simple spring-slider model were given in Stefanou (2019), wherepressure reduction is found again to be needed for stabilization. In these examplestracking was achieved by increasing the effective stress by ∼ –15– imits on the pore pressure increase/decrease could be imposed using techniques fromthe mathematical theory of control. Optimal control could be also designed based onengineering criteria for specific applications.As far as it concerns fluid injections in real scale applications, there are severalexamples in the literature correlating induced seismicity with fluid pressure increase(see above cited works, among many others). Of course, correlation is not causality,especially in complex systems as the one at hand and several interpretations mightexist. Our numerical results show that fluid pressure increase leads to slip accelerationand fluid pressure decrease to slip arrest. This is justified from the physics point ofview (decrease vs increase of effective stress and therefore of friction). Moreover, itseems to be corroborated with fluid injections and production at the low temperaturegeothermal field Laugaland `ı Holtum in the south Iceland Seismic zone. According toFl´ovenz et al. (2015) the decrease of the fluid pressure due to geothermal production(fluid withdrawal) and seasonal variation in the pressure seem to have modulated thenatural seismicity by delaying an impending M w = 6 . SOC could be seen as a spectacular manifestation of order in nature that resultsin sparks of energy relaxation (dissipation). Nonetheless, this does not mean that SOCbehavior cannot be prevented. The necessary and sufficient conditions for SOC thatwere recently proposed by Watkins et al. (2016) and are presented below, leave openthis possibility.
Self-Organized Criticality Control (SOCC) can be of particular importance inmany situations where avalanches due to SOC behavior are unwelcome. SOCC is arelatively new field. Maybe the most popular example of SOCC is the prevention oflarge snow avalanches, by triggering smaller ones (McClung & Schaerer, 1993; Birke-land & Landry, 2002). Cajueiro and Andrade (2010a, 2010b, 2010c) applied andextended this idea for controlling self-organized criticallity in the Abelian sand pilemodel (Dhar & Ramaswamy, 1989) and generalizations of it. Brummitt et al. (2012)studied the suppression of cascade failures in interconnected powergrids, based on thesand pile model of Bak et al. (1988). Again using as model the classical sandpileautomaton of Bak et al., No¨el et al. (2013) proposed a control strategy that deter-mines the grid cell in which a particle should land in order to adjust the probabilityof triggered cascades and mitigate large avalanches. Another example of SOC controlis given by Hoffmann and Payton (2014), who altered the SOC power law statisticsof electrical circuits obeying Kirchoff’s law, by adequately modifying the interconnec-tivity of the circuit network. In this way they proposed mitigation strategies of largecascade events.According to Watkins et al. (2016), a system has to satisfy the following three necessary conditions in order to qualify as SOC:NC1. Non-trivial scaling.NC2. Spatio-temporal power law correlations.NC3. Apparent self tuning to the critical point.These necessary conditions are considered in the logical sense. In other words, a systemcannot exhibit SOC if any of the above three conditions is not fulfilled. –16– n extensive discussion of the meaning of critical point and criticality in the frameof SOC and its connection with existing notions in physics and statistical mechanicsis made in Watkins et al. (2016). Here, critical points are points in the phase portraitof the system that are (Lyapunov) unstable. Indeed, due to the slow driver plate’smovement the system is evolving continuously toward a critical, unstable equilibriumfollowed by abrupt cascade events (non-equilibium states in the mathematical sense,see also Figure 2a). These events may be small, involving few blocks or large, involvingseveral blocks.Watkins et al. (2016) give also the following three sufficient conditions for char-acterizing a system as SOC:SC1. Non-linear interaction, normally in the form of thresholds.SC2. Avalanching.SC3. Separation of time scales.That means that if a system fulfills these conditions, then it can exhibit SOC.The open-loop, uncontrolled GBK model, Σ( P (cid:96) ), satisfies all the sufficient con-ditions and therefore exhibits SOC. In particular, friction introduces the necessarynon-linearities and, due to slip or slip-rate weakening, it takes a maximum value beforeslip initiation (threshold). Avalanches are also observed involving clusters of blocksthat dissipate abruptly the energy of the system (intermittent energy relaxation). Thedriver plate’s slow movement introduces a slow time scale in Σ( P (cid:96) ) (see time-scaleasymptotic analysis in Stefanou (2019)), while the events follow the fast characteristictimes related to the frictional instability during avalanches.In contrast, the controlled, closed-loop system, Σ( P (cid:96) , C ), remains strongly non-linear, satisfying condition SC1, but not conditions SC2 and SC3. Moreover, it doesnot satisfy the necessary condition NC3, because our controller, Σ( C ), is designed toprohibit self-tuning to a critical point. Instead, the closed-loop system is self-tunedtoward desirable stable equilibria. As a result, based on the necessary conditionsof Watkins et al. (2016), the controlled system cannot exhibit SOC. Self-organizedcriticality is controlled.It is worth pointing out that our control approach differs from the aforementionedSOCC approaches in many aspects. First of all, it is based on a totally differentmathematical framework. This framework allows us to derive rigorous mathematicalproofs about the stabilization and controllability of the non-linear system. Moreover, itallows to alter its dynamics by considering also the uncertainties of the physical model(robustness) in a deterministic way. Notice that the numerical examples presented inSection 3 are only for illustrating the mathematical findings and not for proving orverifying them. Another different aspect of our approach is related to the underlyingmodel and the chosen input. Even though sandpile systems and networks exhibit SOCand display rich dynamics, they differ from the GBK frictional model considered herein.GBK is described here by a set of Ordinary Differential Equations and it is prone tothe application of the mathematical tools of control theory. Moreover, we don’t triggerany instabilities to dissipate energy, as it is done in the above cited works, and we donot artificially increase locally the energy of the system or change its interconnectivityby using statistical methods. In this sense our approach is deterministic even though itconsiders the uncertainties of the underlying physics. As a result we go beyond existingSOCC approaches by slowing down the dynamics of the GBK model. In this way, wedo not only dissipate the required energy in a controlled manner, but we also breakthe separation between slow and fast dynamics of the system, which is key ingredientfor SOC as stated above. –17– .2 Restriction of chaotic behavior Chaotic behavior is also restricted. While the evolution of the uncontrolled GBKsystem is chaotic and could be predictable only in a statistical sense, due to complexityand chaos, the evolution of the controlled system is not. The presence of the controllerguarantees global asymptotic stability (see Appendix C). Hence no limit cycles orchaos are possible (Strogatz, 1994). Moreover, the controller is robust, meaning thatit succeeds in altering the dynamics of the system even without knowing its exactproperties. Indeed, only some rough boundaries of the frictional parameters are neededin order to drive and control the system as desired.
Following the seminal work of P. Bak and his colleagues in 1988 (Bak et al.,1988), a broad range of natural, technological, biological and social phenomena wereidentified to show self-organized critical behavior. Some scientists will claim that Self-organized criticality is ubiquitous in nature. However, some others will be skepticand will criticize the universality of self-organized critical behavior (see discussion inWatkins et al., 2016). The latter will be based on the limitations of available data andobservations. Further research is needed before adopting one or the other side. For adebate on the applicability of the SOC concept to earthquakes we refer to Ben-Zion(2008); Sornette (1999); Main (1999); Lomnitz-Adler (1993), among others.Nevertheless, what is beyond any doubt, is that real-world systems do presentcascade failures. The failure of a node in a multi-node system (here, of a block of theGBK system) can trigger accelerating feedback and cause the failure of other nodes inthe system, in a domino-like way. The reason is that systems in nature are inherentlycomplex and allow information to spread in several spatio-temporal scales. This occursin a way that it can be difficult or even impossible to grasp and model in details.Here, we present an example which, with a quite high degree of abstraction,shows that cascade failures could be prevented, even if we are agnostic to the detailsof the exact spatio-temporal correlations and dynamics of the underlying system. Inother words, a complete understanding of the physics behind cascade failures and of theinterconnectivity of the system’s nodes is not an absolute requirement for preventingavalanches. This is proved mathematically and illustrated through some examplesfor the strongly non-linear, complex geophysical system at hand, by assuring robustcontrol.Furthermore, having the possibility of controlling the dynamics of a complex sys-tem in a robust way, can give us useful information about its inherent, but practicallyinaccessible properties. Here, we drive smoothly the system to desired stable equilib-rium states, which, though, are a priori unknown. By back-analyzing the evolution ofthe stabilizing input, one could draw real-time conclusions about the interconnectivity,the dynamics and the evolving hidden characteristics of the system. Therefore, theproposed strategy could help in improving the current understanding in some systemsthat complexity makes opaque.An additional implication of the proposed theory has to do with predictability .Predicting the evolution of complex systems exhibiting self-organized critical behavior(or richer) is a challenging, but controversial topic (see Watkins et al., 2016; Ben-Zion,2008; Sornette, 1999, for an overview). Nevertheless the ability to predict, can haveimportant consequences in many disciplines (cf. earthquakes, tectonics, volcanoes).The distribution and frequencies of cascade events of the model presented herein orof more complex ones is a useful statistical correlation. However, it cannot providewith certainty when and how large exactly the next cascade failure will be. Never-theless, if control is possible, as shown here, then prediction is irrelevant. The more –18– e control a system, the less unpredictable it becomes. Of course, one would needsufficient inputs and monitoring for guaranteeing full control (see controllability and observability notions of the mathematical Theory of Control in Khalil, 2015, amongothers). However, even if only partial control is possible (e.g. due to limited area ofintervention and technological constraints), the space of uncontrolled dynamics willbe reduced, which can lead to improved predictability (e.g. increase the predictabilityhorizon (e.g. Gualandi et al., 2020)) and constrain the size of the next cascade failure.
A direct implication of the present work is inevitably related to earthquakes.Burridge-Knopoff models are frequently used as qualitative analogs of the earthquakephenomenon, either at the level of single fault or of complex fault networks. Of course,models of the Burridge-Knopoff type have several limitations as far it concerns therepresentability of the earthquake phenomenon. These limitations are well identifiedin the early work of Rice (1993) and relevant literature. However, the mathematicalstructure of the generalized Burridge-Knopoff models allows the development of a gen-eral control approach that could be applied in more realistic situations of earthquakerupture, provided that a consistent fault discretization approach is followed (see Rice,1993; Chinnery, 1963; Erickson et al., 2011; Ben-Zion, 2008, among others and rele-vant numerical methods in elastodynamics). Therefore, more realistic cases could betackled using the theoretical developments presented herein (see Appendix A, FigureA1 and Figure 1), but this extends the scope of the present work.Moreover, it is worth emphasizing that we show a way of active stabilization ofthe generalized Burridge-Knopoff without knowing its exact properties. This meansthat detailed knowledge of faults’ frictional parameters, which is practically impossibleto acquire in practice, might not be a sine qua non condition. Notice, that friction isconsidered as the cornerstone for understanding earthquake behavior and it is a majorunknown (Erickson et al., 2011). However, our control approach needs minimal andnot precise information about the frictional characteristics of the fault system, whichcan be easily acquired in practice. Based on this limited information, we show howthe system can be driven to a stable state in a totally controlled and aseismic way.Moreover, our approach guarantees aseismic, slow-slip and smooth energy relaxationand does not require the knowledge of the exact current stress state and tectonicsetting. The system is controlled independently of being far or close to its criticalpoints.Without any doubt, claiming that controlling anthropogenic or natural seismicityis possible, based on the analysis presented herein, is a speculation and further researchis needed. Several theoretical and techno-economical investigations have to be pursuedfurther in order to show into what extend man-made or natural earthquakes can beprevented (or the opposite). For example, some direct limitations of the proposedtheory have to do with the actual technological means for fluid injections in the earth’scrust, the sampling-rate and frequency of observations, the in-situ hydrogeologicaland geomechanical conditions and uncertainty quantification. However, the currentwork sets the mathematical and physical framework for inspiring further research oncontrolling induced seismicity and, maybe, at a later phase, on controlling naturalseismicity as well.
Acknowledgements
I would like to thank J.-P. Avouac for his insightful feedback, constructive com-ments and our fruitful discussion regarding the application and the limitations of theproposed methodology and model. I would like also to thank A. Gualandi for hisconstructive criticism and positive feedback. –19– his work was supported by the European Research Council (ERC) under theEuropean Union Horizon 2020 research and innovation program (Grant agreement757848 CoQuake), http://coquake.eu . Appendix A Equations of motion
The equation of motion of block i is written as follows: m i ˙ v i = n (cid:88) j =1 k cij ( u j − u i ) + n (cid:88) j =1 η cij ( v j − v i )+ k li ( u ∞ − u i ) + η li ( v ∞ − v i )+ n (cid:88) j =1 k cij (cid:0) u j − u i (cid:1) − F ri , (A1)where ˙( . ) is the time derivative, u i and v i are, respectively, the slip (displacement) andslip-rate (velocity) of the block i , u ∞ and v ∞ are, respectively, the displacement andvelocity of the driver plate, which represents the far field tectonic velocity in the caseof faults, and u i is the initial displacement of the block i . k stands for stiffness and η for damping coefficients. The superscript ‘ c ’ denotes the springs and dampers betweenthe blocks and the superscript ‘ l ’ the same elements between the blocks and the driverplate. For instance, k cij is the stiffness coefficient of the spring connecting block i withblock j , while η li is the damping coefficient of the dashpot connecting the block ‘ i ’ withthe driver plate. F ri represents the friction of block i with the rough plane and candepend on slip, rate of slip and other internal state variables (see section Friction andinstability , 2.3). Here F ri = F ri ( u i , v i ) (see Fig. 1d,e).Setting ω ∗ = (cid:113) k ∗ m ∗ , where k ∗ and m ∗ are, respectively, a reference stiffness andmass, the above equations take the dimensionless form:ˆ m i ˆ v (cid:48) i = n (cid:88) j =1 ˆ k cij (ˆ u j − ˆ u i ) + n (cid:88) j =1 ζ ˆ η cij (ˆ v j − ˆ v i )+ ˆ k li (ˆ u ∞ − ˆ u i ) + 2 ζ ˆ η li (ˆ v ∞ − ˆ v i )+ n (cid:88) j =1 ˆ k cij (cid:0) ˆ u j − ˆ u i (cid:1) − ˆ F ri , (A2)where ( . ) (cid:48) is the derivative with respect to the dimensionless time ˆ t = ω ∗ t , ζ = η ∗ m ∗ ω ∗ is the damping ratio, ˆ k c,lx = k ∗− k c,lx , ˆ η c,lx = η ∗− η c,lx , η ∗ = 2 ζm ∗ ω ∗ , ˆ u i = D ∗− u i ,ˆ v i = v i D ∗− ω ∗− , D ∗ a reference displacement and ˆ F ri = ( D ∗ k ∗ ) − F ri .Equations (A2) are written in matrix form as follows:¯ x (cid:48) = (cid:20) ¯¯ O ¯¯ I − ¯¯ K − ζ ¯¯ K (cid:21)(cid:124) (cid:123)(cid:122) (cid:125) ¯¯ G ¯ x + (cid:20) ¯ O ¯Ψ (cid:21)(cid:124)(cid:123)(cid:122)(cid:125) ¯ H . (A3)The first n components of the vector ¯ x represent the dimensionless displacementsand the rest n components the dimensionless velocities of the blocks. ¯¯ O and ¯¯ I are,respectively, the zero and identity matrices of size n × n and ¯ O the zero vector ofsize n . ¯¯ K = ¯¯ M − (cid:0) ¯¯ K l − ¯¯ K c (cid:1) , where ¯¯ M is a diagonal matrix containing the dimen-sionless masses of the blocks, { ¯¯ M } ii = ˆ m i = 1, { ¯¯ K c } ij = ˆ k cij and ¯¯ K l a diagonalmatrix with components { ¯¯ K l } ii = ˆ k li . We say that the system is in equilibriumwhen ¯ x (cid:48) = 0. In the above equation we assumed that the dimensionless dampingcoefficients coincide with the dimensionless stiffnesses. This is a reasonable assump-tion in the absence of more detailed data. The vector ¯Ψ represents the dimension-less forces applied to the blocks due to the initial deformation of the springs, ˆ u i , –20– c k l k l k l k c k c k l Figure A1.
Graphical representation of the ( i, j ) components of the connectivity matrix, nor-malized by its maximum diagonal component, for the Burridge-Knopoff model with 24 blocks, itstwo-dimensional generalization with a 4 × × the displacement, ˆ u ∞ , and velocity, ˆ v ∞ of the driver plate, and the friction forces¯ F r : ¯Ψ = ¯¯ M − (cid:0) − ¯ F r − ¯¯ K ¯ U + ¯¯ K l ¯ U ∞ + 2 ζ ¯¯ K l ¯ V ∞ (cid:1) , where { ¯ F r } i = ˆ F ri , { ¯ U } i = ˆ u i , { ¯ U ∞ } i = ˆ u ∞ i and { ¯ V ∞ } i = ˆ v ∞ i .The matrix ¯¯ K is called here connectivity matrix and contains information on howthe blocks are connected together. Various geometrical configurations, such as thosepresented in Fig. 1, can be described by adequately adjusting the components of theconnectivity matrix. Fig. A1 shows graphically the connectivity matrix for the 1DBurridge-Knopoff model, its 2D generalization in a 4 × × Appendix B Dynamic simulations and SOC
We consider that the driver plate is moving under constant velocity, which isseveral orders of magnitude lower than the velocities of the blocks that are developedduring abrupt sliding. Consequently, we can assume that the driver plate remains stillduring the sliding events. This situation is inspired by the far-field earth’s tectonicmovement, which is several of orders of magnitudes lower (some centimeters per year)than the seismic slip velocities developed during earthquakes that can reach up to onemeter per second.For the simulations, we first calculate the minimum displacement of the driverplate that can trigger the sliding of at least one block. In this way we avoid simulatingthe slow-dynamics (Stefanou, 2019), quasi-static behavior of the system and we onlyintegrate numerically the dynamic equations of motion of the system for determiningits fast dynamic, unstable response. After each dynamic event, the system reachesa new equilibrium (local minimum potential energy state). The slip of the blocksis recorded and a random small overshoot in their displacements is considered (seealso Brown et al., 1991; Rundle & Brown, 1991). The random overshoot is not the –21– ame between the blocks and varies from zero to 20% of its slip during the previousevent. The random overshoot embodies several uncertainties of the system related toits elastic parameters, initial conditions and frictional properties, among others.The frictional properties of the blocks can be uniform or randomly chosen froma distribution. Slip or slip-rate softening is required to render the system unstable andlead to SOC behavior. Here, we use slip weakening friction as in Huang et al. (1992);Stefanou (2019). Simulations with slip-rate weakening would give similar results (cf.Carlson & Langer, 1989; Huang & Turcotte, 1992). In particular, the friction coefficientevolves from its static value ( µ s ) to its kinetic one ( µ k ). In the simulations presentedhere µ s = 0 . µ k = 0 .
5. The friction drop occurs in a characteristic distance equalto ˆ D c = .
01. The damping ratio ζ is set equal to one and k ci k lij = 2.Each slip event can involve sliding of a single block, a cluster of some blocks orof all the blocks of the system. After each slip event the friction coefficient of eachblock is restored to its static friction value and a new period of quiescence takes placeas shown in Fig. 2.The simulation procedure is summarized as follows:1. Quiescence period:
Determine ¯ˆ u ∞ that renders the system unstable by solving¯Ψ = ¯0 and set the driver plate displacement equal to ˆ u ∞ = min ∀ i ( { ¯ˆ u ∞ } i )2. Sudden slip event:
Integrate numerically the dynamic equations of motion (A3)to determine slip ˆ u i , until max ∀ i (ˆ v i ) ≤ threshold. The threshold was set equalto 0 .
02, which is much smaller than the maximum velocity of blocks duringunstable sliding (see Fig. 2b-d) .3.
Healing:
Set block velocities equal to zero and update their positions ˆ u i =ˆ u i + ˜ u i , where ˜ u i is a random overshoot as described above. Restore the frictionof coefficient from µ k to µ s .Repeat 1 to 4 and record events. A sequence of N = 10000 events were simulated forcalculating the frequency-size statistics presented in Fig. 3. The simulation of moreevents (20000) lead to almost identical results. Appendix C Robust state feedback stabilization
After some algebra (A3) (or (1)) is written as follows:¯ x (cid:48) = ¯¯ A ¯ x + ¯¯ B ¯ P + ¯Ψ , (C1)where ¯¯ A = ¯¯ G + ¯¯ S , ¯¯ S = (cid:20) ¯¯ O ¯¯ Os u ¯¯ I s v ¯¯ I (cid:21) , s u and s v are respectively the absolute valuesof the minimum slip and slip-rate softening rates of the frictional law (see Fig. C1),¯¯ B = (cid:20) ¯¯ Oµ min ¯¯ I (cid:21) , µ min = min ∀ i, ∀ u i , ∀ v i µ i > + ¯Ψ .¯Ψ = ¯¯ (cid:101) B ¯ P with ¯¯ (cid:101) B = (cid:20) ¯¯ O ¯¯˜ µ (cid:21) , ¯¯˜ µ a diagonal matrix with diagonal elements such as { ¯¯˜ µ } ii = µ i − µ min and ¯Ψ = − (cid:34) ¯¯ O ¯¯ s (cid:16) ¯¯∆ µ ¯ N + ¯¯ s ¯¯ S ¯ x (cid:17)(cid:35) , where ¯¯ s , ¯¯∆ µ are diagonal matriceswith diagonal elements { ¯¯ s } ii = sgn(ˆ v i ), { ¯¯∆ µ } ii = µ i − µ si , respectively, ¯ N a vector,with elements { ¯ N } i = N i , N i is a reference normal (effective) force applied at block i and sgn( . ) is the sign function. Finally, ¯ P = s ¯ P (equivalently the input ¯ P = ¯¯ s ¯ P ),that needs to be determined for assuring asymptotic stability.Let the scalar function V (¯ x ) = ¯ x T ¯¯Θ¯ x > x ∈ X ⊂ R n ( n isthe number of blocks) and V (¯0) = 0. Under these conditions, ¯¯Θ is positive definite , –22– ¯Θ (cid:31) ¯¯0, (Brauer & Nohel, 1969) (or negative definite if ¯ x T ¯¯Θ¯ x < ∀ ¯ x ∈ X ⊂ R n \
0, i.e.¯¯Θ ≺ ¯¯0). Moreover let: P = −
12 ¯¯ B T ¯¯Θ (¯ x + ¯˜ x ) (C2)be the control law. ¯˜ x is a perturbation (see Fig. 4). We search Θ such that the closed-loop system Σ( P (cid:96) , C ) (Fig. 4) can be asymptotically stable at ¯ x = ¯0. According toLyapunov’s stability theorem (see Lyapunov’s Second Method,
Brauer & Nohel, 1969),if there exists V (¯ x ) > V (cid:48) (¯ x ) is strictly negative ∀ ¯ x ∈ X ⊂ R n \ ¯0, thenthe origin of the system, ¯ x = ¯0, is asymptotically stable. If X extends over the wholereal 2 n -dimensional Euclidean space, then the origin is globally asymptotically stable.Differentiating V (¯ x ) with respect to time, using (C1) and the symmetry of ¯¯Θ we obtain: V (cid:48) (¯ x ) = −
12 ¯ x T ¯¯ Q ¯ x −
12 ¯ x T ¯¯Ξ¯˜ x + ¯ x T ¯¯Θ¯Ψ , (C3)where ¯¯Ξ = ¯¯Θ¯¯ B ¯¯ B T ¯¯Θ and ¯¯ Q satisfies the algebraic Riccati equation:¯¯ A T ¯¯Θ + ¯¯Θ¯¯ A − ¯¯Θ¯¯ B ¯¯ B T ¯¯Θ = − ¯¯ Q. (C4)¯¯ Q is selected to be real, symmetric and positive definite (¯¯ Q (cid:31) ¯¯0). In the numericalexamples ¯¯ Q was taken equal to the identity matrix but any other matrix could beselected as well provided that it satisfies the aforementioned conditions. Notice thatin the frame of linear systems, which though is not our case, the matrix ¯¯ Q can beselected optimally in order to satisfy some engineering criteria (see Linear QuadraticControl problem). For a given ¯¯ Q exists a unique positive definite and symmetric ¯¯Θsatisfying the algebraic Riccati equation.Therefore, a sufficient condition for the system to be asymptotically stable ( V (cid:48) (¯ x ) <
0) is the third and second terms of the right hand side of (C3), to be negative or zero,i.e. ∀ ¯ x ∈ X ⊂ R n \ ¯0: ¯ x T ¯¯Θ¯Ψ = ¯ x T ¯¯Θ¯Ψ (cid:124) (cid:123)(cid:122) (cid:125) Ω + ¯ x T ¯¯Θ¯Ψ (cid:124) (cid:123)(cid:122) (cid:125) Ω ≤ −
12 ¯ x T ¯¯Ξ¯˜ x (cid:124) (cid:123)(cid:122) (cid:125) Ω ≤ . (C6)Using (C2) we obtain:Ω = ¯ x T ¯¯Θ¯Ψ = ¯ x T ¯¯Θ (cid:101) ¯¯ B ¯ P = − y T (cid:101) ¯¯ B ¯¯ B T ¯ y, (C7)where we set ¯ y = ¯¯Θ¯ x . ¯ y = ¯0 if and only if ¯ x = ¯0 due to the positive definiteness of¯¯Θ. By simple inspection, the matrix (cid:101) ¯¯ B ¯¯ B T is positive semidefinite , because it is theproduct of diagonal matrices with positive or zero diagonal elements (see definitions ofthe diagonal matrices (cid:101) ¯¯ B and ¯¯ B ). A matrix ¯¯ D is called positive (negative) semidifinite¯¯ D (cid:23) ¯¯0 ( ¯¯ D (cid:22) ¯¯0) if and only if ¯ x T ¯¯ D ¯ x ≥ x T ¯¯ D ¯ x ≤ ∀ ¯ x ∈ X ⊂ R n \ ¯0. Therefore,Ω ≤ x ∈ X ⊂ R n \ ¯0, as ¯ x = (cid:80) nk =1 α k ¯ ω ( k ) , with ω ( k ) being the eigenvector k of ¯¯Θ and α k real coefficients. The eigenvalues, λ ( k ) of ¯¯Θ are all strictly positive, dueto its positive definiteness, i.e. λ ( k ) >
0. Therefore:Ω = ¯ x T ¯¯Θ¯Ψ = n (cid:88) k =1 λ ( k ) ¯Ψ T ¯ x. (C8) –23– igure C1. Schematic representation of maximum slip-rate weakening slope s v ≥ s u ≥
0. Friction is always opposite to slip velocity.
After some algebra, ¯Ψ T ¯ x = (cid:16) − ¯¯∆ µ ¯ N − ¯¯ s ¯¯ S ¯ x (cid:17) T | ¯ v | , (C9)where {| ¯ v |} i = | ˆ v i | and | . | the absolute value. As shown in Fig. C1, s u ≥ s v ≥ − ¯¯∆ µ ¯ N − ¯¯ s ¯¯ S ¯ x is negative. Consequently,Ω ≤ x + ¯˜ x ) T ¯¯Ξ (¯ x + ¯˜ x ) − (¯ x − ¯˜ x ) T ¯¯Ξ (¯ x − ¯˜ x ) ≥ B and has eigenvalues ξ ( i ) ≥ ≤
0, (C6), becomes: n (cid:88) k =1 ξ ( k ) x k ˜ x k ≥ ≤
0. In any other case, the above inequality defines a subset of the Euclidean spacefor ¯ x , in which asymptotically stability is guaranteed.As a result, the closed-loop system has V (cid:48) (¯ x ) < P = − s ¯¯ B T ¯¯Θ¯ x c , (C12)where ¯ x c = ¯ x + ¯˜ x . Appendix D Robust tracking
Let ¯ r d = ¯ r d ( t ) be a vector describing the desired displacements of the blocks. Wewant to minimize the error ¯¯ C ¯ x − ¯ r d , where ¯¯ C = (cid:2) ¯¯ I ¯¯ O (cid:3) .Here we use the approach of integral action (Khalil, 2015) and we augment thesystem with the equation: ¯ e (cid:48) = ¯¯ C ¯ x − ¯ r d . (D1) –24– et also ¯ x eq , ¯ e eq and ¯ P eq be such that:¯0 = ¯¯ A ¯ x eq + ¯¯ B ¯ P eq + ¯Ψ eq ¯0 = ¯¯ C ¯ x eq − ¯ r d . (D2)¯ x eq , ¯ e eq and ¯ P eq exist for the physical system at hand.Applying the transformation ¯ z = ¯ x − ¯ x eq , ¯ ξ = ¯ e − ¯ e eq , Eqs. (C1, D1) become:¯ w (cid:48) = ¯¯ A a ¯ w + ¯¯ B a ∆¯ P + ¯Ψ a , (D3)where ¯ w = (cid:20) ¯ z ¯ ξ (cid:21) , ¯¯ A a = (cid:20) ¯¯ A ¯¯ O ¯¯ C ¯¯ O (cid:21) , ¯¯ B a = (cid:20) ¯¯ B ¯¯ O (cid:21) , ¯Ψ a = (cid:20) ¯Ψ − ¯Ψ eq ¯ O (cid:21) and ∆¯ P = ¯ P − ¯ P eq .The above system has the same form with (C1) and the same analysis as abovecan be carried out for determining the controller that assures its stabilization, leadingto: ¯∆ P = − s ¯¯ B Ta ¯¯Θ a ¯ w (D4)or, equivalently, ¯ P = − s ¯¯ B Ta ¯¯Θ a ¯ x a , (D5)where ¯ x a = (cid:20) ¯ x ¯ e (cid:21) and ¯¯Θ a is the solution of the Riccati equation. (C4), but for ¯¯ A a ,¯¯ B a and ¯¯ Q a instead of ¯¯ A , ¯¯ B and ¯¯ Q , respectively. Similarly, ¯¯ Q a is selected to be real,symmetric and positive definite (¯¯ Q a (cid:31) ¯¯0).By regulating the input pressure as in (D5), the system can be driven to thedesired position ensuring robust tracking. Indeed, due to asymptotic stability of (D3)under (D4) or (D5), lim t →∞ ¯ ξ (cid:48) = lim t →∞ ¯ e (cid:48) = ¯0 and, therefore, lim t →∞ ¯¯ C ¯ x = ¯ r d . Appendix E Numerical implementation
The numerical integration of the ordinary differential equations presented hereinwas performed using
SciPy (Virtanen et al., 2020) and the
LSDOA implicit algorithm(Hindmarsh, 1983; Petzold, 1983). The designed controller was programmed in
Python3 (van Rossum, 1995) and the algebraic Riccati equation was solved using the
PythonControl Systems Library ( Python Control Systems Library , 2020). The programs areavailable in
Jupyter notebooks (Kluyver et al., 2016) upon request.
References
Bak, P., & Chen, K. (1989, sep). The physics of fractals.
Physica D: Nonlinear Phe-nomena , (1-3), 5–12. Retrieved from https://linkinghub.elsevier.com/retrieve/pii/0167278989901668 doi: 10.1016/0167-2789(89)90166-8Bak, P., & Tang, C. (1989, nov). Earthquakes as a self-organized critical phe-nomenon. Journal of Geophysical Research: Solid Earth , (B11), 15635–15637. Retrieved from http://doi.wiley.com/10.1029/JB094iB11p15635 doi: 10.1029/JB094iB11p15635Bak, P., Tang, C., & Wiesenfeld, K. (1988). Self-organized criticality. Physical Re-view A , (1), 364–374. doi: 10.1103/PhysRevA.38.364Barbot, S. D. (2019). Slow-slip, slow earthquakes, period-two cycles, full andpartial ruptures, and deterministic chaos in a single asperity fault. Tectono-physics , (August), 228171. Retrieved from https://doi.org/10.1016/j.tecto.2019.228171 doi: 10.1016/j.tecto.2019.228171Barbot, S. D., Lapusta, N., & Avouac, J.-p. (2012). Under the Hood of the Earth-quake. Science , (May), 707–710. doi: 10.1126/science.1218796 –25– arras, F., Aldam, M., Roch, T., Brener, E. A., Bouchbinder, E., & Molinari, J. F.(2019). Emergence of Cracklike Behavior of Frictional Rupture: The Originof Stress Drops. Physical Review X , (4), 41043. Retrieved from https://doi.org/10.1103/PhysRevX.9.041043 doi: 10.1103/PhysRevX.9.041043Becker, T. W. (2000). Deterministic chaos in two state-variable friction sliders andthe effect of elastic interactions. In Geophysical monograph series (Vol. 120,pp. 5–26). Retrieved from doi: 10.1029/GM120p0005Ben-Zion, Y. (2008). Collective Behavior of Earthquakes and Faults.
Reviews ofGeophysics , , 1–70. doi: 10.1029/2008RG000260.1.INTRODUCTIONBirkeland, K. W., & Landry, C. C. (2002). Power-laws and snow avalanches. Geo-physical Research Letters , (11), 49–1–49–3. doi: 10.1029/2001GL014623Brauer, F., & Nohel, J. (1969). The Qualitative Theory of Ordinary DifferentialEquations: An Introduction . New York: Dover Publications.Brown, S. R., Scholz, C. H., & Rundle, J. B. (1991, feb). A simplified spring-blockmodel of earthquakes.
Geophysical Research Letters , (2), 215–218. Retrievedfrom http://doi.wiley.com/10.1029/91GL00210 doi: 10.1029/91GL00210Brummitt, C. D., D’Souza, R. M., & Leicht, E. A. (2012). Suppressing cascades ofload in interdependent networks. Proceedings of the National Academy of Sci-ences , (12), E680–E689. doi: 10.1073/pnas.1110586109Burridge, R., & Knopoff, L. (1967). Model and theoretical seismicity. Bulletin of theSeismological Society of America (1967) , (3), 341–371.Byerlee, J. D. (1978). Friction of rocks. Pure and Applied Geophysics PAGEOPH , (4-5), 615–626. doi: 10.1007/BF00876528Cajueiro, D. O., & Andrade, R. F. (2010a). Controlling self-organized criticalityin complex networks. European Physical Journal B , (2), 291–296. doi: 10.1140/epjb/e2010-00229-8Cajueiro, D. O., & Andrade, R. F. (2010b). Controlling self-organized criticalityin sandpile models. Physical Review E - Statistical, Nonlinear, and Soft MatterPhysics , (1), 1–10. doi: 10.1103/PhysRevE.81.015102Cajueiro, D. O., & Andrade, R. F. (2010c). Dynamical programming approachfor controlling the directed Abelian Dhar-Ramaswamy model. Physical ReviewE - Statistical, Nonlinear, and Soft Matter Physics , (3), 1–12. doi: 10.1103/PhysRevE.82.031108Cappa, F., Scuderi, M. M., Collettini, C., Guglielmi, Y., & Avouac, J.-P. (2019,mar). Stabilization of fault slip by fluid injection in the laboratory andin situ. Science Advances , (3), eaau4065. Retrieved from http://advances.sciencemag.org/http://advances.sciencemag.org/lookup/doi/10.1126/sciadv.aau4065 doi: 10.1126/sciadv.aau4065Carlson, J. M., & Langer, J. S. (1989, dec). Mechanical model of an earthquakefault. Physical Review A , (11), 6470–6484. Retrieved from papers3://publication/uuid/81EFCE69-2B62-4348-9999-B37488DC1AAEhttps://link.aps.org/doi/10.1103/PhysRevA.40.6470 doi: 10.1103/PhysRevA.40.6470Cauchie, L., Lenglin´e, O., & Schmittbuhl, J. (2020, may). Seismic asperity sizeevolution during fluid injection: case study of the 1993 Soultz-sous-Forˆetsinjection. Geophysical Journal International , (2), 968–980. Retrievedfrom https://academic.oup.com/gji/article/221/2/968/5721255 doi:10.1093/gji/ggaa051Chinnery, M. (1963). The stress changes that accompany strike-slip faulting. Bul-letin of the Seismological Society of America , (5), 921–932.Collins-Craft, N. A., Stefanou, I., Sulem, J., & Einav, I. (2020, may). A CosseratBreakage Mechanics model for brittle granular media. Journal of the Mechan-ics and Physics of Solids , (0), 103975. Retrieved from https://linkinghub.elsevier.com/retrieve/pii/S0022509620302106 doi: 10.1016/j.jmps.2020.103975 –26– ornet, F. H. (2016). Seismic and aseismic motions generated by fluid injec-tions. Geomechanics for Energy and the Environment , , 42–54. Re-trieved from http://dx.doi.org/10.1016/j.gete.2015.12.003 doi:10.1016/j.gete.2015.12.003Cornet, F. H. (2019, oct). The engineering of safe hydraulic stimulations for EGSdevelopment in hot crystalline rock masses. Geomechanics for Energy and theEnvironment , (1), 1. Retrieved from https://linkinghub.elsevier.com/retrieve/pii/S2352380819300280 doi: 10.1016/j.gete.2019.100151de Arcangelis, L., Godano, C., Grasso, J. R., & Lippiello, E. (2016). Statisticalphysics approach to earthquake occurrence and forecasting. Physics Reports , (March), 1–91. Retrieved from http://dx.doi.org/10.1016/j.physrep.2016.03.002 doi: 10.1016/j.physrep.2016.03.002Dhar, D., & Ramaswamy, R. (1989, oct). Exactly solved model of self-organized crit-ical phenomena. Physical Review Letters , (16), 1659–1662. Retrieved from https://link.aps.org/doi/10.1103/PhysRevLett.63.1659 doi: 10.1103/PhysRevLett.63.1659Di Toro, G., Han, R., Hirose, T., De Paola, N., Nielsen, S., Mizoguchi, K., . . . Shi-mamoto, T. (2011). Fault lubrication during earthquakes. Nature , (7339),494–498. doi: 10.1038/nature09838Dieterich, J. H. (1972, jul). Time-dependent friction as a possible mechanismfor aftershocks. Journal of Geophysical Research , (20), 3771–3781. Re-trieved from http://doi.wiley.com/10.1029/JB077i020p03771 doi:10.1029/JB077i020p03771Dieterich, J. H. (1979). Modeling of rock friction: 1. Experimental results andconstitutive equations. Journal of Geophysical Research , (B5), 2161.Retrieved from http://doi.wiley.com/10.1029/JB084iB05p02161 doi:10.1029/JB084iB05p02161Dieterich, J. H. (2013, mar). Constitutive Properties of Faults With Sim-ulated Gouge. In Mechanical behavior of crustal rocks (Vol. 24, pp.103–120). Retrieved from http://onlinelibrary.wiley.com/store/10.1029/GM024p0103/asset/ch8.pdf?v=1{\&}t=i5e5o4tp{\&}s=fb2dfdc01a8ff33cbc56a49a81bc18e23eca16c3http://doi.wiley.com/10.1029/GM024p0103 doi: 10.1029/GM024p0103Erickson, B. A., Birnir, B., & Lavall´ee, D. (2011). Periodicity, chaos and lo-calization in a Burridge-Knopoff model of an earthquake with rate-and-state friction.
Geophysical Journal International , (1), 178–198. doi:10.1111/j.1365-246X.2011.05123.xFl´ovenz, ´O. G., ´Ag´ustsson, K., ´Arni Gudnason, E., & Kristj´ansd´ottir, S. (2015).Reinjection and Induced Seismicity in Geothermal Fields in Iceland. Proceed-ings World Geothermal Congress (April), 1–15.Gardner, M. (1969, jan). Mathematical Games.
Scientific American , (1),116–120. Retrieved from doi: 10.1038/scientificamerican0169-116Grigoli, F., Cesca, S., Rinaldi, A. P., Manconi, A., L´opez-Comino, J. A., Clinton,J. F., . . . Wiemer, S. (2018, jun). The November 2017 M w 5.5 Pohangearthquake: A possible case of induced seismicity in South Korea. Science , (6392), 1003–1006. Retrieved from doi: 10.1126/science.aat2010Gualandi, A., Avouac, J., Michel, S., & Faranda, D. (2020). The Predictable Chaosof Slow Earthquakes. Science Advances , Accepted .Guglielmi, Y., Cappa, F., Avouac, J.-P., Henry, P., & Elsworth, D. (2015, jun). Seis-micity triggered by fluid injection-induced aseismic slip.
Science , (6240),1224–1226. Retrieved from doi: 10.1126/science.aab0476Hindmarsh, A. C. (1983). ODEPACK, A systematized collection of ODE solvers. In –27– . S. Stepleman & E. Al. (Eds.), Imacs transactions on scientific computation (Vol. 1, pp. 55–64). Amstredam.Hoffmann, H., & Payton, D. W. (2014). Suppressing cascades in a self-organized-critical model with non-contiguous spread of failures.
Chaos, Solitonsand Fractals , , 87–93. Retrieved from http://dx.doi.org/10.1016/j.chaos.2014.06.011 doi: 10.1016/j.chaos.2014.06.011Hofmann, H., Zimmermann, G., Farkas, M., Huenges, E., Zang, A., Leonhardt, M.,. . . Kim, K. Y. (2019). First field application of cyclic soft stimulation atthe Pohang Enhanced Geothermal System site in Korea. Geophysical JournalInternational , (2), 926–949. doi: 10.1093/gji/ggz058Huang, J., Narkounskaia, G., & Turcotte, D. L. (1992). A cellularautomata, slid-erblock model for earthquakes II. Demonstration of selforganized criticalityfor a 2D system. Geophysical Journal International , (2), 259–269. doi:10.1111/j.1365-246X.1992.tb00575.xHuang, J., & Turcotte, D. L. (1992). Chaotic seismic faulting with a mass-springmodel and velocity-weakening friction. Pure and Applied Geophysics PA-GEOPH , (4), 569–589. Retrieved from http://link.springer.com/10.1007/BF00876339 doi: 10.1007/BF00876339Ito, K., & Matsuzaki, M. (1990). Earthquakes as self-organized critical phenom-ena. Journal of Geophysical Research , (B5), 6853. Retrieved from http://doi.wiley.com/10.1029/JB095iB05p06853 doi: 10.1029/JB095iB05p06853Jensen, H. J. (1998). Self-Organized Criticality . Cambridge University Press.Retrieved from doi: 10.1017/CBO9780511622717Kanamori, H., & Brodsky, E. E. (2004). The physics of earthquakes.
Reports on Progress in Physics , (8), 1429–1496. Retrieved from http://stacks.iop.org/0034-4885/67/i=8/a=R03?key=crossref.0eb46da79cd6938ce542994e8554673e doi: 10.1088/0034-4885/67/8/R03Kenigsberg, A. R., Rivi`ere, J., Marone, C., & Saffer, D. M. (2020, mar). Evolutionof Elastic and Mechanical Properties During Fault Shear: The Roles of ClayContent, Fabric Development, and Porosity. Journal of Geophysical Research:Solid Earth , (3), 1–16. Retrieved from https://onlinelibrary.wiley.com/doi/abs/10.1029/2019JB018612 doi: 10.1029/2019JB018612Keranen, K. M., Savage, H. M., Abers, G. A., & Cochran, E. S. (2013). Potentiallyinduced earthquakes in Oklahoma, USA: Links between wastewater injectionand the 2011 Mw 5.7 earthquake sequence. Geology , (6), 699–702. doi:10.1130/G34045.1Khalil, H. k. (2015). Non-linear control: Global edition . Harlow: Pearson.Kluyver, T., Ragan-Kelley, B., P´erez, F., Granger, B., Bussonnier, M., Frederic,J., . . . Willing, C. (2016). Jupyter Notebooks – a publishing format for re-producible computational workflows. In F. Loizides & B. Schmidt (Eds.),
Positioning and power in academic publishing: Players, agents and agendas (pp. 87–90).Kwiatek, G., Saarno, T., Ader, T., Bluemle, F., Bohnhoff, M., Chendorain, M., . . .Wollin, C. (2019, may). Controlling fluid-induced seismicity during a 6.1-km-deep geothermal stimulation in Finland.
Science Advances , (5), eaav7224.Retrieved from http://advances.sciencemag.org/lookup/doi/10.1126/sciadv.aav7224 doi: 10.1126/sciadv.aav7224Lomnitz-Adler, J. (1993, oct). Automaton models of seismic fracture: Con-straints imposed by the magnitude-frequency relation. Journal of Geo-physical Research: Solid Earth , (B10), 17745–17756. Retrieved from http://doi.wiley.com/10.1029/93JB01390 doi: 10.1029/93JB01390Lyapunov, A. M. (1892). The general problem of the stability of motion (DoctoralThesis (in Russian)). University of Kharkov.Main, I. (1999, apr). Earthquake prediction: Concluding Remarks.
Nature (April), –28– –7. Retrieved from doi: 10.1038/nature28133Marone, C. (1998). The effect of loading rate on static friction and the rate of faulthealing during the earthquake cycle.
Nature , (6662), 69–72. doi: 10.1038/34157McClung, D., & Schaerer, P. (1993). The avalanche handbook . Seattle: The moun-tainers.McGarr, A., Simpson, D., & Seeber, L. (2002). Case Histories of Induced and Trig-gered Seismicity.
International Handbook of Earthquake and Engineering Seis-mology , (A), 0–12.Narkounskaia, G., Huang, J., & Turcotte, D. L. (1992). Chaotic and self-organizedcritical behavior of a generalized slider-block model. Journal of StatisticalPhysics , (5-6), 1151–1183. doi: 10.1007/BF01049013No¨el, P. A., Brummitt, C. D., & D’Souza, R. M. (2013). Controlling self-organizingdynamics on networks using models that self-organize. Physical Review Letters , (7), 19–23. doi: 10.1103/PhysRevLett.111.078701Parisio, F., Vilarrasa, V., Wang, W., Kolditz, O., & Nagel, T. (2019, dec). The risksof long-term re-injection in supercritical geothermal systems. Nature Commu-nications , (1), 4391. Retrieved from doi:10.1038/s41467-019-12146-0Petersen, M., Mueller, C., Moschetti, M., Hoover, S., Rubinstein, J. L., Llenos, A.,. . . Anderson, J. (2015). Incorporating Induced Seismicity in the 2014 UnitedStates National Seismic Hazard ModelResults of 2014 Workshop and Sensitiv-ity Studies (Tech. Rep.). U.S. Geological Survey.Petzold, L. (1983). Automatic Selection of Methods for Solving Stiff and NonstiffSystems of Ordinary Differential Equations.
SIAM Journal on Scientific andStatistical Computing , (1), 136–148. doi: 10.1137/0904010 Python Control Systems Library. (2020). Retrieved from https://python-control.readthedocs.io/en/0.8.3/
Raleigh, C. B., Healy, J. H., & Bredehoeft, J. D. (1976). An experiment in earth-quake control at Rangely, Colorado.
Science (New York, N.Y.) , (4233),1230–7. Retrieved from doi: 10.1126/science.191.4233.1230Rattez, H., Stefanou, I., & Sulem, J. (2018, jun). The importance of Thermo-Hydro-Mechanical couplings and microstructure to strain localization in 3Dcontinua with application to seismic faults. Part I: Theory and linear sta-bility analysis. Journal of the Mechanics and Physics of Solids , , 54–76. Retrieved from http://linkinghub.elsevier.com/retrieve/pii/S0022509617309626https://linkinghub.elsevier.com/retrieve/pii/S0022509617309626 doi: 10.1016/j.jmps.2018.03.004Rattez, H., Stefanou, I., Sulem, J., Veveakis, E., & Poulet, T. (2018, jun). Theimportance of Thermo-Hydro-Mechanical couplings and microstructureto strain localization in 3D continua with application to seismic faults.Part II: Numerical implementation and post-bifurcation analysis. Jour-nal of the Mechanics and Physics of Solids , , 1–29. Retrieved from http://linkinghub.elsevier.com/retrieve/pii/S0022509617309638 doi: 10.1016/j.jmps.2018.03.003Rattez, H., & Veveakis, M. (2020). Weak phases production and heat generationcontrol fault friction during seismic slip. Nature Communications , (1), 1–8.Retrieved from http://dx.doi.org/10.1038/s41467-019-14252-5 doi: 10.1038/s41467-019-14252-5Reid, H. F. (1910). The Mechanics of the Earthquake, The California Earthquake ofApril 18, 1906. In Report of the state investigation commission (Vol. 2). Wash-ington: Carnegie Institution of Washington. –29– ice, J. R. (1993). Spatio-temporal complexity of slip on a fault.
Journal of Geo-physical Research , (B6), 9885. Retrieved from http://doi.wiley.com/10.1029/93JB00191 doi: 10.1029/93JB00191Rubinstein, J. L., & Mahani, A. B. (2015). Myths and Facts on Wastewater In-jection, Hydraulic Fracturing, Enhanced Oil Recovery, and Induced Seis-micity. Seismological Research Letters , (4), 1060–1067. Retrieved from http://srl.geoscienceworld.org/lookup/doi/10.1785/0220150067 doi:10.1785/0220150067Ruina, A. (1983). Slip instability and state variable friction laws. Journal of Geo-physical Research , (B12), 10359. doi: 10.1029/JB088iB12p10359Rundle, J. B., & Brown, S. R. (1991). Origin of rate dependence in frictional sliding. Journal of Statistical Physics , (1-2), 403–412. doi: 10.1007/BF01329869Schmittbuhl, J., Vilotte, J.-P., & Roux, S. (1996, jun). Velocity weakening friction:A renormalization approach. Journal of Geophysical Research: Solid Earth , (B6), 13911–13917. Retrieved from http://doi.wiley.com/10.1029/96JB00653 doi: 10.1029/96JB00653Scholz, C. H. (2019). The Mechanics of Earthquakes and Faulting (Third ed.).Cambridge University Press. Retrieved from doi: 10.1017/9781316681473Segall, P., & Lu, S. (2015, jul). Injection-induced seismicity: Poroelastic and earth-quake nucleation effects.
Journal of Geophysical Research: Solid Earth , (7),5082–5103. Retrieved from http://doi.wiley.com/10.1002/2015JB012060 doi: 10.1002/2015JB012060Sornette, D. (1999, apr). Earthquake debate - Didier Sornette. Nature (April). Re-trieved from doi: 10.1038/nature28132Stefanou, I. (2019, aug). Controlling Anthropogenic and Natural Seismicity:Insights From Active Stabilization of the SpringSlider Model.
Journal ofGeophysical Research: Solid Earth , (8), 8786–8802. Retrieved from https://onlinelibrary.wiley.com/doi/abs/10.1029/2019JB017847 doi:10.1029/2019JB017847Stefanou, I., & Alevizos, S. (2016). Fundamentals of bifurcation theory andstability analysis. In J. Sulem, I. Stefanou, E. Papamichos, & E. Ve-veakis (Eds.), Modelling of instabilities and bifurcation in geomechanics,alert geomaterials doctoral school 2016.
Aussois, France. Retrieved from
Strogatz, S.-H. (1994).
Non linear dynamics and chaos . Perseus Books.Terzaghi, K. (1925).
Erdbaumechanik auf bodenphysikalischer grundlage . Wien: F.Deuticke.Turcotte, D. L. (1999). Self-organized criticality.
Reports on Progress in Physics , (10), 1377–1429. doi: 10.1088/0034-4885/62/10/201van Rossum, G. (1995). Python tutorial. CWI Report CS-R9526 , CS-R9526 , 1–65.Retrieved from http://oai.cwi.nl/oai/asset/5007/05007D.pdf
Vardulakis, A. I. (1991).
Linear Multivariable Control: Algebraic Analysis and Syn-thesis Control . Chichester, New York, Brisbane, Toronto, Singapore: John Wi-ley & Sons, Inc.Vardulakis, A. I. (2012).
Introduction to the mathematical theory of the theory of sig-nals, systems and control . Tziola.Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau,D., . . . van Mulbregt, P. (2020, mar). SciPy 1.0: fundamental algorithmsfor scientific computing in Python.
Nature Methods , (3), 261–272. Re-trieved from doi:10.1038/s41592-019-0686-2 –30– atkins, N. W., Pruessner, G., Chapman, S. C., Crosby, N. B., & Jensen, H. J.(2016, jan). 25 Years of Self-organized Criticality: Concepts and Contro-versies. Space Science Reviews , (1-4), 3–44. Retrieved from http://dx.doi.org/10.1007/s11214-015-0155-xhttp://link.springer.com/10.1007/s11214-015-0155-x doi: 10.1007/s11214-015-0155-xWolfram, S. (1983, jul). Statistical mechanics of cellular automata. Reviews of Mod-ern Physics , (3), 601–644. Retrieved from https://link.aps.org/doi/10.1103/RevModPhys.55.601 doi: 10.1103/RevModPhys.55.601doi: 10.1103/RevModPhys.55.601