J P Dalmont

S Maugeais

Piano strings with reduced inharmonicity

Even on modern straight pianos, the inharmonicity of the lower strings is rather large especially for the first octave. Consequently, the timber of these strings can sometimes sound awful and chords on the first octave be highly dissonant. The idea of the present study is to show how this defect can be rectified using an inhomogeneous winding on the whole string in order to minimize inharmonicity. The problem is solved using an optimisation procedure considering a non uniform linear density. Results show that the inharmonicity of the first partials could be highly reduced by a non uniform winding limited to a quarter of the string.

Introduction

So-called harmonic strings are largely used in music because a uniform string without stiffness, and stretched between two fixed points, naturally have harmonic eigenfrequencies. When considering the string's stiffness the eigenfrequencies are no longer harmonic. The consequence on instruments like harpsichord and early pianos is limited. However, the development of piano making during the 19th century saw a tendency of increasing string tension by a factor 4, and the mass in the same proportion (cf. [START_REF] Chaigne | The making of pianos : a historical view[END_REF]). A consequence is that the inharmonicity cannot be considered negligible anymore, which reflects on the tuning of the instrument and the timbre. For the lower strings a solution has been found in order to increase the mass of the string without increasing too much its bending stiffness: the wound strings. Nevertheless, the inharmonicity remains rather large especially for the first octave of medium grand piano, and even worse on an upright piano: according to Young the inharmonicity of the bass strings on a medium piano is twice that of a grand piano, and that of a straight piano twice that of a medium ( [START_REF] Young | Inharmonicity of piano bass strings[END_REF], [START_REF] Giordano | Evolution of music wire and its impact on the development of the piano[END_REF]). Our study focuses on straight piano because we consider that designing strings with reduced inharmonicity would improve quite a lot the musical quality of these instruments.

The idea of the present study is to show how this defect could be rectified by using an inhomogeneous 42 winding on the whole string, in order to minimise in-43 harmonicity. The string is thus considered to be inho-44 mogeneous that is with a non uniform linear density 45 ([8]). From a theoretical point of view, the problem 46 translates into finding an "optimal" non uniform lin-47 ear density for a stretched string with uniform stiff-48 ness. Here, the optimality condition amounts to being 49 as harmonic as possible. This problem is solved using 50 an optimisation procedure, initialised with the char-51 acteristics of a real string. The diameter of the opti-52 mised string is allowed to vary between the diameter 53 of the core (supporting the winding) and about twice 54 the diameter of the reference string. An area with 55 uniform winding will be kept to reduce the amount 56 of work during the manufacturing of the true string. 57 Moreover, it is proposed to limit the non uniform 58 winding to one side of the string. The model chosen for the string is linear and only 61 involves tension, bending stiffness and mass per unit 62 length. According to Chabassier, this model is suffi-63 cient for low frequencies (cf. [START_REF] Chabassier | Modeling and numerical simulation of a piano[END_REF], remark I.1.2, there 64 is no need to add a shear term), small amplitudes. 65 Moreover, even if this not completely true (see [START_REF] Conklin | Design and tone in the mechanoacoustic piano. part iii. piano strings and scale design[END_REF]), it 66 is considered that the increase of stiffness due to the 67 wrapping can be neglected. So, the stiffness is that of 68 the core (cf. [1] I.1.5) and is therefore constant along 69 the string. Finally, the model is taken without any 70 losses as only the eigenfrequencies are of interest. The displacement equation is then given in the Fourier domain by (cf. [START_REF] Ducasse | On waveguide modeling of stiff piano strings[END_REF])

-µω 2 y = T ∂ 2 y ∂x 2 -EI ∂ 4 y ∂x 4 (1)
where µ is the mass per unit length (function of x), T is the tension, E is the string core's Young modulus, I = Ar 2 /4 with A the core section, r the core radius and ω is the pulsation.

The string being simply supported at both ends, the boundary conditions are given by

y| x=0 = y| x=L = ∂ 2 y ∂x 2 x=0 = ∂ 2 y ∂x 2 x=L = 0.

Solution with constant µ

When µ(x) = µ 0 is constant, it is possible to find an explicit solution of (1). The eigenfrequencies of the oscillator are given by (cf. [START_REF] Valette | Mechanics of the vibrating string[END_REF] §3.4)

f n = nf 0 1 + Bn 2 with f 0 = 1 2L T µ 0 and B = π 2 EI T L 2 .
The inharmonicity factor comes from

√ 1 + Bn 2 .
This equation shows the influence of string's stiffness on inharmonicity. A shorter string with small tension and high stiffness has a higher B, and therefore inharmonic eigenfrequencies. For a typical piano string, the inharmonicity is minimum for the second octave and in the range of 10 -4 . It increases with the frequency in upper octaves but also for the first octave ( [START_REF] Rigaud | A parametric model and estimation techniques for the inharmonicity and tuning of the piano[END_REF]). For a grand piano for the first note A0, B is less than 10 -4 which leads to an inharmonicity of 20 cents for the 16th harmonic but for a straight piano B can reach 10 -3 which leads to an inharmonicity of 200 cents for the 16th harmonic (see [START_REF] Conklin | Design and tone in the mechanoacoustic piano. part iii. piano strings and scale design[END_REF]). With such values of B, the sound of the lower string is awful and chords on the first octave cause a lot of beatings inducing a high roughness. Therefore, we consider that the reduction of the harmonicity factor would have a beneficial influence on the sound of the medium and straight piano, especially for the first two octaves.

Optimisation

The goal is to find a density function that corresponds to a harmonic string, i.e. a function µ that is a minimum for the "inharmonicity" function C (cf. [START_REF] Gaudet | The evolution of harmonic indian musical drums: A mathematical perspective[END_REF]) defined by

C : µ → nmax n=2 ω n (µ) nω 1 (µ) -1 2
where the ω n s are the eigenfrequencies of the equation This choice of inharmonicity function is quite natural for an optimisation problem as it is a quadratic function of the higher frequencies.

In general, it is not possible to work directly with 106 the function µ, and a discretisation of the space is 107 needed so that an approximation in finite dimension 108 can be used. When mass is non uniform, solutions of equation ( 1) have to be approximated by a numerical method. The problem with non-constant µ is thus solved with a classical FEM in space (similar to the more complex setting of [START_REF] Chabassier | Modeling and numerical simulation of a piano[END_REF], §II.1) using Hermite's polynomials (cf. [START_REF] Ern | Aide-memoire of finite elements[END_REF], 1.7) and a uniform discretisation of [0, L] for fixed N and 0 < i < N + 1 with h = L N +1 , x i = hi. The projection of the operators "multiplication by µ", T ∂ 2 ∂x 2 and EI ∂ 4 ∂x 4 then define three 2N × 2N matrices M(µ), T and E so that (1) can be approximated by

ω 2 M(µ)U = (T + E)U (2)
with

U = t (y(x 1 ), y (x 1 ), • • • , y(x N ), y (x N ))
, 111 where the prime denotes the spatial derivative.

112

More precisely, the matrices T and E can be com-113 puted using Hermite's polynomials and leads to the 114 formulas 115

