This work presents a framework to inversely quantify uncertainty in the model parameters of the friction model using earthquake data via the Bayesian inference. The forward model is the popular rate- and state- friction (RSF) model along with the spring slider damper idealization. The inverse model is to determine the model parameters using the earthquake data as the response of the RSF model. The conventional solution to the inverse problem is the deterministic parameter values, which may not represent the true value, and quantifying uncertainty in the model parameters increases conﬁdence in the estimation. The uncertainty in the model parameters is estimated by the posterior distribution obtained through the Bayesian inversion. 1 Introduction Earthquakes occur as a result of global plate motion. Some earthquakes stop after only a few hundred meters while others continue rupturing for a thousand kilometers. An earthquake produces P waves, or compressional waves, that travel faster and reach the seismograph ﬁrst, and S waves, or shear waves, that are slower (Fig. 1). Both are transmitted within the Earth and are called body waves. Even slower are surface waves that run along the surface of the earth and do a lot of the damage. The earthquake focus is the point within the Earth where the earthquake originates. The epicenter is a point on the surface directly above the focus. The simplest model for earthquake initiation is to assume that when the stress accumulated in the plates exceeds some failure criterion on a fault plane, an earthquake happens [1]. Evaluating this criterion requires both a measure of the resolved stress on the fault plane and a quantiﬁable model for the failure threshold. The groundbreaking work of [2] arrived at the hypothesis that faulting occurs when the resolved shear stress exceeds the internal friction on some plane in the medium leading to fault slip. The quantiﬁcation of earthquakes from the fault slip is achieved using the Rate- and State-dependent Friction (RSF) model, which is considered the gold standard for modeling earthquake cycles (interseismic loading followed by coseismic relaxation) on mature faults [3–7]. It is given by =0+AlnV V0 +BlnV0 dc ; d dt= 1V dc;(1) whereV=jdd=dtjis the slip rate magnitude, a=dV dtwhich we hypothesize is of the same order as recorded by seismograph, 0is the steady-state friction coefﬁcient at the reference slip rate V0,AandBare A framework to quantify uncertainty in critical slip distance in rate and state friction model for earthquakes empirical dimensionless constants, is the macroscopic variable characterizing state of the surface and dcis a critical slip distance. Here, may be understood as the frictional contact time [3], or the average maturity of contact asperities between the sliding surfaces [8]. The evolution of is assumed to be independent of changes in the normal traction that can accompany the fault slip due to changes in ﬂuid pressure. The model accounts for the decrease in friction (slip-weakening) as the slip increases, and the increase in friction (healing) as the time of contact or slip velocity increases. The two effects act together such that A>B leads to the strengthening of the fault, stable sliding and creeping motion, and A<B leads to weakening of the fault, frictional instability, and accelerating slip. In this way, the model is capable of capturing repetitive stick-slip behavior of faults and the resulting seismic cycle [4, 6]. Figure 1: P wave arrives ﬁrst, followed by the S wave and then by surface waves.The critical slip distance, dc, is the distance over which a fault loses or regains its frictional strength after a perturbation in the loading conditions [9]. In principle, it determines the maximum slip accel- eration and radiated energy during an earthquake insofar that it inﬂuences the magnitude and time scale of the associated stress breakdown process (e.g., fracture energy) [10]. Regardless of the impor- tance, it is paradoxical that the values of dcreported in the literature range from a few to tens of microns as determined in typical laboratory experiments with bare surfaces and gouge layers [10], to 0.1–5 m as determined in numerical and seismological esti- mates based on geophysical observations [11], and further to several meters as determined in high- velocity laboratory experiments [12]. Note that among these studies, the critical length scale pa- rameter in the constitutive friction laws ( dc) may differ from the slip-weakening distance inferred from the traction evolution curves obtained for nat- ural or laboratory faults. The latter, as usually de- rived from scenarios where perturbations are large (e.g., velocity steps of large magnitudes or tips of dynamic rupture nucleation zones), is also referred to as the equivalent or effective slip-weakening dis- tance (deq 0ord0, see a review by [13]). Moreover, in most numerical simulations of dynamic rupture propagation with prescribed friction laws, dcis im- posed a priori and its value is often assumed to be constant and uniform on the fault plane. Understanding the physics that controls the critical slip distance and explains the gap between observations from experimental and natural faults is thus one of the crucial problems in both the seismology and laboratory communities. With that in mind, we provide a framework in which the earthquake data is used alongside RSF to quantify uncertainty in critical slip distance. While the resolution and coupled ﬂow and geomechanics [14 –25] associated with subsurface activity in the realm of energy technologies and concomitant earthquake quantiﬁcation is a hot topic, in this work, we focus on the effect of stress perturbations in the absence of pore pressure variable. In section 2, we explain the spring slider damper idealization to infer the inﬂuence of critical slip distance on RSF without recourse to complicated elastodynamic equations. In section 3, we explain the Bayesian inference framework to inversely quantify uncertainty in the estimation of critical slip distance. In section 4, we present conclusions and outlook for future work. 2 A framework to quantify uncertainty in critical slip distance in rate and state friction model for earthquakes 2 Spring slider damper idealization to study earthquake response We ﬁrst rewrite Eq. (1) as V=V0exp 1 A 0Bln V0 dc ; _= 1V dc; =_V dc9 >>= >>;(2) As shown in Fig. 2, we model a fault by a slider spring system [26 –28]. The slider represents either a fault or a part of the fault that is sliding. The stiffness krepresents elastic interactions between the fault patch and the ductile deeper part of the fault, which is assumed to creep at a constant rate. This simple model assumes that slip, stress, and friction law parameters are uniform on the fault patch.  Figure 2: Spring Slider Damper Idealization of Fault Behavior The friction coefﬁcient of the block is given by = = lkV  whereis the normal stress, the shear stress on the interface, lis the remotely applied stress acting on the fault in the absence of slip, - kis the decrease in stress due to fault slip [1] and is the radiation damping coefﬁcient [29]. We consider the case of a constant stressing rate _ l=kVlwhereVlis the load point velocity. The initial stress may be smaller or larger than steady state friction owing to coseismic slip on the fault patch or on adjacent parts of the fault. The expression neglects inertia, and is thus only valid for low slip speed in the interseismic period. The stiffness is a function of the fault length land elastic modulusEaskE l. Withk0=E l, we get _k0(VlV)k00_V; k0(_Vl_V)k00V (3) wherek00= . Once the phenomenological form of _andis known, we use the following to get _Vand_a, _V=V A _B _ ; _a=_V A _B _ +V A B +B 2_9 = ;(4) 2.1 Forward Model Response to a Standard Impulse We follow the steps outlined in Algorithm 1 to arrive at the temporal variations of acceleration and fault friction coefﬁcient. We initialize the friction coefﬁcient and state variable and obtain the slip rate and rate of change of the state variable. We then use these values to obtain time derivatives of acceleration and slip rate. These time derivatives are required as we employ the integrated feature of the scientiﬁc Python package SciPy [30]. The inﬂuence of critical slip distance on system response to a load point perturbation of the form Vl= 1 + exp (t=20)sin(t=10) 3 A framework to quantify uncertainty in critical slip distance in rate and state friction model for earthquakes Algorithm 1 Rate and state friction model with radiation damping term Initialize=0,=ref Use Eq. (2) to get V,_and Use Eq. (3) to get _and Use Eq. (4) to get _Vand _a 0.600.62Friction 010State Variable 2.55.0Velocity 0 10 20 30 40 50 Time [sec]−250Accelerationdc= 1/uni03BCμ (a) 0.6000.605Friction 510State Variable 1.01.5Velocity 0 10 20 30 40 50 Time [sec]−10Accelerationdc= 5/uni03BCμ (b) 0.5990.6000.601Friction 9.9510.00State Variable 1.01.1Velocity 0 10 20 30 40 50 Time [sec]−10Accelerationdc= 10/uni03BCμ (c) 0.5950.600Friction 1015State Variable 01Velocity 0 10 20 30 40 50 Time [sec]0.02.5Accelerationdc= 20/uni03BCμ (d) Figure 3: System response for different values of critical slip distance. Units of displacement, velocity and accelearation are m,m=s andm=s2 is shown in Fig. 3. The code to generate the plots has been given in Appendix A. This code is a part of the GitHub repository https://github.com/karthikncsu/ Bayesian-inference-using-earthquake-data . The ballpark values are taken from [1] and [29]. Elastic modulus E= 51010Pa, Critical fault length l= 3102m, Normal stress = 200106Pa, Radiation damping coefﬁcient = 20106Pa=(m=s), 4 A framework to quantify uncertainty in critical slip distance in rate and state friction model for earthquakes A= 0:011andB= 0:014. The effective stiffness and damping is obtained as k0=E l=51010 31022108[1=m]104[1=m]102[1=m]; k00= =2107 2108= 0:1[s=m]1107[s=m ] 3 Bayesian Inversion Framework 0 5 10 15 20 25 30 35 time (Sec)-0.4-0.3-0.2-0.100.10.20.3acceleration (g) Figure 4: Seismograph reading (measured as normalized to g= 9:8m=s2) of 1940 EL Centro earthquake of magnitude Mw= 6:9that occurred in the Imperial Valley in southeastern Southern California near the USA-Mexico border. It was the ﬁrst major earthquake to be recorded by a strong-motion seismograph located next to a fault rupture, and led to a total damage of $6million [31]. Fig. 3 shows the response of the forward (rate and state friction) model for a given point load perturbation (Vl), critical slip distance parameter ( dc), and empirical constants ( 0,V0,AandB). The response of the model is the acceleration of the slider, computed using algorithm 1. In an inverse problem, the acceleration response of the model is known and the goal is to ﬁnd the parameter, critical slip distance parameter ( dc). The earthquake acceleration data shown in Fig. 4 can be considered as the acceleration response of the model. To deﬁne the inverse problem, considered the relationship between acceleration ( ai(t)) and the model response by the following statistical model a(ti) =f(ti;;;A;B;d c) +i (5) whereiis the error in the statistical model. Here the a(ti)andiare the random variables. The earthquake data over time a(t1);:::;a (tn)are thenobservations for a(ti)andf(A;ti;;;A;B;d c)is the acceleration response of the model over time obtained using the Algorithm 1. The goal of the inverse problem is to determine the model parameter ( dc) from the Eq. (5)and conventional method to determine the model parameter that mimizes the norm of the errors using the least squares ﬁt solution as shown below dc;0= arg min dcnX i=1(i)2= arg min dcnX i=1(a(ti)f(ti;;;A;B;d c))2(6) The critical slip distance parameter ( dc) obtained using the least-squares ﬁt solution, Eq. (6)is deterministic value. The values estimated using the least square ﬁt are not the true values due to inherent noise in the 5 A framework to quantify uncertainty in critical slip distance in rate and state friction model for earthquakes data and in most cases, the noise in the data makes it difﬁcult to ﬁnd the true value. Instead, ﬁnding a probability distribution for the model parameters encompasses the true model parameter values and increases the conﬁdence in the prediction. Using the Bayes theorem [32], the distribution for the model parameters is given by the posterior distribution (dcja(t1);:::;a (tn)) =(a(t1);:::;a (tn)jdc)0(dc)R dc(a(t1);:::;a (tn)jdc)0(dc)ddc(7) Here(dcja(t1);:::;a (tn))is the posterior, (a(t1);:::;a (tn)jdc)is the likelihood and 0(dc)is the prior distribution for the model parameters. Assuming the iN(0;2)as unbiased, independent and identical normal distribution with standard deviation , the likelihood function is expressed as (a(t1);:::;a (tn)jdc) =nY i=1(a(ti)jdc) =nY i=11 p 2e1 2a(ti)f(V;ti;;;A;B;dc) 2 (8) In above equation, the f(V;ti;;;A;B )is calculated using the forward problem given by algorithm 1. The information of the model parameters can be included in the posterior distribution through the prior, 0(dc). In this study, the prior is assumed to be uniform distribution and the prior is a constant value inside the uniform distribution limits. Algorithm 2 Bayesian inference Input data: earthquake data, a(t1);:::;a (tn) Generate grid dc Use Algorithm 1 to get RHS of Eq. (5) Use (8) to get (a(t1);:::;a (tn)jdc)for each grid point. Integrate Eq. (7) to obtain the posterior distribution The goal of the inverse problem is to calculate the posterior distribution Eq. (7), which represents the uncertainty in the critical slip distance parameter ( dc) due to the noise in the earthquake data. Direct evaluation of the posterior distribution using quadrature rules is expensive and often requires adaptive methods to ﬁnd the posterior distribution. Alternatively, sampling methods like Markov chain Monte Carlo (MCMC) methods [32–40] can be used to generate samples from the posterior distribution. 4 Conclusions and Future Work This work presents a framework to inversely quantify uncertainty in the critical slip distance parameter of the rate and state friction (RSF) model via the Bayesian inference using the earthquake data. The forward model is to determine the acceleration, using the RSF model, for the given model parameters. In case of an inverse problem, the acceleration data is known and the goal is to ﬁnd the model parameters. Using conventional methods such as least-squares methods, a deterministic value of the critical slip distance parameter can be obtained from the inverse problem. However, the deterministic parameter value estimated using the conventional methods does not represent the true values due to the noise in the earthquake data, and quantifying uncertainty in the model parameters increases the conﬁdence in the prediction. The uncertainty in the model parameter is estimated by the posterior distribution obtained from the Bayes theorem. The future work will be to demonstrate a simulation to quantify uncertainty in the critical slip distance parameter using the earthquake data via sampling methods such as Markov chain Monte Carlo. 6 A framework to quantify uncertainty in critical slip distance in rate and state friction model for earthquakes A Python Program to Generate Forward Model Response 1from scipy.misc import derivative 2import numpy as np 3from scipy import integrate 4import matplotlib.pyplot as plt 5from math import exp,log,pi,sin,cos 6 7fig = plt.figure() 8fig.suptitle(’$d_c=20 \mu m$’) 9 10def friction(t,y): 11 k = 1e-2 12 a = 0.011 13 b = 0.014 14 Dc = 20. 15 mu_ref = 0.6 16 V_ref = 1. 17 k1 = 1e-7 # radiation damping term 18# k1 = 0 19 20# if t < 5: 21# V_lp = 1. 22# else: 23# V_lp = 10. 24 25 amp_t = 20 26 trig_t = .1 27 tp = t/amp_t 28 tt = t/trig_t 29 30 amp = 1 31 V_lp = 1 + amp *exp(-tp) *sin(tt) 32 dV_lp = -amp/amp_t *exp(-tp) *sin(tt) + amp/trig_t *exp(-tp) *cos(tt) # time derivative of V_l 33 34 # Just to help readability 35 #y[0] is mu (friction) 36 #y[1] is theta 37 #y[2] is velocity 38 #y[3] is acceleration 39 40 n = len(y) 41 dydt = np.zeros((n,1)) 42 43 # compute v 44 temp_ = V_ref *y[1] / Dc 45 temp = 1/a *(y[0] - mu_ref - b *log(temp_)) 46 v = V_ref *exp(temp) 47 48 # time derivative of theta 49 dydt[1] = 1. - v *y[1] / Dc 50 51 # double derivative of theta 52 ddtheta = - dydt[1] *v/ Dc 53 54 # time derivative of mu 55 dydt[0] = k *V_lp - k *v 56 57 # time derivative of velocity 58 dydt[2] = v/a *(dydt[0] - b/y[1] *dydt[1]) 59 60 # double derivative of mu 61 ddmu = k *dV_lp - k *dydt[2] 7 A framework to quantify uncertainty in critical slip distance in rate and state friction model for earthquakes 62 63 # time derivative of acceleration 64 dydt[3] = dydt[2]/a *(dydt[0] - b/y[1] *dydt[1]) + v/a *(ddmu - b/y[1] *ddtheta + b/y[1] *dydt [1]/y[1]) 65 66 # radiation damping 67 dydt[0] = dydt[0] - k1 *dydt[2] 68 dydt[2] = v/a *(dydt[0] - b/y[1] *dydt[1]) 69 ddmu = ddmu - k1 *dydt[3] 70 dydt[3] = dydt[2]/a *(dydt[0] - b/y[1] *dydt[1]) + v/a *(ddmu - b/y[1] *ddtheta + b/y[1] *dydt [1]/y[1]) 71 72 return dydt 73 74r = integrate.ode(friction).set_integrator(’vode’, order=5,max_step=0.001,method=’bdf’,atol=1 e-10,rtol=1e-6) 75 76# Time range 77t_start = 0.0 78t_final = 50. 79delta_t = 1e-2 80num_steps = int(np.floor((t_final-t_start)/delta_t)+1) 81 82# Initial conditions 83mu_t_zero = 0.6 84V_ref = 1. 85Dc = 10. 86mu_ref = 0.6 87theta_t_zero = Dc/V_ref 88v = V_ref 89start_acc = 0 90r.set_initial_value([mu_t_zero, theta_t_zero, V_ref, start_acc], t_start) 91 92# Create arrays to store trajectory 93t = np.zeros((num_steps,1)) 94mu = np.zeros((num_steps,1)) 95theta = np.zeros((num_steps,1)) 96velocity = np.zeros((num_steps,1)) 97acc = np.zeros((num_steps,1)) 98t[0] = t_start 99mu[0] = mu_ref 100theta[0] = theta_t_zero 101velocity[0] = v 102acc[0] = 0 103 104# Integrate the ODE(s) across each delta_t timestep 105k = 1 106while r.successful() and k < num_steps: 107 #integrate.ode.set_f_params(r,velocity,k) 108 r.integrate(r.t + delta_t) 109 110 # Store the results to plot later 111 t[k] = r.t 112 mu[k] = r.y[0] 113 theta[k] = r.y[1] 114 velocity[k] = r.y[2] 115 acc[k] = r.y[3] 116 k += 1 117 118# Make some plots 119 120#cjm_t,cjm_mu = np.loadtxt(’cjm_step.tim’,skiprows=2,unpack=True) 121 122ax1 = plt.subplot(411) 123ax1.plot(t, mu,color=’r’, linewidth=0.5) 8 A framework to quantify uncertainty in critical slip distance in rate and state friction model for earthquakes 124ax1.set_xlim(t_start, t_final) 125#ax1.plot(cjm_t+10.,cjm_mu,color=’k’) 126#ax1.set_xlabel(’Time [sec]’) 127ax1.set_xticklabels([]) 128ax1.set_ylabel(’Friction’) 129ax1.grid(’on’) 130 131ax2 = plt.subplot(412) 132ax2.plot(t, theta, ’r’, linewidth=0.5) 133ax2.set_xlim(t_start, t_final) 134#ax2.set_xlabel(’Time [sec]’) 135ax2.set_ylabel(’State Variable’) 136ax2.set_xticklabels([]) 137ax2.grid(’on’) 138 139ax3 = plt.subplot(413) 140ax3.plot(t, velocity, ’r’, linewidth=0.5) 141ax3.set_xlim(t_start, t_final) 142#ax3.set_xlabel(’Time [sec]’) 143ax3.set_ylabel(’Velocity’) 144ax3.set_xticklabels([]) 145ax3.grid(’on’) 146 147ax4 = plt.subplot(414) 148ax4.plot(t, acc, ’r’, linewidth=0.5) 149ax4.set_xlim(t_start, t_final) 150ax4.set_xlabel(’Time [sec]’) 151ax4.set_ylabel(’Acceleration’) 152ax4.grid(’on’) 153 154#ax5 = plt.subplot(515) 155#ax5.plot(t, forcing, ’r’) 156#ax5.set_xlim(t_start, t_final) 157#ax5.set_xlabel(’Time [sec]’) 158#ax5.set_ylabel(’Forcing’) 159#ax5.grid(’on’) 160 161plt.show() Listing 1: Python example