We present an algorithmic framework to solve an inverse problem using Bayesian inference and Markov chain Monte Carlo sampling. The input of the inverse problem is the acceleration of the slipping seismogenic fault and the output is the probability distribution of the critical slip distance parameter of the rate and state model for fault friction. 1 Introduction Earthquakes are triggered after critically stressed seismogenic faults slip [1] after being activated due to natural causes like earth tides [2], rainfall [3], snowfall [4], typhoons [5], as well as human activity [6]. The component linking fault activation and fault slip is the popularly used rate- and state- dependent friction (RSF) model for friction evolution [7–9] given by =0+Aln V V0 +Bln V0 dc _= 1V dc(1) whereV=jdd=dtjis the slip rate magnitude, 0is the steady-state friction coefﬁ- cient at the reference slip rate V0,AandBare empirical dimensionless constants, is the macroscopic variable characterizing state of the surface and ﬁnally dcis a critical slip distance. The critical slip distance is the distance over which a fault loses or regains its frictional strength after a perturbation in the loading condi- Bayesian inference and Markov chain Monte Carlo based estimation of critical slip distance parameter in rate and state model for fault friction tions [10]. In principle, it determines the maximum slip acceleration and radiated energy to such an extent that it inﬂuences the magnitude and time scale of the associated stress breakdown process [11]. Regardless of the importance, it is para- doxical that the values of dcreported in the literature range from a few to tens of microns as determined in typical laboratory experiments [11], to 0:15mas de- termined in numerical and seismological estimates based on geophysical observa- tions [12], and further to several meters as determined in high-velocity laboratory experiments [13]. Moreover, in most numerical simulations of dynamic rupture propagation with prescribed friction laws, dcis imposed a priori and its value is often assumed to be constant and uniform on the fault plane. 1.1 The inverse problem Understanding the physics that controls the critical slip distance and explains the gap between observations from experimental and natural faults is one of the cru- cial problems in both the seismology and laboratory communities [14]. Further- more, inferring the causality between fault activation and fault slip is important for seismic design and monitoring of underground structures [15], bridges [16], rein- forced concrete buildings [17], carbon sequestration projects [18] and enhanced geothermal systems [19]. In lieu of that, an inverse problem is formulated with the fault slip acceleration being the input and critical slip distance being the output. We consider the relationship between acceleration ai(t)and the model response by the following statistical model a(ti) =F(ti;dc) +i (2) whereiis the error in the statistical model. Here the a(ti)andiare the random variables. The synthetic data over time a(t1);:::;a (tn)are thenobservations for a(ti)andF(ti;dc)is the acceleration response of the model over time obtained us- ing the forward model. The goal of the inverse problem is to determine the model parameterdcfrom the Eq. (2) using Bayesian inference [20] with Markov chain Monte Carlo (MCMC) [21–27] to generate samples from the posterior distribu- tion. 2 Bayesian Inversion and MCMC Framework The following components are present XdataD Xmodel M 2 Bayesian inference and Markov chain Monte Carlo based estimation of critical slip distance parameter in rate and state model for fault friction Xparameters XpriorP(jM) XlikelihoodP(Dj) XposteriorP(jD) A key part of Bayesian inversion is to express the posterior using Bayes’ theo- rem [28] given by P(jD) =P(Dj)P(jM)R P(Dj)P(jM)d(3) A tractable property of Bayesian inversion is that small changes in data lead to small changes in the posterior even when the inverse problem is ill-posed in the classical sense, so Bayesian inversion is stable. The posterior is however quite complicated with no closed form expression, so much of the contemporary re- search focuses on realizing the aforementioned advantages with Bayesian inver- sion without requiring access to the full posterior. 2.1 Prior P(jM)is often referred to as the prior. This describes the probability of having a particular set of values for our given model M before conditioning on our data. Because this is independent of the data, this term is often interpreted as represent- ing our “prior beliefs” about what the parameterss should be based on previous measurements, physical concerns, and other known factors. In practice, this has the effect of essentially “augmenting” the data with other information. 2.2 Likelihood P(Dj)is often referred to as the likelihood. The data likelihood is often given from insight into the physics of how data is generated (simulator). 2.3 Posterior Finally,P(jD)represents our posterior. This quantiﬁes our belief in parame- ters after combining our prior intuition P(jM)(prior) with current observations P(Dj)(likelihood). The posterior will be some compromise between the prior and the likelihood, with the exact combination depending on the strength and prop- erties of the prior and the quality of the data used to derive the likelihood. 3 Bayesian inference and Markov chain Monte Carlo based estimation of critical slip distance parameter in rate and state model for fault friction 2.4 Markov chain Monte Carlo The calculation of these quantities reduces to the evaluation of the following inte- gralZ h()P(jD)d under various choices of the function h. The idea behind MCMC is that one can produce a simulated sample f(1);(2);:::;(M)gfrom this posterior density leading to an evaluation of the integral: M1MX j=1h((j))M!1!Z h()P(jD)d The formal deﬁnition of MCMC is that it is a method to sample a given multivari- ate distribution by construcing a suitable Markov chain with the property that its limiting, invariant distribution, is the target distribution . 2.4.1 Metropolis-Hastings algorithm Algorithm 1: Metropolis-Hastings algorithm 1Specify an initial value (0); 2forj=1,2,...,M do 30q((j)); 4 ((j);0) =( minh (0)q(0;(j)) ((j))q((j);0);1i if ()q((j);0)>0 1 Otherwise; 5(j+1)= 0if Unif (0;1) ((j);0) (j)Otherwise 6Returnf(1);(2);:::;(M)g The Metropolis-Hastings (M-H) method is a general MCMC method to produce sample variates from a given multivariate density. As shown in Algorithm 1, it is based on a candidate generating density that is used to supply a proposal value and a probability of move that is used to determine if the proposal value should be taken as the next item of the chain. The probability of move is based on the ratio of the target density (evaluated at the proposal value in the numerator and the current value in the denominator) times the ratio of the proposal density (at the current value in the numerator and the proposal value in the denominator). The method is extremely general and powerful, it being possible in principle to view almost any MCMC algorithm, in one way or another, as a variant of the M-H algorithm. 4 Bayesian inference and Markov chain Monte Carlo based estimation of critical slip distance parameter in rate and state model for fault friction 3 Bayesian applied to our problem We deﬁne the probability distribution of the critical slip distance, dc, as a posterior distribution for the given acceleration data and the predictive model. The para- metric posteriors distribution arrived from the Bayes theorem is as shown below (dcja(t1);:::;a (tn)) =(a(t1);:::;a (tn)jdc)0(dc)R dc(a(t1);:::;a (tn)jdc)0(dc)ddc(4) Herenis the number of data points, (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(ti;dc) 2 (5) The data points a(t1);:::;a (tn)in Eq. (5) are generated from the forward model explained next 3.1 The forward model  Figure 1: Spring slider damper idealization of fault activation Firstly, the RSF model Eq. (1) is rewritten as V=V0exp1 A 0BlnV0 dc (6) to obtain _Vas _V=V A _B _ (7) 5 Bayesian inference and Markov chain Monte Carlo based estimation of critical slip distance parameter in rate and state model for fault friction To obtain _, we deploy the spring slider damper model as shown in Fig. 1. The stiffnesskrepresents elastic interactions between the fault patch and the ductile deeper part of the fault, and is a function of the fault length land elastic modulus E askE l. This simple model assumes that slip, stress, and friction law parameters are uniform on the fault patch. The friction coefﬁcient of the block is given by = = lkV (8) 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. In lieu of Eq. (8), we obtain the time derivative of the friction coefﬁcient as _=k0(VlV)k00_V (9) wherek0=E landk00= . From Eqs. (1), (7) and (9), we get _V=A V+k001 k0(VlV)B  1V dc (10) In lieu of Eqs. (6) and (9), the forward model is written as V=V0exp 1 A 0Bln V0 dc F(t;dc) =A V+k001 k0(VlV)B  1V dc (11) 3.2 System response for ballpark values We employ the load point velocity variation of the form Vl=Vref(1 + exp (t=10)sin(10t)) (12) The ballpark values are taken from [1] and [29] as follows: elastic modulus E= 51010Pa, critical fault length l= 3102m, normal stress = 200106Pa, radiation damping coefﬁcient = 20106Pa=(m=s),A= 0:011andB= 0:014. The effective stiffness and damping is obtained as k0102[1=m]andk00= 107[s=m ]respectively. The reference values are Vref= 1m=s ,0= 0:6and ref=0= 0:6. As shown in Appendix A, we employ SciPy [30] to arrive at acceleration time series, which is plotted for critical slip distance parameter values ranging over three orders of magnitude in Fig. 2. 6 Bayesian inference and Markov chain Monte Carlo based estimation of critical slip distance parameter in rate and state model for fault friction 0 10 20 30 40 50 Time [sec]−10−5051015Accelerationdc= 0.1 /uni03BCμ 0 10 20 30 40 50 Time [sec]−6−4−202468Accelerationdc= 1/uni03BCμ 0 10 20 30 40 50 Time [sec]−1.00−0.75−0.50−0.250.000.250.500.75Accelerationdc= 10/uni03BCμ 0 10 20 30 40 50 Time [sec]−0.075−0.050−0.0250.0000.0250.0500.075Accelerationdc= 100 /uni03BCμ Figure 2: System response for different values of critical slip distance over three orders of magnitude. Units of acceleration are m=s2 4 Conclusions and Outlook We presented a Bayesian inference and MCMC based framework to arrive at a probability distribution of the critical slip distance parameter in the RSF model with the fault slip acceleration as the input to the inverse problem. We deployed the popularly used SSD model to generate data sets of the fault slip acceleration. While the solution of inversion problem with the fault slip acceleration as the input and a parameter of RSF model is a step forward, the eventual goal is the solution of an inverse problem with the seismograph reading of the earthquake as the input and parameters of the appropriate forward model as the output. It can be argued that the fault slip acceleration is of the same order of magnitude as the seismo- graph reading, but such an abstraction would be a gross over-simpliﬁcation of the problem at hand. The slipping fault generates a stress ﬁeld which spawns off seismic waves which then reach the seismograph typically mounted on the earth’s surface. Any forward model which encompasses the domain going all the way to the earth’s surface [31] and considers the impact of pore pressure perturbations on fault activation [32] would be more appropriate for the subsequent inverse prob- lem. The computational framework popularly used to study the impact of pore 7 Bayesian inference and Markov chain Monte Carlo based estimation of critical slip distance parameter in rate and state model for fault friction pressure perturbations on geotechnics is the ﬁxed stress split staggered solution algorithm [31–40]. Furthermore, the effect of wave propagation needs to be con- sidered for faults slipping at moderate to high speeds to accurately capture the induced seismicity. A two-way coupled poroelastodynamics forward model in the presence of faults with the parameters identiﬁed for a subsequent inversion using Bayesian and MCMC with the seismograph readin as the input is the way forward. Be that as it may, it is important to identify the problem statement at hand without necessarily resorting to computational overkill. A Python code to generate acceleration time series 1d e f f r i c t i o n ( t , y ) : # y [ 0 ] i s mu , y [ 1 ] i s t h e t a , y [ 2 ] i s v e l o c i t y 2 dy dt = np . z e r o s ( ( l e n ( y ) , 1 ) ) 3 power = 1 / a *( y [0] − mu_ref −b *l o g ( V_ref *y [ 1 ] / Dc ) ) 4 v = V_ref *exp ( power ) # Compute v 5 dy dt [ 1 ] = 1. − v *y [ 1 ] / Dc # Compute t im e d e r i v a t i v e of t h e t a 6 dy dt [ 0 ] = kprime *V_l − kprime *v 7 dy dt [ 2 ] = v / a *( dy dt [0] − b / y [ 1 ] *dy dt [ 1 ] ) 8 dy dt [ 0 ] = dy dt [0] − k1 *dy dt [ 2 ] # Compute ti me d e r i v a t i v e of mu 9 dy dt [ 2 ] = v / a *( dy dt [0] − b / y [ 1 ] *dy dt [ 1 ] ) # Compute ti me d e r i v a t i v e of v e l o c i t y 10 r e t u r n d yd t 11 12r = i n t e g r a t e . ode ( f r i c t i o n ) . s e t _ i n t e g r a t o r ( ’ vode ’ , o r d e r =5 , max_step = 0 . 0 0 1 , method= ’ bdf ’ , a t o l =1e −10 , r t o l =1e −6) # Use i n t e g r a t e f u n c t i o n a l i t y w i t h i n Python package SciPy 13r . s e t _ i n i t i a l _ v a l u e ( [ 0 . 6 , Dc / V_ref , V_ref ] , t _ s t a r t ) # I n i t i a l c o n d i t i o n s on mu , t h e t a , v e l o c i t y 14k = 1 15w h i l e r . s u c c e s s f u l ( ) and k < num_steps : 16 r . i n t e g r a t e ( r . t + d e l t a _ t ) # I n t e g r a t e t h e ODE( s ) a c r o s s each d e l t a _ t t i m e s t e p 17 t [ k ] , mu[ k ] , t h e t a [ k ] , v e l [ k ] = r . t , r . y [ 0 ] , r . y [ 1 ] , r . y [ 2 ] 8 Bayesian inference and Markov chain Monte Carlo based estimation of critical slip distance parameter in rate and state model for fault friction 18 acc [ k ] = ( v e l [ k ] − v e l [ k − 1 ] ) / d e l t a _ t # Compute t h e a c c e l e r a t i o n 19 k += 1 Listing 1: Python code for forward response of spring slider damper idealization