The critical slip distance in rate and state model for fault friction in the study of potential earthquakes can vary wildly from micrometers to few me- ters depending on the length scale of the critically stressed fault. This makes it incredibly important to construct an inversion framework that provides good estimates of the critical slip distance purely based on the observed ac- celeration at the seismogram. The framework is based on Bayesian inference and Markov chain Monte Carlo. The synthetic data is generated by adding noise to the acceleration output of spring-slider-damper idealization of the rate and state model as the forward model. 1. Introduction The quantification of fault slip is achieved using the Rate- and State- dependent Friction (RSF) model for friction evolution, which is considered the gold standard for modeling earthquake cycles on faults [1–6]. It is given by µ=µ0+Aln V V0 +Bln V0θ dc , dθ dt= 1−θV dc,(1) where V=|dd/dt|is the slip rate magnitude, a=dV dtwhich we hypothe- size is of the same order as recorded by seismograph, µ0is the steady-state friction coefficient at the reference slip rate V0,AandBare empirical di- mensionless constants, θis the macroscopic variable characterizing state of the surface and 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 conditions [7]. Regardless of the importance, 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 [8], to 0.1−5mas determined in numerical and seismological estimates based on geophysical observations [9], and further to several meters as determined in high-velocity laboratory experiments [10]. 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 [11]. With that in mind, we provide a framework in which synthetic earthquake data is used to quan- tify uncertainty in critical slip distance. While the resolution and coupled flow and poromechanics [12–26] associated with subsurface activity in the realm of energy technologies and concomitant earthquake quantification is a hot topic, in this work, we focus on the effect of a standard trigonometric perturbation with exponentially decreasing amplitude. 1.1. The forward model Figure 1: Spring Slider Damper Idealization of Fault Behavior As shown in Fig. 1, we model a fault by a slider spring system [27–29]. 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. The friction coefficient of the block is given by µ=τ σ=τl−kδ−ηV σ 2 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 [30] and ηis the radiation damping coefficient [31]. We consider the case of a constant stressing rate ˙ τl=kVl where Vlis the load point velocity, which we take to be Vl=V0(1 + exp ( −t/10)sin(10t)) (2) The stiffness is a function of the fault length land elastic modulus Eas k≈E l. With k′=E lσ, we get ˙µ≈k′(Vl−V)−k′′˙V (3) where k′′=η σ. Once the phenomenological form of ˙ µis known, we rewrite Eq. (1) as V=V0exp 1 A µ−µ0−Bln V0θ dc , ˙θ= 1−θV dc, ˙V=V A ˙µ−B θ˙θ(4) The ballpark values are: ✓Elastic modulus E= 5×1010Pa ✓Critical fault length l= 3×10−2m ✓Normal stress σ= 200 ×106Pa ✓Radiation damping coefficient η= 20×106Pa/(m/s) ✓A= 0.011 ✓B= 0.014 ✓V0= 1µm/s ✓θ0= 0.6 ✓µ0=µ0= 0.6 from which the effective stiffness and damping are obtained as ✓k′=E lσ=5×1010 3×10−2×2×108[1/m]≈10−2[1/µm] ✓k′′=η σ=2×107 2×108[s/m] = 10−7[s/µm ] 2. Bayesian inversion framework To define the inverse problem, considered the relationship between ac- celeration ( ai(t)) and the model response by the following statistical model a(ti) =f(ti, θ, µ, A, B, d c) +ϵi (5) 3 Figure 2: Noise is modeled as normal distribution where ϵiis the noise. Assuming the ϵi∼N(0, σ2) as unbiased, indepen- dent and identical normal distribution with standard deviation σallows us to conveniently generate the synthetic data as shown in Fig. 2. The syn- thetic data over time a(t1), ..., a (tn) are the nobservations for a(ti) and f(A, ti, θ, µ, A, B, d c) is the acceleration response of the model over time obtained using the forward model. The goal of the inverse problem is to determine the model parameter ( dc) from the Eq.(5) . Using the Bayes the- orem [32], the distribution for the model parameter is given by the posterior distribution π(dc|a(t1), ..., a (tn)) given by π(dc|a(t1), ..., a (tn)) =π(a(t1), ..., a (tn)|dc)π0(dc)R dcπ(a(t1), ..., a (tn)|dc)π0(dc)ddc(6) where π(a(t1), ..., a (tn)|dc) is the likelihood given by π(a(t1), ..., a (tn)|dc) =nY i=1π(a(ti)|dc) =nY i=11 σ√ 2πe−1 2a(ti)−f(V,ti,θ,µ,A,B,d c) σ2 (7) andπ0(dc) is the prior distribution. The denominator is a normalizing in- tegral. The information of the model parameter 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 4 the uniform distribution limits. The quantity that makes evaluation of the posterior difficult is the integral term in the denominator. Direct evaluation of the integral using quadrature rules is expensive and often requires adaptive methods. Alternatively, sampling methods like Markov chain Monte Carlo (MCMC) methods [33–36] can be used to evalute the integral. The MCMC method generates the parameter samples from a proposal dis- tribution, and the sample is rejected or accepted by evaluating the posterior distribution for the corresponding parameter sample. Adaptive Metropolis MCMC algorithm [35] is used as the sampling method. The MCMC simu- lation starts from a random value for the model parameter generated from uniform bounded distribution. The MCMC sampling method uses the ob- tained predicted acceleration and a new value for the model parameter is generated using the Markov chain. The errors are assumed to be a normal distribution with a standard deviation and at each iteration of the MCMC sampling method, a new value of standard deviation is generated, associated with the model parameter using inverse-gamma distribution. In the next iteration, the predictive model is solved again for the new parameter and the obtained predicted acceleration is used by the MCMC sampling method to generate the next parameter sample. A total of nMCMC samples are generated from the MCMC sampling method. 2.1. Results Figs. 3- 4 show a gradual improvement in the estimated model parameter as the number of samples is increased. It is evident that the framework has a lot more difficulty converging to the true value if the initial guess is far away from the true value. In reality, it is not really known what the true value is, and arbitrary initial guesses will lend a solution that is way off target regardless of the number of sampling points. We observe that the Bayesian framework settles to a value with maximum probability after a certain number of sampling points regardless of the number of samples following it. In all of the simulations, we observe in all the simulations that the Bayesian framework does an excellent job of dropping to the ballpark of the expected true value within the first few samples. We see that the framework is not heavily reliant on the initial guess for accuracy for the set of simulations we are rendering. 5 (a) 100 samples (b) 500 samples (c) 1000 samples (d) 2000 samples Figure 3: Output of the Bayesian MCMC framework for model parameter for different sample numbers with true value of dc= 1µmand an initial guess of 10 µm 3. Conclusions and Outlook The rate and state model for fault friction evolution is a piece of the puzzle in which forward simulations are used to arrive at the seismic impact of fault slip on earthquake activity. Typically, the earthquakes are measured with seismograms and geophones on the surface as P-waves and S-waves, and then these readings are used to calibrate the seismic activity for constant monitoring. The acceleration field around the fault slip activity translates to these waves recorded on the surface, and forward simulations with wave propagation bridge that gap. That being said, a seismic recording on the sur- face cannot easily be backtraced to acceleration around the fault, and finally to the source of the fault slip. That is precisely inverse modeling, and the Bayesian framework coupled with MCMC allows us to put such a framework 6 (a) 100 samples (b) 500 samples (c) 1000 samples (d) 2000 samples Figure 4: Output of the Bayesian MCMC framework for model parameter for different sample numbers with true value of dc= 1µmand an initial guess of 2 µm in place. The accelerations are field quantities, and the inverse estimation of the field around the fault from the time series at the seismogram and/or geophone is not a trivial task. With that in mind, we test the robustness of the Bayesian/MCMC framework to inversely estimate a value instead of a field. The thing is that any inversion framework works on data as the input, and since subsurface data is not available other the sensors at the wells, the forward simulations are used to generate this data with all the computational physics put in place. This generated data is then used to inversely estimate the acceleration field from the data at the surface. Although such forward simulations to generate data are infeasible in real-time scenarios where we need estimates of what is happening in the subsurface from the reading on the surface almost immediately, the framework robustness would eventually 7 lend itself to that scenario. The scenario is that the recording at the seismo- gram and/or geophone would be fed into the Bayesian/MCMC fraemwork as an input, and the framework would provide an estimate of the accelera- tion field around the fault as the output. In this work, we work on arriving at estimates of critical slip distance in the rate and state model by running a spring-slider-damper as the forward model instead of using a full-fledged forward simulator. In the future, we will be deploying the coupled flow and geomechanics simulator as the forward model, and gradually arrive at the aforementioned scenario.