The burgeoning need to sequester anthropogenic CO 2for climate mitigation and the need for energy sustenance leading upto enhanced geothermal energy production has made it incredibly critical to study potential earthquakes due to ﬂuid activity in the subsurface. These earthquakes result from reactivation of faults in the subsurface due to pore pressure perturbations. In this work, we provide a framework to model fault slip due to pore pressure change leading upto quantifying the earthquake magnitude. 1 Introduction 0 5 10 15 20 25 30 35 time (Sec)-0.4-0.3-0.2-0.100.10.20.3acceleration (g) (a) 5 10 15 20 25 30 T (s)0.10.20.30.40.50.60.70.80.911.1displacement (m) (b) Figure 1: 1940 El Centro earthquake of magnitude Mw= 6:9. (a) Accelerogram (measured as normalized to acceleration due to gravity g= 9:8m=s2) and (b) Displacement Response Spectrum Earthquakes occur as a result of global plate motion. Some plate boundaries glide past each other smoothly, while others are punctuated by catastrophic failures. Some earthquakes stop after only a few hundred A fault slip model to study earthquakes due to pore pressure perturbations A P REPRINT metres while others continue rupturing for a thousand kilometres. 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] started with the fact that any stress ﬁeld can be completely described by its principal stresses, which are given by the eigenvectors of the stress tensor and are interpretable as the normal stresses in three orthogonal directions. He then proposed that: (1) the stress state could be resolved by assuming that one principal stress is vertical since the Earth’s surface is a free surface and (2) faulting occurs when the resolved shear stress exceeds the internal friction on some plane in the medium. Internal friction is deﬁned analogously with conventional sliding friction as a shear stress proportional to the normal stress on a plane. One complication to this simple picture was recognized early on. High ﬂuid pressures can support part of the load across a fault and reduce the friction. The importance of the ﬂuid effect on fault friction was ﬁrst recognized by [ 3]. In the course of their work on oil exploration, they observed that pressures in pockets of ﬂuids in the crust commonly exceeded hydrostatic pressure. They connected this observation with studies of faulting and proposed that the pore pressure at a depth can approach the normal stress on faults, resulting in low friction. Fig. 1 shows the seismograph recorded acceleration of the 1940 EL Centro earthquake that 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 [ 4]. In section 2, we explain how the pore pressure enters into the computation of fault slip and we explain the friction model. In section 3, we explain the characteristics of earthquake response using an idealization and most popular empirical relations for quantifying earthquakes. Finally, in section 4, we provide conclusions and outlook for future work. 2 Fault slip model p p+ 0  n  + 0 +  Figure 2: The two sides of the fault surface, which need not be planar, are designated as the ‘ +’ side and the ‘’ side We treat faults as surfaces of discontinuity embedded in the continuum across which displacement is allowed to be discontinuous to recognize the possibility of fault slip [ 5,6]. As shown in Fig. 2, the fault normal vector npoints from the negative side to the positive side of the fault. Slip on the fault is the displacement of the positive side relative to the negative side, (u+u)d=0on; (1) wheredis the fault slip vector. We impose the effective traction on the fault by introducing a Lagrange multiplier (see [ 7] for a clear understanding of Lagrange multipliers), l, which is a force per unit area required to satisfy Eq. (1). A positive value of lnindicates that a tensile effective stress is transmitted across the fault surface. The Kuhn-Tucker conditions of contact mechanics are obeyed such that no penetration occurs and the effective normal traction stays compressive at the contact surface. The frictional stress fon 2 A fault slip model to study earthquakes due to pore pressure perturbations A P REPRINT the fault as f= cln;ln<0; c;ln0;(2) where cis the cohesive strength of the fault and is the coefﬁcient of friction. The shear traction on the fault is computed as =jl(ln)nj We use the Mohr-Coulomb theory to deﬁne the stability criterion for the fault [ 8]. When the shear traction on the fault is below the friction stress,  f, the fault does not slip. When the shear traction is larger than the friction stress, > f, the contact problem is solved to determine the Lagrange multipliers and slip on the fault, such that the Lagrange multipliers are compatible with the frictional stress [5]. A difference in ﬂuid pressure across the fault leads to a pressure jump JpK=p+p, wherep+andp are the equivalent multiphase pressures on the ‘positive’ and the ‘negative’ side of the fault. This pressure jump leads to a discontinuity in the effective stress across the fault, such that the total stress is continuous, 0 nbpn=0 +nbp+n; (3) wherebis the Biot coefﬁcient. Fault stability can be assessed by evaluating the stability criterion on both sides of the fault separately. The side of the fault where the criterion is met ﬁrst determines the fault stability. Equivalently, this can be achieved by deﬁning a fault pressure that is a function of the pressures on the two sides, p+andp. Introducing the fault pressure allows us to uniquely deﬁne the effective normal traction on the fault, 0 n, and determine the fault friction f(Eq. (2)). Since the stability criterion > fis ﬁrst violated with the larger pressure, we deﬁne the fault pressure pas p= max (p;p+): (4) The pressures on either side of the fault are obtained from the in-house coupled multiphase ﬂow and geomechanics simulator [ 5,9] which employs a staggered solution algorithm (see [ 10,11,12,13,14,15,16]) to solve the coupled system of equations. 2.1 Friction model The rate- and state-dependent friction model [17, 18, 19, 20, 21] is =0+AlnV V0 +BlnV0 dc ; d dt= 1V dc;(5) whereV=jdd=dtjis the slip rate magnitude, 0is the steady-state friction coefﬁcient at the reference slip rateV0,AandBare 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 [17], or the average maturity of contact asperities between the sliding surfaces [ 22]. 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 increase. The two effects act together such that A>B leads to strengthening of the fault, stable sliding and creeping motion, andA < 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 [ 18,20]. Eq. (5) can be written as V=V0exp1 A 0BlnV0 dc ; d dt= 1V dc(6) 3 A fault slip model to study earthquakes due to pore pressure perturbations A P REPRINT 3 Spring slider idealization to study earthquake response Following [23, 24, 25], we model a fault by a slider spring system. 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, that 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. 3.1 Without radiation damping term The friction coefﬁcient of the block is given by = = lk (7) wherekis the spring stiffness, is the normal stress, the shear stress on the interface, lis the remotely applied stress acting on the fault in the absence of slip, and - kis the decrease in stress due to fault slip. We consider the case of a constant stressing rate _ l=kVl (8) whereVlis 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. Expression (7)neglects inertia, and is thus only valid for low slip speed in the interseismic period. The stiffness is a function of the fault lengthland shear modulus G[25] kG l(9) In lieu of the above, we write d dtk0(VlV); dV dt=Vd dt1 A 0BlnV0 dc V Ad dtB d dt V A k0(VlV)B d dt(10) where k0=G l(11) Algorithm 1 Rate and state friction w/o radiation damping term Initialize=0,=ref Use Eq. (6) to get V,d dt Use Eq. (10) to getd dtanddV dt 3.2 With radiation damping term The friction coefﬁcient of the block is given by = = lkV (12) In lieu of the above, we write dV dt= 1 +Vk00 A1V A k0(VlV)B d dt d dtk0(VlV)k00dV dt(13) 4 A fault slip model to study earthquakes due to pore pressure perturbations A P REPRINT Algorithm 2 Rate and state friction with radiation damping term Initialize=0,=ref Use Eq. (6) to get V,d dt Use Eq. (13) to getdV dtandd dt 3.3 Plots Units of velocity and fault slip are m=s andmrespectively. k0= 0:1[1=m],k00= 1[s=m ],dc= 10[m], 0= 0:6,V0= 1[m=s ]and the load point velocity taking a jump as follows Vl= 1[m=s ]; t< 5[s]; 10[m=s ]; t5[s] 0.590.600.61Friction 510State Variable 510Velocity 0 10 20 30 40 50 Time [sec]0250Fault slipA<B (a) Slip weakening behavior without radiation damping 0.590.60Friction 510State Variable 510Velocity 0 10 20 30 40 50 Time [sec]0250Fault slipA<B (b) Slip weakening behavior with radiation damping Figure 3:A= 0:005andB= 0:01 3.4 Quantifying earthquake The total amount of fault slip during the event is correlated to seismic magnitude. The seismic moment M0 is a measure of the size and energy release of an earthquake [ 26].M0is deﬁned as the integral of fault slip over the fault area Stimes the shear modulus as M0=GZ jdjdSGjdjS (14) Although fault surface area can seldom be measured, several empirical relations make it possible to estimate it for faults for which the average slip is known. [ 27] successfully explained several empirically determined scaling laws of earthquakes using a model in which the surface area of slip is proportional to the square of average slip S=G2 4ec2c3jdj (15) wherec22104GPa2m1for fault slip values ranging from 1mto105m[28],c31=2ande2[29]. From Eqs. (14) and (15), M0G3 81014jd2j (16) 5 A fault slip model to study earthquakes due to pore pressure perturbations A P REPRINT 0.600.62Friction 510State Variable 510Velocity 0 10 20 30 40 50 Time [sec]0250Fault slipA>B (a) Slip strengthening behavior without radiation damping 0.600.61Friction 510State Variable 510Velocity 0 10 20 30 40 50 Time [sec]0250Fault slipA>B (b) Slip strengthening behavior with radiation damping Figure 4:A= 0:01andB= 0:005. From [30], the seismic moment magnitude Mwis deﬁned as Mw=log10M0 1:56:06 (17) whereM0is deﬁned in Nm. For a shear modulus of 20GPa and fault slip of jdj= 2m, the seismic moment magnitude is M081030 8101441:21017 Mwlog10(1:21017) 1:56:06 = 5:3 4 Conclusions and outlook We provided a framework in which the pore pressure changes are likely incorporated into the quantiﬁcation of earthquakes triggered as a result of ﬂuid injection and/or withdrawal in the subsurface. The rate and state model is the standard bearer for the friction evolution with fault slip. It is important to note that the slider spring model would be replaced by the output of the coupled ﬂow and geomechanics simulator in our future works. The slider spring model, though, provides an excellent gauge for the potential earthquake response of the subsurface to stress perturbations. It is also important to note that the constant stress rate loading Eq. (8)is valid for low slip speeds in the interseismic period. It can be replaced by some other phenomenological models to study earthquake response for high slip speeds in the interseismic period. It is also important to note that the computation of fault shear stress can involve the use of the radiation damping term, especially in the context of enhanced geothermal energy production [ 31]. We shall study all these models as a part of future work.