Physics-Based Approach to Deep Interseismic Creep: Implications for North Tabriz Fault Behavior Using MCMC

: Many studies assume that the distribution of a fault slip rate remains constant with time when examining surface deformation. However, recent numerical simulations have shown that dynamic rupture can penetrate regions with increased friction and diffuse from the lock-to-creep transition, contradicting this assumption. Bruhat and Segall (2017) introduced a new method to account for the downward penetration of interseismic slip into the locked zone. This study builds upon their work by applying their model to strike-slip fault environments and incorporating creep coupling to viscoelastic ﬂow in the lower crust and upper mantle. In this study, using Bruhat’s (2020) model, the interseismic deformation rates on the North Tabriz Fault are investigated. This study utilizes elastic and viscoelastic probabilistic models to ﬁt horizontal surface rates. By employing this updated approach, a physics-based solution for deep interseismic creep is developed, revealing potential slow vertical propagation. The improved ﬁt of horizontal deformation rates on the North Tabriz Fault is examined, leading to reasonable estimations of earthquake rupture depth and seismic displacement. The best-ﬁt solutions suggest a half-space relaxation time of approximately 156 years, with a diffusion rate of less than 1 m/year and around 0.419 m/year, indicating minimal creep diffusion.


Introduction
In the initial phase of investigating interseismic deformation, fault depiction relied on a solitary screw dislocation within an elastic half-space, indicating the fault's immobilization at a certain depth with a consistent reduced velocity [1].This simplified model enabled researchers to use kinematic inversions of geological surface rates to estimate the locking depth, a parameter presumed to define the seismic zone's depth range and that was utilized for decades.However, the screw dislocation model's physical realism is constrained by infinite stress concentration at the dislocation tip.To address this, the development of more realistic models incorporating transition zones between fully locked faults and freely creeping areas becomes imperative.These zones involve the gradual smoothing of creeplocked slip distribution, crucial for mitigating stress singularity.Consequently, many inversions for interseismic slip rate include smoothing or linear transition from locked to creeping faults (e.g., [2]) to alleviate stress singularity, regardless of the implementation speed.Nonetheless, comprehension of the governing physical properties in these transition zones remains limited [3].
Fully dynamic earthquake cycle simulations suggest that the transition from locking to creep within fault systems can experience significant changes over the course of the seismic cycle.The velocity-amplifying region's enhanced dynamic attenuation behavior is believed to potentially aid the propagation of dynamic rupture during the post-earthquake shift Eng.Proc.2023, 56, 26 2 of 8 from locking to creep [4].Pseudodynamic simulations, incorporating thermal stress effects, have additionally confirmed that in the interseismic periods between major earthquakes, the extent of slow slip events gradually expands into the locked zone [5].A novel method for characterizing interseismic slip rate, proposed by [6], introduces the concept of slip penetration upslope within the locked zone.This approach treats deep interseismic slip as a loaded crack with a constant slip rate at its downslope end, offering analytical expressions for stress drop, slip, and slip rate along the fault.These expressions permit the extension of non-singular slip rate distributions through Chebyshev polynomials, thus simplifying the inversion of fault interface-related physical properties and initiating a connection between kinematic inversions and physics-based earthquake cycle simulations.
The current study applies the deep interseismic creep propagation model to the north Tabriz strike-slip fault.Unlike ref. [6], who solely considered creep propagation within a fully elastic environment, ref. [3] incorporates long-term deformation resulting from viscoelastic flow in the lower crust and upper mantle, along with potential cumulative effects from previous earthquake cycles.Consequently, surface predictions undergo significant changes.Typically, fully elastic models tend to indicate a locking depth exceeding the seismic depth.However, incorporating viscoelastic effects (e.g., [7,8]) allows for a more reasonable fit to interseismic deformation rates, resulting in shallower locking depths (e.g., [9], Section 12.4.2).While physically motivated models have recently emerged (e.g., [10][11][12][13][14]), they remain scarce and computationally demanding.From the study of ref. [3], a simple and efficient method for implementing kinematic inversions, particularly within the Bayesian framework, is introduced.This method offers an improved physical representation of deep interseismic creep.
The study conducted by [3] investigates the application of probabilistic elastic and viscoelastic models to accurately fit horizontal surface rates.Building upon the model proposed in [6], the study incorporates the interaction between fault creep and viscoelastic flow, resulting in an enhanced understanding of deep interseismic creep.The updated approach provides a physics-based solution that reveals potential slow vertical propagation.This study assesses the improvement in capturing the appropriate horizontal deformation rate for the North Tabriz Fault.
Rizza et al. in 2013 [15] have determined slip rate values for the North Tabriz fault through the utilization of two geodetic observation techniques.Based on data collected from GPS stations, the slip rate is reported as 7.3 mm/year, whereas radar interferometric analysis indicates a rate of 6.0 mm/year.In a separate study [16], the estimated the slip rate for this fault is 7.3 mm/year using a block model approach.Additionally, ref. [17] conducted kinematic finite element modeling to evaluate land surface deformations in the Iranian plateau.The modeling incorporated various data inputs, including fault geometry, slip rate estimations based on geological and geomorphological studies, GPS velocity field during significant earthquake intervals, and stress directions and velocity values at the boundary points of the modeling region.This modeling study estimates a sliding rate of 8.5 mm/year for the northern fault of Tabriz.
The proximity of the Tabriz metropolis to the North Tabriz fault underscores the significance of investigating the mechanical interaction of this fault to analyze seismic risk in the city.In this study, a Monte Carlo Markov chain tool and Bayesian statistical method are employed to estimate the values of eight parameters (full rupture depth, the downdip limit of the transition region (elastic thickness), long-term rate (slip rate), coseismic displacement, relaxation time, the depth of uniform creep, locking depth, and free parameter to account for rigid block motion (α)) for the North Tabriz fault).A total of 1,000,000 simulations were conducted to obtain these estimates.