T = T h 2            . . . -6 5 h 10 -h 10 -h 2 30 -6 5 -h 10 12 5 0 h 10 -h 2 30 0 4h 2 15 . . .            , E = EI h 4           . . . -12 6h -6h 2h 2 -12 -6h 24 0 6h 2h 2 0 8h 2 . . .           116
and as well the matrix M(µ) is computed for func-117 tions µ constants on each interval ]x i , x i+1 [ by

118 M =            . . . µ i-1/2 9 70 -µ i-1/2 13h 420 µ i-1/2 13h 420 -µ i-1/2 h 2 140 µ i-1/2 9 70 µ i-1/2
13h 420

(µ i-1/2 + µ i+1/2 ) 13 35 (-µ i-1/2 + µ i+1/2 ) 11h 210 -µ i-1/2 13h 420 -µ i-1/2 h 2 140 (-µ i-1/2 + µ i+1/2 ) 11h 210 (µ i-1/2 + µ i+1/2 ) h 2 105 . . .            119
where µ i+1/2 denotes the value of µ on ]x i , x i+1 [.

120

Equation ( 2) is a generalised eigenvalue problem 121 that can be solved using generalised Schur decompo-122 sition. In the present study, the gradient algorithm is used 125 to find a minimum of the "inharmonicity" function C. 126 It works well for low stiffness and/or small numbers 127 n max . When converging, it gives a solution having 128 arbitrarily low inharmonicity for the first n max har-129 monics.

130

It uses the derivative of functions ω n with respect to µ for n ∈ {1, • • • , n max }. The value of this derivative is computed using perturbation theory and is approximated by the formula

grad µ ω n = U n .U n t U n M(µ)U n
where U n is an eigenvector of equation ( 2) associated 131 to the eigenvalue -ω 2 n , and U n .U n is the Hadamard 132 (entrywise) product of U n with itself.

The stopping condition is dictated by the inharmonicity of each partial, as it should be at least as good as that of the uniform string and at most ≤ ε, for a fixed constant ε which, in practice, is taken to be 10 -3 .

Examples and results with different strategies

Gradient algorithm looks for solutions in a vector space, but most of the elements of this space are physically irrelevant. It is therefore necessary to reduce the search space and add conditions to find useful solutions. In particular, they must verify at least that the density remains bigger than that of the core, as it would otherwise weaken the string. On the other hand it can be interesting to limit the inhomogeneous part in order to make the manufacturing easier. These considerations lead to different strategies which are described below.

In the following the chosen nominal string is a C1 straight piano string corresponding to a frequency of 32.7Hz : L = 1.035m, µ 0 = 180g/m, T = 825N and EI = 0.028N m 2 , which leads to an inharmonicity coefficient B = 3.13 10 -4 . Here the string is considered to be strictly uniform on all its length and it is considered that the stiffness is that of the core only. In practice it is likely that the inharmonicity coefficient might be significantly higher ( [START_REF] Conklin | Design and tone in the mechanoacoustic piano. part iii. piano strings and scale design[END_REF]).

Minimum of constraints

This corresponds to the case were the only constraint is that the density is at least equal to that of the wire. Results depend on the number of harmonics which are taken into account in the optimisation. On be noticed that the fluctuations are more important 187 near the ends off the string which suggests that fluc-188 tuations in the middle of the string could be avoided. 189 It can be seen on figure 2 that the diameter at the 190 ends is much larger than that of the uniform string. 191 In practice this will be probably difficult to manage. 192

Non negative density fluctuations 193

The question is now whether it would be possible to 194 only add some masses to a uniform string in order 195 to apply local corrections. On the numerical point of 196 view this is obtained by replacing the gradient by its 197 positive part. The convergence is much slower than 198 for the previous case (4730 steps are needed for a max-199 imum inharmonicity of ε = 10 -3 against 332). How-200 ever, the comparison of figure 2 and figure 3 shows 201 that it should be possible to avoid valleys without in-202 creasing too much the amplitude of the hills.

One sided density fluctuation

Owing to the symmetrical nature of the problem, for all the previous examples the algorithm converges to a symmetrical solution. However there is a priori no reason to keep a symmetrical string. Moreover, to avoid the hammer hitting the string on a non uniform region it would be good to limit the non uniform region to only one end of the string. The first mass is about 5 g and is centered at 3 cm 240 from the end. The third one is 0.4 g and is centered 241 at 24 cm from the end. In practice, the second mass 242 is probably useless. So, it is surprising to realise that 243 for n max = 10 a single mass is probably sufficient to 244 significantly improve the harmonicity of a string. It is noticeable that the present result is not far 246 from what can be found in [START_REF] Sanderson | Method for making wound strings for musical instruments characterized by reduced inharmonicity[END_REF] in which a local over 247 winding is used to compensate for the inharmonicity 248 induced by the bare ends. The present study shows that it might be possible 251 to build piano strings with an optimised non uniform 252 density leading to a reduced inharmonicity of the first 253 partials. Moreover the non uniform part of the string 254 could be limited to a short portion of the string near 255 the end (25 cm in the given example). Moreover, it 256 appears that a single mass a few centimeters from 257 the end might significantly improve the harmonicity. 258 Now, many questions arise. In practice, how to take 259 into account the effective inharmonicity of the string 260 including the influence of the sound board and that of 261 the winding? On the perceptive point of view, what 262 is the minimum number of harmonics to consider in 263 order to obtain a significant improvement? To answer 264 these questions, the next step is to find a process to 265 optimise the harmonicity on an actual string. Obvi-266 ously, the final design will be the result of exchanges 267 between piano tuners, pianists and engineers. Finally 268 we are deeply convinced that such harmonic strings 269 will make it possible to highly increase the musical 270 quality of the first two octaves of straight pianos but

71 2 . 1

 21 Mathematical model 72 Let us consider the Euler-Bernoulli model without 73 losses for a stiff string of length L.

  74

Figure 1 :

 1 Figure 1: Sketch of the string and notations

( 1 )

 1 sorted by increasing magnitude, n max being the number of harmonics considered in the optimisation.

109 3 . 1

 31 Numerical implementation 110

123 3 . 2

 32 Optimisation algorithm 124

figure 2 Figure 2 :

 22 figure 2 the density for different n max = 10, 15, 18, is shown. The convergence is considered to be obtained when the maximum of harmonicity is < 10 -3 . The number of steps needed to obtain the convergence increases with the number of harmonics: for n max = 10, 15, 18 the number of steps are respectively n step = 332, 848, 1188. It appears that the density fluctuations show a number of valleys equal to the number of harmonics involved in the optimisation. Physically it can be interpreted as a way to slower the waves up to the maximum frequency considered. Another important observation is that density fluctuations are allthe more important as the number of harmonics considered is high: the amplitude varies roughly as the square of the harmonic number. This is an important limitation, because this limits the number of harmonics on which inharmonicity can be minimised. On the given example, for n max = 18 the density is localy multiplied by about 2 which means that the diameter of the string is locally increased by 40% . It can also

Figure 3 :

 3 Figure 3: Result of convergence with non negative fluctuation of density and n max = 10. Top: density profile; bottom: harmonicity as a function of harmonic number.

4. 3 Figure 4 :

 34 Figure 4: Result of convergence with no fluctuation on one half of the string in the middle and n max = 10. Top: density profile; bottom: harmonicity as a function of harmonic number.

  Figure 5 shows that similar results are obtained when the fluctuation are concentrated on one side. As one could expect, this tends to increase the amplitude of the fluctuations. Results are given on figure 5 for n max = 10 and a = L/4. It is noticeable that the result tends to a point mass near the end combined with two periods of a "damped sinusoidal" variation of density.

4. 5

 5 One sided non negative density fluctuation Now, the next step is to consider only non negative density fluctuations on a small part of the string. The results given on figure 5 converge to what can be interpreted as small masses. The second mass being rather negligible it seems that two masses are sufficient to correct the inharmonicity of the first ten harmonics. 239

Figure 5 :

 5 Figure 5: Result of convergence with and without positive density fluctuations on one quarter of the string and n max = 10. Top: density profile; bottom: harmonicity as a function of harmonic number.


Ioannis Stefanou
email: ioannis.stefanou@enpc.fr

Controlling Anthropogenic and Natural Seismicity: Insights From Active Stabilization of the Spring-Slider Model

Keywords: induced seismicity, stability, control theory, earthquake control

Introduction

Given the current intense human activity for conventional and unconventional energy production (e.g. oil and gas), renewable energies (e.g. geothermal) and potential environmental friendly methods related to climate change (e.g. CO 2 sequestration), avoiding anthropogenic seismicity is a challenging topic. Even though evidence and proof of the origin of reported seismicity and its relation with human activities will always have a degree of uncertainty, it is nowadays generally accepted that humans can induce or trigger earthquakes. Anthropogenic earthquakes are usually of small to moderate magnitude, i.e. less than M w = 4 (moment magnitude). However, they usually exceed the acceptable limits set by the authorities and can cause damage. Moreover, there are several cases where induced earthquakes had important magnitudes (no distinction is made here between triggered and induced earthquakes and both terms are used interchangeably for simplicity). Some examples of induced earthquakes, among several others, are the M w = 5.3 at Trinidad, Colorado earthquake (EQ) in US due to wastewater injection [START_REF] Rubinstein | Myths and Facts on Wastewater Injection, Hydraulic Fracturing, Enhanced Oil Recovery, and Induced Seismicity[END_REF], the M w = 5.7 at Prague, Oklahoma EQ in US prob-ably due to wastewater injection [START_REF] Keranen | Potentially induced earthquakes in Oklahoma, USA: Links between wastewater injection and the 2011 Mw 5.7 earthquake sequence[END_REF][START_REF] Mcgarr | Maximum magnitude earthquakes induced by fluid injection[END_REF] and the M w = 3.6 at Basel geothermal project in Switzerland [START_REF] Cornet | Seismic and aseismic motions generated by fluid injections[END_REF][START_REF] Deichmann | Earthquakes Induced by the Stimulation of an Enhanced Geothermal System below Basel (Switzerland)[END_REF]. In fact the number and the importance of induced seismicity events were such that the United States Geological Survey (USGS) incorporated them in the 2014 United States National Seismic Hazard Model [START_REF] Petersen | Incorporating Induced Seismicity in the 2014 United States National Seismic Hazard Model-Results of 2014[END_REF].

Here we address the general question of earthquake control: Is it possible to control and avoid anthropogenic and natural earthquakes and how? We address this question using a mathematically rigorous framework and keeping complexity to a minimum degree. In particular, we investigate the conditions under which a simplified, seismogenic fault system can be stabilized by automatic fluid pressure adjustment/control. The concept that pore pressure increase due to fluid injection can stimulate fault slip is well established, nowadays [START_REF] Frohlich | Two-year survey comparing earthquake activity and injectionwell locations in the Barnett Shale, Texas[END_REF][START_REF] Healy | The Denver Earthquakes[END_REF][START_REF] Hubbert | Role of Fluid Pressure in Mechanics of Overthrust Faulting[END_REF][START_REF] Raleigh | An experiment in earthquake control at Rangely, Colorado[END_REF]. However, the question of EQ control and EQ risk mitigation is still open and challenging.

Evidence of EQ control from field experiments is very limited. We refer to the seminal field experiment of [START_REF] Raleigh | An experiment in earthquake control at Rangely, Colorado[END_REF] in the 70's in Rangely, Colorado, US, where earthquakes could be turned off and on by varying the pore pressure. Another example of EQ control is in Dale, New York, US [START_REF] Fletcher | Earthquakes related to hydraulic mining and natural seismic activity in western New York State[END_REF], where earthquakes of magnitude -1 to 1.4 formed a cluster about 650m across near the bottom of a 426m injection well. The earthquake activity was arrested when the top hole pressure dropped below 5 MPa. More recent field experiments involve the well monitored tests by [START_REF] Cappa | Stabilization of fault slip by fluid injection in the laboratory and in situ[END_REF]; [START_REF] Guglielmi | Seismicity triggered by fluid injection-induced aseismic slip[END_REF] performed at 252m depth within the LSBB underground laboratory in France. In these field experiments aseismic slip was systematically preceding seismic slip, giving concrete evidence that slip can be also aseismic.

Here we follow a theoretical approach in order to get useful insight of the controllability of a fault system and the possibility of injecting fluids in a way that guarantees aseismic slip (definitions of the terms seismic and aseismic slip, as used here, are given in Section 2.2). Our analysis is based on the classical spring-slider model. We adopt a general frictional law that depends on slip and rate of slip. Additionally to rock elasticity (spring) we consider also radiational damping. This is performed by accounting for the rock viscosity through a damper in Kelvin-Voigt configuration (see also [START_REF] Wang | Multistable slip of a one-degree-of-freedom spring-slider model in the presence of thermal-pressurized slip-weakening friction and viscosity[END_REF]).

We consider that the rock is saturated and that our system is well oriented for slip in the ambient stress-field. Both physics and geometry of the system are kept as simple as possible. The role of heterogeneities, of pore fluid diffusivity and special hydrological conditions are not considered in the present work. Consequently, we don't focus on EQ rupture and propagation in details, but only on average (over the fault's length) using the spring-slider model. The above mentioned aspects, as well as the observability of the real system and other techno-economical aspects of EQ control, exceed the scope of the present article and they are explored in the frame of the ongoing ERC project "Controlling earth-Quakes -CoQuake" (see http://coquake.com).

We give particular emphasis on the study of the stability of the system which is constantly driven by the far-field tectonic velocity. The term stability is used here only in the sense of Lyapunov stability (i.e. the system remains close to its equilibrium state under small perturbations from it; for a rigorous mathematical definition of Lyapunov stability we refer to [START_REF] Lyapunov | The general problem of the stability of motion[END_REF][START_REF] Stefanou | Fundamentals of bifurcation theory and stability analysis[END_REF]). Due to the far field tectonic velocity and the general rheology considered for friction, the system is non-autonomous (i.e. it depends explicitly on time) and the classical methods of Lyapunov stability cannot be directly applied. For this purpose, we use a double-time scale asymptotic analysis that a) eliminates the secular terms (i.e. growing, unbounded terms in time), b) provides the steady-state of the system describing its slow, creep-like motion and c) allows to derive the (in-)stability conditions, i.e. determine when slow, creep-like slip is possible and when an earthquake takes place. Based on these results we prove mathematically that the system is stabilizable by fluid pressure control. This is a major result. An opposite conclusion would mean that EQ control is impossible, implying inevitable risks for on-going, large-scale industrial applications.