Tectonic Setting (Study Area)
Iran is positioned within the oblique collision zone of the Arabian and Eurasian tectonic plates, experiencing significant intracontinental deformation that shapes its prominent topography.The collision between these plates drives a convergence rate of 22 ± 2 mm/year, with the northwestern part of Iran, contributing 8 ± 2 mm/year to this rate [18].This convergence is manifested through processes like shortening in the Zagros and Alborz Mountains, as well as internal deformation from strike-slip faults in central Iran.The region's complex geodynamic system is closely linked to the interaction between the Arabian, Anatolian, and Eurasian plates, incorporating features like the North Anatolian Fault, Eastern Anatolian Fault, and Caucasus Mountains, acting as boundaries for the Zagros Mountains.The intricate fault system facilitates the transfer of northward motion from the Arabian Plate to the Anatolian Plateau.The Zagros Mountains' collision zone inclination enables displacement release between Caucasus shortening and the North Tabriz fault's dextral strike-slip movement [19].Notably, ref. [20] segmented the northern Tabriz fault into several sections based on historical earthquake surface ruptures from 1780, 1721, and 1786, highlighting key segments east and west of Tabriz city.The geometry of the fault displays misaligned sections in linear steps, with two significant segments showing evidence of past surface faulting and right-lateral movement [21].Furthermore, the presence of clear evidence of surface faulting resulting from past earthquakes is discernible through morphogeostructural complexities and paleoseismological investigations [22].

Data
In this study, a 1 × 1 grid was utilized for the inversion of the complete fault plane located to the north of Tabriz.The geometric attributes that were taken into account are provided in Table 1.
Table 1.In the simplest scenario, the North Tabriz fault was characterized using average values.The table includes the fault segment name in the first column, UTM coordinates of the fault segment's starting point in the second and third columns, strike and dip of the fault in the fourth and fifth columns, and fault length and width in the sixth and seventh columns, respectively.Furthermore, in this investigation, the initial values required for the Monte Carlo Markov chain procedure were determined by averaging data from previous studies, as outlined in Table 2. Considering the historical record of significant earthquakes in Tabriz, it was determined that the most recent major earthquake took place in 1786.As a result, the value for the t component, representing the time elapsed since the last earthquake, was set at 237 years.The prior limits for the eight parameters (full rupture depth, the downdip limit of the transition region (elastic thickness), long-term rate (slip rate), coseismic displacement, relaxation time, the depth of uniform creep, locking depth, and free parameter to account for rigid block motion (α)) mentioned are consolidated and presented in Table 3, based on previous research findings.

Methodology 4.1. Forward Problem
In the initial phase of this investigation, the boundary element method proposed by [23] is employed to perform a forward problem solution on the North Tabriz Fault.By incorporating various characteristics of the fault, such as its length, width, dip, strike, and slip rate, we are able to extract the GPS velocity data along the profile that are perpendicular to the fault.The primary objective of employing the forward problem approach is to obtain accurate GPS velocities on the perpendicular profile to the fault, thereby addressing the necessary corrections associated with the sphericity of the Earth, the Euler pole correction, and the three-dimensional corrections.
To simulate the displacements resulting from slip rates in the fault investigated here, we utilize the analytical model proposed by [24], commonly known as the Okada model, which is based on the theory of dislocation.The Okada model employs two categories of input parameters: physical and geometrical.The physical parameters employed in this model consist of the Lame coefficients µ and λ specific to the study area, which should be known approximately.In the absence of precise values, global average values obtained from sensitivity analyses of the Okada model can be utilized.The geometrical parameters considered in the Okada model encompass length, width, locking depth, dip, strike, slip rate, the coordinates of the starting point of the fault, and observation point coordinates.By incorporating the fault's geometry and the physical characteristics of the study area, this model converts fault displacement or slip rate into a displacement field or velocity field resulting from these parameters.

A Physics-Based Approach to Deep Interseismic Creep
In the subsequent phase of this study, the approach proposed by [3] was employed to estimate the interseismic deformation rates in the North Tabriz fault system.This method incorporates a viscoelastic earthquake cycle model, accounting for the upward propagation of deep interseismic creep.By employing this approach, we determined a total of 8 parameters associated with the North Tabriz fault system.

Slip Rate Inversion
The surface velocities, denoted as V horz , are derived from the combined influence of the viscoelastic seismic cycle V EQcycle , as well as the elastic and viscoelastic responses arising from interseismic creep, represented by V elcreep and V vecreep , respectively.
In Equation ( 1), the parameter α represents the dissimilarity between the reference frames of the fault model (typically antisymmetric in fault cases) and the measured velocities.The covariance matrix of the data, denoted as Σ, incorporates the statistical relationships among the variables.
V elcreep is connected to the slip rate vector .s, which is not known, by means of Green's functions G in a homogenous elastic half-space.In this study, we aim to find the values of several parameters through inversions, including the elastic thickness H, rupture depth D, present-day position of the locking depth (characterized by the top of the creep zone) d, long-term plate motion rate v ∞ , viscoelastic relaxation time t R , and maximum earthquake displacement ∆u, which is associated with the earthquake recurrence time T. Additionally, to account for the possibility of a constant creep in the transition region between the elastic and viscoelastic zones, we perform an inversion to determine the lower depth of the creep zone, H creep .It is important to note that the fault slides at a constant speed .δ ∞ between H and H creep .
The slip distribution of earthquakes can be characterized as follows.From the Earth's surface down to the complete rupture depth D, the seismic slip distribution is equivalent to the maximum seismic displacement, ∆u.In order to ensure that the total slip along the fault corresponds to the maximum seismic displacement at the end of a seismic cycle, the slip distribution between the top of the creep zone, d, and the lower depth of the creep zone, H creep , is defined as the complement of the slip distribution at the end of the cycle.This implies that the slip distribution, T = ∆u v ∞ , is integrated not only up until the current time but also up until the completion of the entire cycle.Similarly, the propagation speed of the up-dip rupture front, v up , is constrained so that the slip in the elastic region matches ∆u at the end of the cycle.Moreover, it is essential that the creeping zone reaches the lower boundary of the earthquake zone slope by the conclusion of the earthquake cycle.
The primary objective of this investigation is to develop inverse techniques for evaluating various models of interseismic deformation, some of which take into account the propagation of deep creep.In order to account for viscoelastic deformation, several parameters including the maximum rupture depth (D), elastic thickness (H), viscoelastic relaxation time (t R ), seismic displacement (∆u), long-term plate movement rate (v ∞ ), and α are incorporated into the calculations.When examining models involving deep interseismic slip, we also perform an inversion for the locking depth (d) and infer the propagation rate of v up .For these inversions, Markov chain Monte Carlo (MCMC) methods are employed.MCMC algorithms efficiently estimate the solution with the maximum likelihood and facilitate the construction of posterior distributions.Depending on the specific inversions conducted, prior knowledge about other model parameters will be integrated.The limits of the prior information are summarized in Table 3.