Following these mathematical developments and using the classical mathematical theory of control we design a stabilizing controller (compensator ) and we investigate the conditions for which the controller can stabilize the system, even when complete knowledge on the evolution of the frictional properties of the system is not available. The ability of the controller to stabilize the system and avoid unwanted seismicity, is then illustrated through two scenarios of fluid injection at 5km depth. The first one refers to an injection project, where an injection under constant pressure rate is planned. It is shown that the controller automatically stabilizes the system and avoids the anticipated earthquake event of M w = 3.8. The second scenario concerns a fault system with higher rup-ture area, able to give earthquakes of M w = 5.8. In this second scenario, we drive the system from its initial, high energy, unstable state to a new, stable one. This is done by automatically adjusting the fluid pressure such as to assure a constant slow slip rate. In this way no seismicity is observed and the earthquake is avoided. Finally, the robustness of the designed controller to successfully mitigate seismicity is challenged by perturbing the friction properties of the system.

The paper is organized as follows. In Section 2 we present the basic modeling assumptions of the spring-slider analog and we study it using double-scale asymptotics in time. The steady-state motion is approximated by a power series of the orders of the slow time scale, for which we calculate explicitly the dominant one. The necessary and sufficient conditions for stable and unstable steady-state motion are determined. This asymptotic analysis renders the dominant order of the system time invariant and allows us to study it further in Section 3 using the classical mathematical theory of control. The linearization of the system is investigated in paragraph 3.2 and its stabilizability in paragraph 3.3. A stabilizing controller is designed in paragraph 3.4. Asymptotic tracking and robustness of the controller are discussed next. Finally, in Section 4 we present two numerical examples, i.e. the aforementioned scenarios of EQ control.

2 Steady-state and stability conditions for the spring-slider analog Consider the spring-slider model of Figure 1a. This is the classical paradigm and starting point for studying the dynamic instability of earthquake nucleation in a mathematically simplified manner (see [START_REF] Burridge | Model and theoretical seismicity[END_REF]; [START_REF] Reid | The Mechanics of the Earthquake, The California Earthquake of April 18[END_REF]). In this model, the block represents the mobilized rock mass, m, during an earthquake event. With the spring, k, we model the elastic deformation of the surrounding to the fault rock. This allows to account for the progressive elastic energy build-up due to the far field tectonic movement, δ ∞ . The far field tectonic movement is of the order of some cm's per year, contrary to the seismic slip that can rise up to one meter per second. The coefficient of the elastic spring is proportional to the effective shear elastic modulus of the surrounding to the fault rocks, G, and inversely proportional to the fault length, L, i.e. k ∝ G L .

This scaling is retrieved from elasticity theory [START_REF] Palmer | The Growth of Slip Surfaces in the Progressive Failure of Over-Consolidated Clay[END_REF]. With effective shear modulus we mean the apparent/averaged shear modulus over a region that extends at a distance L in the direction and perpendicular to the fault (see Figure 1b). This region includes the damaged area that extends from some meters to some kilometers around the fault [START_REF] Okubo | Dynamics, radiation and overall energy budget of earth[END_REF]. Additionally, we consider the effective viscosity of the surrounding rock, C. This viscosity is represented through an equivalent dashpot of coefficient η, connected in a Kelvin-Voigt configuration as shown in Figure 1. With the term fault we mean the region of the rock system that is under pronounced localized shear deformation and opposes to the movement of the block, with a frictional force F . A fault is usually modeled as a mathematical plane due to its very small thickness, h, compared to the other characteristic lengths of the problem. Nevertheless, its thickness is linked with the softening response of the system, i.e. the reduction of friction in function of slip, rate of slip and other variables related to multiphysical couplings.