Result
In this investigation, the initial step involves employing the forward problem approach.This entails utilizing the specified geometry for the North Tabriz fault, as stated in Table 1, along with the average slip rate obtained from previous research, as presented in Table 2. To solve the forward problem, we adopt the boundary element approach proposed by [23].
By implementing this method, we successfully reconstruct the GPS velocity field along the profile perpendicular to the North Tabriz fault.The resulting vectors of the GPS velocity field, obtained through the application of the forward problem, are depicted in Figure 1.In the subsequent phase of this study, we utilize the acquired GPS ve input data for conducting an inversion process.The inversion is performe a physics-based approach and a Markov chain Monte Carlo (MCMC) too of this inversion, which assess the compatibility of geodetic rates in the presented in Figure 2. In Figure 2, the posterior distributions of various cluding the complete rupture depth (), elastic thickness (), slip rate ( placement (∆) , relaxation time ( ) , recurrence time () , locking dep creep depth ( ), and propagation velocity ( ), are depicted.The ried out through a substantial number of simulations, specifically, 1,000,0 During this inversion process, the estimated range for the rupture d 10 to 15 km.The elastic thickness is determined to range from 22 to 38 km displacement measures between 2 and 5 m.The recurrence time spans years.The slip rate is observed to be between 3 and6 mm per year, while ing depth, where a majority of the microseismicity occurs, is estimated t and 15 km.Furthermore, the propagation velocity and creep propagatio 0 to 8.5 m/year.The relaxation time for the Best fitting solution is deter years, with a range of 140 to 250 years.Figure 2 showcases the posteri representing the estimated parameters obtained through the creep diff along with the resulting propagation velocity.In the subsequent phase of this study, we utilize the acquired GPS velocity vectors as input data for conducting an inversion process.The inversion is performed by employing a physics-based approach and a Markov chain Monte Carlo (MCMC) tool.The outcomes of this inversion, which assess the compatibility of geodetic rates in the Tabriz fault, are presented in Figure 2. In Figure 2, the posterior distributions of various parameters, including the complete rupture depth (D), elastic thickness (H), slip rate (v ∞ ), seismic displacement (∆u), relaxation time (t R ), recurrence time (T), locking depth (d), uniform creep depth H creep , and propagation velocity v up , are depicted.The inversion is carried out through a substantial number of simulations, specifically, 1,000,000 iterations.
During this inversion process, the estimated range for the rupture depth falls within 10 to 15 km.The elastic thickness is determined to range from 22 to 38 km.The coseismic displacement measures between 2 and 5 m.The recurrence time spans from 450 to 750 years.The slip rate is observed to be between 3 and6 mm per year, while the current locking depth, where a majority of the microseismicity occurs, is estimated to be between 12 and 15 km.Furthermore, the propagation velocity and creep propagation rate vary from 0 to 8.5 m/year.The relaxation time for the Best fitting solution is determined to be 156 years, with a range of 140 to 250 years.Figure 2 showcases the posterior distributions, representing the estimated parameters obtained through the creep diffusion inversion, along with the resulting propagation velocity.
The optimal solutions obtained from the analysis suggest a complete rupture depth of 12.7 km, an elastic thickness of 23.1 km, a slip rate of 5.34 mm/year, a coseismic displacement of 4.36 m/year, and a relaxation time of 156 years.Additionally, the locking depth is estimated to be 12.9 km, the uniform creep depth is 20.8 km, and the propagation velocity is determined to be 0.419 m/year.

Conclusions
This study investigates models incorporating viscoelastic flow and deep interseismic creep to elucidate the deformation rate occurring within the Tabriz fault.The approach employed in this research utilizes the methodology developed by [3,6] to effectively characterize the coupling between the viscoelastic half-space and the time-dependent interseismic creep.The initial stage involved conducting a forward problem utilizing the boundary element method to accurately reconstruct the GPS velocities along a profile perpendicular to the fault.Subsequently, the obtained velocity field was utilized in the desired inversion process.This model facilitates the spatial migration of the creeping zone during the seismic cycle, offering valuable insights into the transitional region between the locked zone and the upper portion of the viscoelastic medium.
In this study, the optimal solutions exhibit a significantly low propagation velocity, measuring less than 1 m/year, as depicted in Figure 2. Furthermore, these solutions demonstrate minimal creep propagation.By incorporating these findings into the inversion process, the resulting slip rate on the fault plane aligns closely with the actual slip rate.The work conducted by [22] utilized fault paleo-seismological methods and presented a sliding rate ranging from 3.1 to 6.4 mm/year for the North Tabriz fault.In our current research, the slip rate falls within the range of 3 to 6 mm/year, establishing a strong agreement with the slip rate obtained through paleo-seismological methods.This underscores the superiority of our method over conventional approaches, both in terms of the scientific block model employed and the technical finite element method utilized.

Figure 1 .
Figure 1.GPS velocity vectors obtained from the forward problem using the b method.

Figure 1 .
Figure 1.GPS velocity vectors obtained from the forward problem using the boundary element method.
Eng. Proc.2023, 56, 26 7 of 8 and 15 km.Furthermore, the propagation velocity and creep propagation rate vary from 0 to 8.5 m/year.The relaxation time for the Best fitting solution is determined to be 156 years, with a range of 140 to 250 years.Figure2showcases the posterior distributions, representing the estimated parameters obtained through the creep diffusion inversion, along with the resulting propagation velocity.

Figure 2 .
Figure 2. Marginal posterior distributions for parameters estimated in the propagating creep inversion.Best-fitting solution are indicated by the red lines.

Table 2 .
The initial values required for the Monte Carlo Markov chain.

Table 3 .
The prior limits for the eight parameters required for the Monte Carlo Markov chain.