This region of extreme shearing is usually consisted of ultracataclastic materials and it has a complex structure [START_REF] Ben-Zion | Characterization of Fault Zones[END_REF][START_REF] Brodie | Structural terms including fault rock terms[END_REF] due to various physico-chemical phenomena that take place during pre-and coseismic slip (see [START_REF] Anthony | Influence of particle characteristics on granular friction[END_REF]; Rattez, Stefanou, and Sulem (2018); Rattez, Stefanou, Sulem, Veveakis, andPoulet (2018a, 2018b); [START_REF] Reches | Fault weakening and earthquake instability by powder lubrication[END_REF]; [START_REF] Scuderi | Frictional stability and earthquake triggering during fluid pressure stimulation of an experimental fault[END_REF]; [START_REF] Tinti | On the evolution of elastic properties during laboratory stick-slip experiments spanning the transition from slow slip to dynamic rupture[END_REF], among others). As a result, the apparent friction, F , does not depend only on the extent and the rate of slip, δ, v = δ, but also on the evolution of the microstructural network, the grain size, the presence and pressure of interstitial fluids, the temperature, time (state) and the reactivation of chemical reactions [START_REF] Brantut | Strain Localization and Slip Instability in a Strain-Rate Hardening, Chemically Weakening Material[END_REF]; [START_REF] Veveakis | Chemical reaction capping of thermal instabilities during shear of frictional faults[END_REF]; [START_REF] Veveakis | Failure in shear bands for granular materials: thermo-hydro-chemo-mechanical effects[END_REF], among others).

As far as it concerns the triggering of the dynamic instability studied here, the knowledge of the constitutive description of the apparent friction is central. Depending on the way that the fault's apparent friction evolves with shearing, an earthquake can be nu-cleated or arrested. The transition of the system from a potential unstable state to a stable one can be studied mathematically.

Let the apparent fault friction depend on δ and δ, i.e. F = F (δ, δ). The springslider analog is described by the following equation:

m δ = -F r (δ, δ) + k(δ ∞ -δ) + η(v ∞ -v) (1)
where δ ∞ is the far-field displacement due to the far-field velocity v ∞ . Notice that the above expression considers the dependence of friction on the degree of freedom of the system (slip) and its first derivative (rate of slip) and not on any internal, state variables (see rate and state friction laws). Nevertheless, it is often possible to eliminate these internal variables and express friction as a function of slip and rate of slip only. Later, in paragraph 2.3, the apparent fault friction will depend on the fluid pressure too, which will allow to control seismic slip.

It is worth emphasizing that the above equation is a non-autonomous, non-linear dynamic system, whose stability and steady state cannot be directly studied using the classical Lyapunov methods [START_REF] Brauer | The Qualitative Theory of Ordinary Differential Equations: An Introduction[END_REF]. Stability of this system has been studied in the literature by several researchers who considered it ad-hoc as autonomous by either neglecting the far field velocity or by applying a constant force. Nevertheless, this is a strong assumption that can lead, in general, to incorrect results regarding stability [START_REF] Brauer | The Qualitative Theory of Ordinary Differential Equations: An Introduction[END_REF]. Here we follow a different approach that allows us to asymptotically approximate the steady-state movement of the block and study its (Lyapunov) stability.

Due to the fact that the far-field tectonic movement is many orders of magnitude slower than seismic slip (v ∞ is a very small quantity) the steady-state solution of the springslider can be asymptotically approximated using the double-scale approach presented in the following paragraph.

Double-scale asymptotic analysis

Equation ( 1) can be written as follows:

ÿ = -a(y, ẏ) + c (εt -y) + d (ε -ẏ) , (2)
where ÿ and ẏ denote, respectively, the second and first order derivatives with respect to time, t, of the unknown function y = y(t), ε is a small parameter (ε 1) and a, c, d are of O(ε 0 ). The above equation is dimensionless y expresses normalized displacement and the exact expressions of coefficients are given in Section 4 (Eq.( 41)), where the theory is applied to specific scenarios of earthquake control.

We introduce a new variable τ = εt, which defines a long time scale because it is not negligible when τ is of order ε -1 or larger. In this way the system is expressed in terms of two time scales, a slow one, τ , and a fast one, t:

ÿ = -a(y, ẏ) + c (τ -y) + d (ε -ẏ) . (3)
This double scale approach allows to seek solutions of y, which are functions of both variables t and τ , treated as independent [START_REF] Bender | Advanced Mathematical Methods for Scientists and Engineers I[END_REF]. It is worth emphasizing that expressing y as a function of two variables is an artifice to remove secular effects. The exact solution y(t) is a function of t alone. We assume the following perturbation expansion for the unknown function y(t):

y(t) = y 0 (t, τ ) + εy 1 (t, τ ) + ε 2 y 2 (t, τ ) + O(ε 3 ). ( 4
)
Using the chain rule for differentiation, we obtain:

dy(t) dt = ∂y 0 ∂t + ε ∂y 0 ∂τ + ∂y 1 ∂t + ε 2 ∂y 1 ∂τ + ∂y 2 ∂t + O(ε 3 ) (5)
and

d 2 y(t) dt 2 = ∂ 2 y 0 ∂t 2 + ε 2 ∂ 2 y 0 ∂τ ∂t + ∂ 2 y 1 ∂t 2 + ε 2 ∂ 2 y 0 ∂τ 2 + 2 ∂ 2 y 1 ∂τ ∂t + ∂ 2 y 2 ∂t 2 + O(ε 3 ). ( 6
)
We assume that a(y, ẏ) can be expanded in power series in terms of ε:

a(y, ẏ) = a(y 0 , ∂y 0 ∂t ) + ε ∂a ∂y (y0, ∂y 0 ∂t ) y 1 + ∂a ∂ ẏ (y0, ∂y 0 ∂t ) ∂y 0 ∂τ + ∂y 1 ∂t + O(ε 3 ). ( 7
)
Inserting the above equations into Eq.( 2) and collecting powers of ε, we obtain the following cascade problems:

ε 0 : ∂ 2 y 0 ∂t 2 = -a(y 0 , ∂y 0 ∂t ) + c (τ -y 0 ) -d ∂y 0 ∂t ε 1 : ∂ 2 y 1 ∂t 2 = - ∂a ∂y (y0, ∂y 0 ∂t ) + c y 1 - ∂a ∂ ẏ (y0, ∂y 0 ∂t ) + d ∂y 1 ∂t - ∂a ∂ ẏ (y0, ∂y 0 ∂t ) + d ∂y 0 ∂τ -2 ∂ 2 y 0 ∂t∂τ + d ε 2 : . . . (8)

Slow dynamics and Lyapunov stability

We will first study the first of Eqs.( 8). Setting q(t) = ∂y0(t,τ ) ∂t and r(t) = y 0 (t, τ )

we obtain the following first order, non-linear system of ODE's:

q = -a(r, q) + c (τ -r) -d q ṙ = q. (9)
Notice that in this system, τ acts as a parameter, given that t and τ are independent variables. The above system can be written in vectorial form ż = f ( z), with z = (z 1 , z 2 ) = (q, r), is autonomous and its stability can be studied using the classical tools of Lyapunov stability theory (see [START_REF] Brauer | The Qualitative Theory of Ordinary Differential Equations: An Introduction[END_REF]). This system has an equilibrium point (fixed point) at z * = (q * , r * ) satisfying:

-a(r * , q * ) + c (τ -r * ) = 0 q * = 0. (10)
The above fixed point shows that the system has zero velocity in terms of the fast time variable t, i.e. q * = ∂y * 0 (t,τ ) ∂t = 0, leading to y 0 = y 0 (τ ), which is a function of the slow time variable, only. It should be emphasized that the velocity in terms of the slow time variable τ is not zero.

The eigenvalues, λ, of the Jacobian of f satisfy the characteristic polynomial:

λ 2 + λ ∂a ∂z 1 z * + d + ∂a ∂z 2 z * + c = 0. ( 11
)
The system is unstable when a positive eigenvalue exists and stable when all the eigenvalues are negative. Therefore, if for any τ , exists z * = (0, r * ), such that ∂a + d = 0 is of no particular interest for the physical problem at hand and it is not studied herein. However, we can show that in this particular case the system is stable.

The time evolution of the stable solution of the system is y 0 = y * 0 = y * 0 (τ ). We call here this state steady-state, given that

∂y * 0 ∂t = 0. Eq.(8).2 becomes: ∂ 2 y 1 ∂t 2 = -αy 1 -β ∂y 1 ∂t -β ∂y * 0 ∂τ + d, (12)
where α = ∂a ∂y (y * 0 ,0)

+c, β = ∂a ∂ ẏ (y * 0 ,0)
+d and ∂y * 0 ∂τ can be determined by differentiating Eq.(10.1) in terms of τ :

∂y * 0 ∂τ = c ∂a ∂y | (y * 0 ,0) + c = c α . ( 13
)
The general solution of ( 12) is:

y * 1 (t, τ ) = αd -βc α 2 + C 1 e 1 2 (-β+ √ ∆)t + C 2 e 1 2 (-β- √ ∆)t , (14)
where C 1 and C 2 are constants determined by the initial conditions of the ε 1 problem and ∆ = β 2 -4α.

The asymptotic approximation of the steady-state solution of the problem, y * , is therefore:

y * (t) = y * 0 (εt) + εy * 1 (t, εt) + O(ε 2 ). ( 15
)
The steady state velocity, v * (t), is approximated by the series (see Eq.( 5)):

v * (t) = ε c α + C 1 1 2 -β + √ ∆ e 1 2 (-β+ √ ∆)t + C 2 1 2 -β - √ ∆ e 1 2 (-β- √ ∆)t + O(ε 2 ). ( 16
) Notice that if β > 0, v * (t) approximates asymptotically the solution v * (t) = ε c α + O(ε 2
). Positive β is usually the case, due to the high viscosity of the surrounding to the fault rocks.

Conditions for steady-state slip of the spring-slider analog

Considering a monotonous motion ( δ ≥ 0) of the spring-slider system, the friction can be considered as a non-linear function of slip δ and can be expanded in power series as done in Eq.( 7). Under these assumptions, the multiscale asymptotic approach presented above can be used for the spring-slider system. More specifically, according to Eq.( 16), the slip-rate at steady-state of the spring-slider model is given by:

v ss = v ∞ 1 + 1 k ∂F ∂δ , (17)
where ∂F ∂δ is calculated at a given δ(t) = δ 0 and at zero slip-rate δ = 0. The above relation can give a useful global estimation of velocity, if the frictional properties of the system are known. Inversely, if the time evolution of the block velocity is known (e.g. measured), ∂F ∂δ can be determined:

∂F ∂δ = k v ∞ v ss -1 . ( 18
)
According to the previous paragraph, the steady-state slip is unstable, if and only if for

any δ 0 : ∂F ∂δ (δ0,0) < -k or ∂F ∂ δ (δ0,0) < -η. ( 19
)
The first condition refers to slip-weakening and coincides with the one described in [START_REF] Dieterich | Time-dependent friction and the mechanics of stick-slip[END_REF] (see also [START_REF] Goodman | Introduction to Rock Mechanics[END_REF]; [START_REF] Scholz | The mechanics of earthquakes and faulting[END_REF]). The so-called nucleation length [START_REF] Andrews | Rupture velocity of plane strain shear cracks[END_REF]) can be retrieved setting k = A 2 π λ+µ λ+2µ µ L ,where λ and µ = G are the Lamé constants (µ not to be confused with the friction coefficient in the next paragraphs)

and A the slip area. The second condition shows that the system is unstable when the friction is slip-rate weakening (velocity-weakening). However, for common parameters of frictional velocity weakening [START_REF] Reches | Fault weakening and earthquake instability by powder lubrication[END_REF] and viscosity of rocks [START_REF] Vutukuri | Introduction to Rock Mechanics[END_REF] this instability condition is not critical.

We assume Coulomb friction, F = µN , where N is the effective normal force applied on the block, as shown in Figure 1, and µ is the friction coefficient. The coefficient of friction varies from an initial value µ max (static friction coefficient), to a residual one µ res (kinetic friction coefficient). Figure 2a shows schematically the transition between static and kinetic friction. This transition is made in a characteristic distance D c and depends on the frictional properties of the fault system. In Figure 2b we show the in-

(a) (b)
F s = k(v ∞ -v ss )
is not high enough (F s < F ), then no sliding takes place and the system is stable. Therefore, the necessary condition for instability is F s = F . In this case we say that an (existing) fault is (re)activated. However, it is worth emphasizing that this condition is not a sufficient one for seismic slip (unstable behavior). If the aforementioned instability conditions ( 19) are not satisfied, we have slip, but this a slow, aseismic, creep-like slip with velocity given by Eq.( 18).

From the physical point of view, a dynamic instability takes place if the slip weakening is higher than the negative slope of the effective elastic response of the surrounding to the fault rocks or, in other words, when the elastic unloading of the surrounding rocks cannot be counterbalanced by fault friction. The same happens also when the friction shows velocity weakening that cannot be counterbalanced by the viscosity of the surrounding rock mass. However, the situation changes when the fluid pressure is not constant with time P f = P f (t). The question addressed in the next section is exactly how one should control P f (t) in order to assure stable slip.

3 Controling instabilities

Control system configuration

We assume a general negative-feedback control system as depicted in Figure 3. Σ(P )

is the system to be controlled, the spring-slider in our case, Σ(C) the controller we need to design, y(t), the output of the controlled system Σ(P, C), i.e. the displacement of the block δ(t), u(t) the input of Σ(P ), u c (t) the input of the controller Σ(C) and y c (t) its output. u 1 (t) and u 2 (t) are inputs to the system (e.g fluid pressure, long and short range perturbations, applied slip velocity etc.). We seek the controller Σ(C) that can stabilize the spring-slider model by modifying (controlling) the applied fluid pressure P f (t) (input). The problem is challenging due to friction and the consequent non-linearities it introduces. Additionally, the pore pressure multiplies the friction coefficient and does not allow us to write the mathematical system in canonical forms frequently used in the mathematical theory of control [START_REF] Vardulakis | Linear Multivariable Control: Algebraic Analysis and Synthesis Control[END_REF][START_REF] Vardulakis | Introduction to the mathematical theory of the theory of signals, systems and control[END_REF]. However, the target of the present work is to stabilize the system and stay in the vicinity of an evolving steady-state. This justifies the linearization of the equations in terms of slip, slip-rate and fluid pressure.

Linearization

We want to stabilize the fast dynamics of the system described by Eq.( 2). Therefore, we need to assure that the system is stable in the fast time scale described by the order ε 0 in the cascade system of equations ( 8). Setting y 0 = ỹ + y * 0 in Eq.( 8).1 and the same equation expressed at the reference state y * 0 , we obtain:

ÿ = -a(ỹ + y * 0 , ẏ, φ) + a(y * 0 , 0, φ * 0 ) -cỹ -d ẏ, (20)
where ỹ represents a perturbation of the reference steady-state solution, φ = φ(t) the dependence of the a(y, ẏ, φ) function in terms of an input (here the fluid pressure), φ(t), which takes the reference value φ * 0 at steady state. Expanding a(ỹ+y * 0 , ẏ, φ) in power series around ỹ = 0, ẏ = 0 and φ = φ * 0 (steady-state) we obtain:

ÿ + a 1 ẏ + a 0 ỹ = b 0 u, (21)
where

u(t) = φ(t) -φ * 0 , a 0 = α = ∂a ∂y (y * 0 ,0,φ * 0 ) + c, a 1 = β = ∂a ∂ ẏ (y * 0 ,0,φ * 0 ) + d and b 0 = -∂α ∂φ (y * 0 ,0,φ * 0 )
. In the following, the tilde over y is dropped for simplicity in notation.

Stabilizability

The above equation describes the behavior of the system for small perturbations from its steady-state. It is linear both in terms of y and the input u and its characteristic polynomial is:

D(s) = s 2 + a 1 s + a 0 , (22)
while for the input (right-hand-side) is:

N (s) = b 0 . ( 23
)
The transfer function of the system is:

P (s) = N (s) D(s) = b 0 s 2 + a 1 s + a 0 . ( 24
)
The poles of the system are at s 0 satisfying D(s 0 ) = 0. As expected, the poles have nonnegative real part when α or β are negative (instability). Since the system has no decoupling zeros in the closed right half complex plane, i.e. in C + = {s ∈ C, Re(s) ≥ 0}, the system is stabilizable. Decoupling zeros are called the common roots of N (s) and D(s) that are not roots of its transfer function P (s). Consequently, a stabilizing compensator (stabilizing controller) can be designed. This is a major result for the application at hand as it shows that earthquakes can be controlled, at least from the mathematical point of view.

Proper stabilizing controller

According to [START_REF] Vardulakis | Linear Multivariable Control: Algebraic Analysis and Synthesis Control[END_REF] the system Σ(C), with transfer function

C(s) = Y (s) X(s) ∈ R pr (s)
, is a stabilizing compensator, if and only if, the characteristic polynomial of the closed system Σ(P, C):

D c (s) = X(s)D(s) + Y (s)N (s), (25)
has all its roots in C -= {s ∈ C, Re(s) < 0}. Following the procedure described in [START_REF] Vardulakis | Linear Multivariable Control: Algebraic Analysis and Synthesis Control[END_REF][START_REF] Vardulakis | Introduction to the mathematical theory of the theory of signals, systems and control[END_REF] it is possible to determine Y (s) and X(s) and, therefore, design the desired stabilizing controller. Let Λ 5 = {λ 1 , λ 2 , λ 3 , λ 4 , λ 5 } the set of roots of D c (s) such that λ i = λ ∈ C -. The polynomials X(s) and Y (s) are determined by solving the following linear system:

ωT M4 = dT , (26)
where M4 is the Wolovich resultant of rank 6 [START_REF] Antoniou | On the computation and parametrization of proper denominator assigning compensators for strictly proper plants[END_REF]:

M 4 =                       a 0 a 1 1 0 0 0 b 0 0 0 0 0 0 0 a 0 a 1 1 0 0 0 b 0 0 0 0 0 0 0 a 0 a 1 1 0 0 0 b 0 0 0 0 0 0 0 a 0 a 1 1 0 0 0 b 0 0 0                       (27)
and dT = -λ 5 5λ 4 -10λ 3 10λ 2 -5λ 1 (28) ωT contains the coefficients of X(s) and Y (s):

ωT = χ 0 ψ 0 χ 1 ψ 1 χ 2 ψ 2 χ 3 ψ 3 (29)
such that:

X(s) = χ 3 s 3 + χ 2 s 2 + χ 1 s + χ 0 and Y (s) = ψ 3 s 3 + ψ 2 s 2 + ψ 1 s + ψ 0 . ( 30
)
The transfer function of the closed system Σ(P, C) is then:

H cl (s) =    P (s) 1+C(s)P (s) P (s)C(s) 1+C(s)P (s) -C(s)P (s) 1+C(s)P (s) C(s) 1+C(s)P (s)    (31)
and therefore the Laplace transform of the output of the closed system is:

   Y (s) Y c (s)    = H cl (s)    U 1 (s) U 2 (s)    , (32)
where U 1 (s) and U 2 (s) are the Laplace transforms of the inputs u 1 (t) and u 2 (t), respectively (see Figure 3).

Asymptotic tracking

We want to control the response of the system y(t) in order to asymptotically track a reference input u 2 (t) as t → ∞. In other words, we want the error function u c (t) = u 2 (t) -y(t) tend to zero for t → ∞. Consider the transfer function between the error function u c (t) and the reference input u 2 (t):

S(s) = U c (s) U 2 (s) = 1 1 + C(s)P (s) = X(s)D(s) X(s)D(s) + Y (s)N (s) , (33)
where U c (s) and U 2 (s) are respectively the Laplace transforms of u c (t) and u 2 (t). Let u 2 (t) be the step function, such that u 2 (t) = K for t ≥ 0 and zero for t < 0 or u 2 (t) = Kt for t ≥ 0 and zero for t < 0. Then according to [START_REF] Vardulakis | Introduction to the mathematical theory of the theory of signals, systems and control[END_REF], if Σ(P, C) is asymptotically stable, then lim t→∞ u c (t) = 0, if and only if, C(s) has a double pole at s = 0 or equivalently if and only if χ 0 = 0 and χ 1 = 0. Solving the linear system (26)

under these constraints we determine the stabilizing controller, which has the following coefficients:

χ 0 = 0, χ 1 = 0, χ 2 = -(a 1 + 5λ), χ 3 = 1 ψ 0 = - λ 5 b 0 , ψ 1 = 5λ 4 b 0 , ψ 2 = 5a 0 λ + a 0 a 1 -10λ 3 b 0 , ψ 3 = 5a 1 λ + a 2 1 -a 0 + 10λ 2 b 0 (34)
and therefore:

C(s) = Y (s) X(s) = ψ 3 s 3 + ψ 2 s 2 + ψ 1 s + ψ 0 (χ 3 s + χ 2 )s 2 . ( 35
)

Approximate frictional parameters

Knowing a priori the frictional parameters of a fault system is practically impossible. Various geophysical methods can give only approximate estimations. Given a controller designed as described in the previous paragraphs, we will investigate here the tolerance in the frictional parameters that guarantee stabilization. Let Σ(P ) be the system for which a stabilizing controller, Σ(C), was designed and Σ(P ) the real system that having different frictional parameters than Σ(P ). The characteristic polynomial of the closed system Σ(P , C) is:

D c (s) = X(s)D (s) + Y (s)N (s), (36)
Let also the frictional parameters of Σ(P ) be a 0 = a 0 + ∆a 0 = ∂a ∂y (y * 0 ,0,φ * 0 )

+ c, a 1 = a 1 +∆a 1 = ∂a ∂ ẏ (y * 0 ,0,φ * 0 ) +d and b 0 = b 0 -∆b 0 = -∂α ∂φ (y * 0 ,0,φ * 0 )
. Then Eq.( 36) becomes:

D c (s) = D c (s) + X(s)∆D(s) + Y (s)∆N (s), (37)
where ∆D(s) = ∆a 1 s + ∆a 0 and ∆N (s) = ∆b 0 . This system is stable when D c is a stable polynomial, i.e. when it has all its roots in C -. The stability of the polynomial can be explored using the Hurwitz matrix or approximately using asymptotic methods, which gives the required tolerances for ∆a 1 , ∆a 0 and ∆b 0 in function of λ or numerically for a given system.

4 Numerical examples

Geological setting and scaling laws

We consider a fault system at 5km depth. This is a common depth for many energyrelated human activities in the earth's crust and in the range of modern drilling technology. At 5km the normal to the fault stress is about σ n = 100MPa (σ n = W A ) and the fluid pressure p f = 50MPa (p f = P f A ), leading to an effective normal stress of σ n = σ n -p f = 50MPa. Notice that these values vary considerably with the tectonic configuration, i.e. they depend on whether the fault system is in extensional, compressional or strike-slip setting. The density of the rock is taken equal to ρ = 2500kg/m 3 , its apparent shear modulus equal to G = 30GPa and its apparent viscosity C = 10 5 MPa s.

Viscosity is not a well constrained quantity. The chosen value is several orders of magnitude lower than the viscosity of rocks at ambient temperature, in order to account for the earth's crust temperature gradient and the high deviatoric stresses during (pre-)seismic slip. We refer to [START_REF] Vutukuri | Introduction to Rock Mechanics[END_REF] for experimental results on viscosity and its dependence on temperature and high stresses.

The shear stress drop caused by an earthquake varies considerably over the rupture area. This is due to material heterogeneities, fault roughness, geometrical factors, multiphysical couplings, locking etc.. Nevertheless, interesting conclusions can be drawn if one considers the average shear stress drop over the whole rupture area, ∆τ . According to seismological inversions of actual earthquakes, ∆τ varies between 0.1 and 10 MPa.

Here, we take ∆τ = 5MPa. It is worth emphasizing that this is a spatial average of the shear stress drop over the whole area of the fault and that stress drop can be much higher locally. According to [START_REF] Kanamori | The physics of earthquakes[END_REF] the expected (average) slip is:

D = θ -1 G -1 L∆τ, ( 38
)
where θ is a geometric constant of order unity. The seismic moment is defined as

M 0 =
GDA, where A is the rupture area, which here is assumed circular, A ≈ L 2 . Consequently:

M 0 = θ -1 L 3 ∆τ. ( 39
)
The seismic magnitude M w is defined as follows:

M w = 2 3 log 10 M 0 -6.07 (M 0 in Nm) (40)
From the above scaling equations it is clear that the magnitude of an earthquake is logarithmically related to the length of the fault L, or in other words to the size of the fault's rupture area.

Assuming µ max = 0.6 and µ res = 0.5 such that ∆τ = σ n ∆µ = 5MPa the springslider model provides useful insights regarding earthquake instability. It models earthquake nucleation and seismic slip in terms of average quantities. However, it ignores the spatial rupture process and propagation, which will be studied in details in a future work.

Based on field measurements, we take D c = 10mm. Notice that this value is much higher than the D c measured in the laboratory [START_REF] Kanamori | The physics of earthquakes[END_REF]. Expressing Eq.( 1)

in the form of Eq.( 3) we obtain:

a = σ n ρL µ(δ) T 2 D , c = G ρL 2 T 2 , d = C ρL 2 T, (41)
where T and D are, respectively, any reference/characteristic time and length scales leading to Eq.( 2) such that ε = v∞ (Figure 2). We neglect velocity weakening as at low slip velocities it is several orders of magnitude lower [START_REF] Reches | Fault weakening and earthquake instability by powder lubrication[END_REF] than the apparent viscosity of the surrounding rock and, therefore, it does not influence instability (Eq.( 12)). Finally, we assume a far field tectonic movement of v inf = 1cm/year.

In Figure 4 we present the response of the spring-slider for L = 5km. The equations of the system were integrated using Wolfram Mathematica 11.2. According to Eq.( 38) the seismic slip is 0.83m and, as expected, it coincides with the final displacement obtained from the spring-slider model. The seismic moment is M 0 = 6.25×10 17 Nm, which corresponds to an earthquake of magnitude M w = 5.8. For smaller fault length L=500m, the seismic moment is M 0 = 6.25 × 10 14 Nm and the magnitude M w = 3.8. Figure 5 shows the evolution of slip and slip-rate for this fault length as obtained by the spring-slider model. Piece-wise linear or exponential evolution of µ has minor influence on the response of the system. tion, unconventional oil and gas production etc. (see [START_REF] Rubinstein | Myths and Facts on Wastewater Injection, Hydraulic Fracturing, Enhanced Oil Recovery, and Induced Seismicity[END_REF]).

Of course, in all these projects, avoiding seismic events is of paramount importance. We assume that the system is not on the verge of instability and that the actual shear stress due to the far field tectonic movement is lower than the critical one. Critical we call the average shear stress required to render the system unstable and provoke the earthquake event of M w = 3.8 described in the previous paragraph (Figure 5). Knowing the exact state of stress in a fault system is not trivial. In this scenario we assume that the real (average) shear stress state along the fault system of L=500m is at 90% of the critical shear stress, i.e. τ real = 0.9τ crit . At this stress state the system becomes unstable when the fluid pressure increases for ∆p f,crit = 5 MPa. However, due to inaccurate measurements and other uncertainties, the project coordinators have considered an average shear stress level lower than the real one. Based on this wrong estimation, they have decided to inject fluid leading to a maximum fluid pressure increase of ∆p f,applied = 10 MPa, which is double than ∆p f,crit and will cause a considerable seismic event without any control system. The injection program will last one week.

To mitigate this risk the fluid pressure increase is continuously regulated by the stabilizing compensator designed in section 3. We choose λ = -0.1. The connectivity of the controller assures negative-feedback to the fault system as shown in Figure 3. We use the friction parameters at δ = 0 for setting the controller (Eq.( 29)). More specifically, µ = 0.6, dµ dδ = -0.01 and a 0 = -3.5 × 10 -4 , a 1 = 1.6 × 10 -1 , b 0 = 4.8 × 10 -4

for L = 500 m (see Eq.( 41)).

Before the reactivation of the fault no slip is observed (locked). Nevertheless, when the critical fluid pressure is reached the system becomes unstable. The controller suc-manuscript submitted to JGR: Solid Earth ceeds in stabilizing the system assuring zero slip. No seismicity is observed and the M w = 3.8 earthquake is prevented.

Figure 6 depicts the total pressure change, which is automatically adjusted by the controller. Immediately after fault reactivation, the controller inhibits any further increase of the applied pressure, which according to the injection program should reach 10 MPa. Then it slowly reduces the pressure in order to keep the system stable. This reduction is barely perceptible, but it is necessary for the stabilization of the system. Notice that the far field tectonic velocity is always acting and it has a destabilizing effect as we are at the verge of instability. If the controller is deactivated the system will loose stability and an earthquake of M w ≈ 3.8 will take place. In this case we would like to drive safely the system from its unstable state to a stable one. This the objective of the next paragraph. normalized by the effective stress at that depth. Immediately after the fault reactivation at ∆p f = 5MPa, the controller cancels out the applied pressure and progressively reduces the fluid pressure (suction/pumping) in order to stabilize the system, which is now unstable.

Scenario #2: Controlling earthquakes

In this paragraph we present an example of active fluid pressure control for avoiding earthquakes, like the ones presented in the previous paragraphs. In particular we focus on avoiding the M w = 5.8 earthquake described in paragraph 4.1 (Figure 4).

Contrary to the previous scenario, in this scenario we don't want to avoid slip, but to induce it in a controllable way. Our target is to assure slow, aseismic slip such as to move the system from its initial unstable state to a stable one. In this way we will dissipate in a controllable manner the energy surplus and we will avoid the sudden energy release that leads to earthquakes.

The earthquake control operation is performed in a time window of ten minutes.

The technological feasibility of such an intervention is not discussed herein. We assume the worst case scenario, i.e. the system is on the verge of instability. However, it is constantly controlled by the stabilizing compensator we designed in Section 3. In this scenario, we choose to make the system slip for 0.83 m, which coincides with the seismic slip distance of the uncontrolled system (see Figure 4). Consequently, the target velocity that we want to apply is v c 8.3cm/min.

As far it concerns the controller, we choose λ = -0.1 and we calculate its parameters using the friction parameters at δ = 0 (Eq.( 29)). More specifically, µ = 0.6, dµ dδ = -0.01 and a 0 = -4.0 × 10 -5 , a 1 = 1.6 × 10 -3 , b 0 = 4.8 × 10 -5 for L=5km (see Eq.( 41)).

In Figure 7 we show the response of the controlled system by integrating the differential equations and connecting the controller as shown in Figure 3. We observe that no abrupt sliding takes place and that the system is successfully controlled (cf. Figure 4). The seismic event is avoided.

Figure 8 shows how the controller alters the fluid pressure from its initial hydrostatic value in order to drive the system from its initial unstable state to a stable one in an aseismic way. We observe that the controller reduces the pressure (pumping, suction), but in a non-monotonous way. In the beginning of sliding, the fluid pressure is decreased rapidly reaching a minimum value of approximately ∆p f = -9 MPa. This phase corresponds to the unstable phase of the system (critical distance, D c ). Then the fluid pressure is progressively increased and when we reach the target distance of the operation it recovers its initial hydrostatic value, i.e. ∆p f = 0 MPa. Notice that at this stage, the fault system is stable even in the absence of the controller ( dµ dδ = 0, see Eq.( 19)).

In Figure 9 we compare the kinetic energy of the controlled system with the uncontrolled one. The kinetic energy in the controlled system is four orders of magnitude lower than the one developed during the seismic event. Moreover, it is practically constant during the earthquake control operation.

Figure 10 shows the energy dissipation at the fault. The controlled system dissipates energy under almost constant rate, which depends on the chosen operation time.

Of course, this is not the case for the uncontrolled system where high dissipation rates are reported due to fast, seismic slip.

Finally, Figure 11 illustrates the drop of the elastic energy of the system. Again, the controlled system manages to reduce its potential energy smoothly (linearly in this example) and avoid its sudden release as in the case of the uncontrolled system. The above simulations were performed with a predetermined friction law that follows the exponential relation that was given in paragraph 4.1. However, it has to be emphasized that the parameters of the controller were kept constant during the simulation and they were not updated by taking into account the exact evolution of the frictional (left) and the uncontrolled one (right). In the controlled system the energy is dissipated progressively, while in the uncontrolled system the energy is dissipated abruptly (earthquake).

properties. This shows that the system is controlled even if the frictional properties of the system are not exactly known (see paragraph 3.6). In order to explore further the robustness of our approach we keep the same parameters as before for the controller and we add a sinusoidal perturbation to the friction coefficient as shown in Figure 12. Our stabilizing controller manages to control the system and avoid the earthquake instability despite the fluctuations of the frictional properties. In Figure 13 we show the fluid pressure change that the controller automatically adjusts to assure stability. The system is driven again from its unstable state to a stable one. The perturbation in the friction coefficient is reflected in the fluctuations of the fluid pressure.

Conclusions

The current paper presents a theoretical work focusing on exploring the possibility of preventing earthquakes by controlling fluid injection pressure. Our analysis is based on the classical spring-slider model (frictional slider), which we actively stabilize. We adopt not involve inertia effects related to the fast time scale at dominant order, ε 0 (Eqs. ( 8)).

This state represents the slow, creep-like, aseismic motion of the system and has a simple mathematical expression. More precisely, it is proportional to the far field tectonic velocity and inversely proportional to a term that involves the first derivative of the apparent friction with respect to slip and the apparent elastic stiffness of the surrounding the fault rocks (see Eqs.( 17) and ( 18)). The stability of this slow-slip motion is then investigated using Lyapunov's first method of stability and the conditions for which steadystate motion is stable are determined. Notice, that due to the far field tectonic velocity and the general rheology considered for friction, the system is non-autonomous. With this double-scale methodology, we manage to alleviate the explicit dependence of the system on the fast time scale, rendering it time-invariant at the dominant order.

Based on these results and the abovementioned time-invariance, the application of well established tools of the classical mathematical theory of control is possible. More specifically, control theory is used for stabilizing the system and make it remain in the vicinity of an evolving steady-state. This justifies the linearization of the equations in terms of slip, slip-rate and fluid pressure.

We show mathematically that the system is stabilizable by controlling fluid pressure. In other words, it is possible to stabilize and control the system when it is unstable, by appropriately adjusting the fluid pressure. This is a major result. The contrary would mean that earthquake control is impossible. Moreover, the opposite conclusion would imply that ongoing, large-scale industrial applications, involving injection of large quantities of fluids in the earth's crust, have high degree of risk, which cannot be mitigated.

Assuming a general negative-feedback control system, we designed a proper stabilizing controller. By monitoring slip, the designed controller succeeds in adjusting the fluid pressure and assures stable, aseismic slip, even in the absence of complete knowledge of the exact frictional properties of the system (robustness). Moreover, it succeeds in automatically controlling the pore fluid pressure and impose a prescribed slip or slip velocity (asymptotic tracking). These features of the controller allow a) to mitigate the seismic risk related to induced seismicity and b) to drive the system from an unstable state of high energy to a stable one of lower energy with constant slip velocity.

These features are illustrated through two scenarios of earthquake control. In these scenarios, fluids are injected under controlled pressure at 5km depth. More specifically, the first scenario refers to a fictitious injection project, where injection under constant pressure rate was planned. It is shown that the controller automatically stabilizes the system and avoids the anticipated earthquake event of M w = 3.8. In particular, the controller inhibits any further fluid pressure increase when the system enters to the unstable regime. This is done automatically by the controller by simply monitoring slip evolution. In this way the earthquake event is avoided. As the time increases, the controller progressively reduces the fluid pressure by pumping, in order to prevent seismic slip. Fluid removal reminds us the hypothetical scheme for EQ control proposed by [START_REF] Raleigh | An experiment in earthquake control at Rangely, Colorado[END_REF] and the related experiment at Rangely, Colorado in USA. However, with our approach, fluid pressure is automatically regulated in real-time, based on a designed negativefeedback control system. Consequently, more complex situations can be treated and stability can be actively assured without prescribing the fluid pressure history in advance.

The next scenario concerns the mitigation of a M w = 5.8 event by imposing constant slip velocity. For simplicity we assume that we are at the verge of unstable, seismic slip. The controller reduces the fluid pressure by pumping, but this is done in a nonmonotonous way. In the beginning of (imposed, desired) sliding, the fluid pressure is decreased rapidly reaching a minimum value of approximately ∆p f = -9 MPa. This phase corresponds to the unstable phase of the system (slip smaller that the critical distance, D c ). Then, the fluid pressure is progressively increased by the controller and, finally, it recovers its initial hydrostatic value. At this stage, the fault system is stable. During the controlled injection, the kinetic energy is kept constant and it is four orders of magnitude lower than the one that would develop abruptly, if the system was not controlled (seismic event). Regarding energy dissipation, the controller assures an almost constant dissipation rate, which depends on the chosen operation time. The same happens with the elastic energy, which is decreased smoothly over the duration of the EQ control operation. Of course, this is not the case for the uncontrolled system, where high dissipation rates and fast elastic energy drop are observed due to fast, sudden, seismic slip. Notice, that in the case of faults that are not at the verge of instability, the controller would increase the pore fluid pressure to enhance slip and then, once the system enters to the unstable regime, it would automatically reduce the pore-fluid pressure in order to stabilize the system and drive it to a stable state of lower energy, as mentioned above.

The above two scenarios show the ability of the designed controller to prevent earthquake events. Of course, these are academic examples and intensive research is needed before eventual real-scale applications. Notice that both physics and geometry of the fault system are kept as simple as possible. The role of heterogeneities, of pore fluid diffusivity and special hydrological conditions are not considered as well. Moreover, fault rupture and slip propagation is considered only on average over the whole fault's length. The investigation of the above mentioned limitations, as well as the observability of the real system and other techno-economical aspects of EQ control, exceed the scope of the present article and they are explored in the frame of the ongoing ERC project "Controlling earth-Quakes -CoQuake" (http://coquake.com).

Figure 1 .

 1 Figure 1. Spring-slider analog.

Figure 2 .

 2 Figure 2. Schematic representation of the evolution with slip of (a) the friction coefficient and (b) of the Coulomb friction force with and without fluid pressure, P f . The system becomes more ductile (lower slope) when the fluid pressure increases.

Figure 3 .

 3 Figure3. Negative feedback control system Σ(P, C). Σ(P ) is the system to be controlled with the controller Σ(C).

v ref 1 ,

 1 with v ref = D T , and y = δ D . The denominators ρL and ρL 2 are derived by considering that the mobilized mass, m, is equal to ρL 3 . Regarding friction evolution with slip, i.e. µ(δ), two cases are explored here. The first one is a piece-wise linear function of δ and the second the exponential expression µ(δ) = µ res 1 -∆µ µres e -δ Dc

Figure 4 .

 4 Figure 4. Evolution of slip and slip velocity in function of time for a fault length equal to 5km.

4. 2 Figure 5 .

 25 Figure 5. Evolution of slip and slip velocity in function of time for a fault length equal to 500m.

Figure 6 .

 6 Figure 6. Evolution of fluid pressure change. The total injection operation time was set to 1 week. On the left we present the fluid pressure change and on the right the fluid pressure change

Figure 7 .Figure 8 .

 78 Figure 7. Evolution of slip (left) and slip velocity (right) in function of time normalized withthe control operation time (here 10 minutes). The earthquake event is avoided and the system is successfully controlled.

Figure 10 .

 10 Figure 10. Comparison of energy dissipation due to friction between the controlled system

Figure 11 .Figure 12 .Figure 13 .

 111213 Figure 11. Comparison of the elastic energy drop between the controlled (left) and uncontrolled (right) systems. The designed controller avoids the sudden energy release that happens during an earthquake.

  Figure 9. Comparison of the kinetic energy of the controlled system (left) with the uncontrolled one (right). The designed controller manages to reduce four orders of magnitude the kinetic energy and keep it practically constant during the earthquake control operation.

	0.25 0.30 E kin [TJ]					4000 5000 E kin [TJ]
	0.20					3000
	0.05 0.10 0.15					1000 2000
	0.2	0.4	0.6	0.8	1.0	t/t operation	10	20	30	40	50	t [sec]
	300 000 400 000 500 000 W friction [TJ]					300 000 400 000 500 000 W friction [TJ]
	200 000					200 000
	100 000					100 000
	0.2 0.4 0.6 0.8 1.0	t/t operation	10	20	30	40	50	t [sec]

a general frictional law that accounts for slip and slip rate for the fault behavior. The fault is considered fully saturated and ideally oriented for slip in the ambient stress-field.

The rocks that surround the fault are considered as a Kelvin-Voigt material.

We describe the dynamics of the system using two time scales, i.e. a slow and a fast one. This scale separation technique allows to asymptotically approximate the solutions of the system as power series of a small parameter, ε, that expresses the ratio between the fast and the slow time scales. We define as steady-state, the motion that does -24-manuscript submitted to JGR: Solid Earth Acknowledgments I would like to thank A.I.Vardulakis for his insightful reviews and help in control theory. This work was supported by the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program (Grant agreement no. 757848 Co-Quake). This paper contains no data.