STRESS TRANSFERE AND STATISTICAL ANALYSIS OF 2012 AHAR-VARZEGHAN SEISMIC SEQUNCE , NORTHWESTERN IRAN

In 11 of August 2012, two destructive earthquakes with Mw= 6.4 and 6.2 occurred between cities Ahar and Varzeghan (Northwest Iran). They had a close epicentral distance of 6 Km and also had a short time lag of 11 minutes. Following that, a high-rate of aftershock activity began where during the first month more than 2000 events (M≥0.7) affected several villages in the area. The seismic released energy induced significant damage and losses in an extensive zone. Right after the seismic doublet occurrence, a surface rupture with a primarily east-west orientation was observed. The idea of having an almost vertically dipped fault plane for the first shock is more consistent with the trace of the upper edge on the surface and the focal mechanism solutions which propose a steady dipping EW. Previous studies propose different geometries for the generating faults of the second earthquake. In this study, we associate the surface rupture with the first mainshock and both nodal plane explaining the relationship between the two main seismic events are discussed after Coulomb failure stress calculation due to the first shock. Then the stress transfer because of the doublet is analyzed in order to determine its consistency with the statistical modeling prediction for the aftershock population and spatial distribution. For statistical modeling a temporal version of Epidemic Type Aftershock Sequence (ETAS) is applied on one-year seismicity including events with minimum magnitude of 2.5.


Introduction
The northwestern Iran (Fig. 1a) is historically associated with destructive earthquakes mainly related to North Tabriz fault (NTF) which is the largest active fault in this area created due to the collision of the Arabia-Eurasia plateau.This fault stands at about 50 km to the southeast of Ahar and Varzeghan cities expanding along 150+ km in NW-SE direction.NTF in most parts has a nearly vertical dip and right lateral strike-slip character (Berberian and Arshadi 1976).During past two decades since the Iranian seismic network is covering this area, there has been no important seismic activity rather than dispersed seismicity on the northern region of NFT.GPS vectors at the NTF indicates a mainly north-northeast tectonic movement.Close to Ahar, the average slip (between 1999 and 2009) along the segment is of about 11 mm/year (Djamour et al. 2011).Right after the doublet a surface rupture was occurred in the area (Fig. 1b).It was mapped along 8 km by mean of an optical satellite images cross-correlation analysis (Copley et al. 2014).Field studies also revealed a longer ruptured distance up to 13.The east-west right-lateral strike slip solutions for both events' focal mechanism may lead to associate the rupture with the first, second or both shocks.At the same time, doublets with similar focal mechanisms are less expected to happen so close in both time and space (less than 10 km and 1 hour) and even less likely in an intraplate zone.In this work the spatio-temporal statistical characterization of the 2012 Ahar-Varzeghan seismic sequence is investigated.After this, we then analyze the static stress transfer because of the 2012 AV-doublet and study its correlation with spatial distribution of re-located series.

Data Analysis
Identifying the changes in the seismic behavior after the occurrence of the doublet where the seismic activity rate increases reaching up to 250 events per day needs a selection of data where they are homogeneosly registred in terms of time and magnitude.Here in order to select such cataloge we study the evolution of seimic rate and magnitude frequency since the AV-doublet occurrence.
Figure 2 shows that about three months later after the beginning of the series, a new jump in the activity rate (128 events per day) takes place following the largest aftershock with Mw=5.6 on 7th Nov 2012.In Figure 2, it is shown that it takes about one year for the seismic rate to cool down into a more stable state.Indeed, a oneyear time period might contain the essential information of seismic activity.This period starting from 1st Aug 2012 to the end of Jul 2013, includes 4558 events (M ≥0.7) inside the circle in Figure 1 (R=50 km) around the epicentral location of the AV-doublet.To assure an adequate analysis of the series we performed a completeness analysis of the catalog.The best minimum or cut-off magnitude Mc=1.9 is calculated using Wiemer and Wyss (2000) method and applying magnitude frequency law of Gutenberg and Richter (1956).Earthquake magnitudes are assumed to follow the Gutenberg-Richter law (GR) which describes the number of earthquakes with magnitude equal or greater than M as in Eq. (1): The b-value generally scatters around 1 but for rocks, however under high stress it might have a lower value (Scholz 2015;Schorlemmer et al. 2005).A b=0.75 is calculated for the whole catalog using the maximum likelihood estimation and Mc=1.9.However, the IRSC online database only gives the epicenter error of events with M -2.5.Thus, even though the obtained completeness magnitude is less than 2.5, the final selected catalog in this study will be 826 events with 1 -2.5 occurred during one year inside an epicentral area of about 2500 Km 2 and we named it AV1213.

Re-localization
Enhancing the original hypocenter of the earthquake is essential for further verifications of static stress changes due to the AV-doublet.Thus, AV1213 with 826 events was re-located using Double Difference re-localization method (Waldhauser and Ellsworth 2000).The initial condition for generating the event-pairs modifies the precision of the process and the percentage of relocalized events gets smaller under strongly confined coupling condition.A minimum requirement of having 8 double-difference equations for pairing earthquakes leads to 405 remaining events.The velocity model of Rezaeifar et al. (2016) for northwestern Iran with a mean Vp/Vs of order 1.76 is applied.Figure 3 shows how the 376 re-located events are concentrated in two main depth ranges; 4-9 and 15-22 km.

ETAS Modeling
The AV1213 sequence starts with two relatively big shocks.The upcoming earthquake population on one side consists of events triggered by earthquakeearthquake stress transfer.On the other side we also have events which might not be exclusively attributed to this internal stress loading.Instead, they are associated with external processes tend to change the state of seismic activity during upcoming months.The causing aseismic process in the case of northwestern Iran might be tectonic plate motions and stress buildup at north of Tabriz fault.A point process model with a conditional intensity function would characterize the seismic population in a way that at any time t, any previous earthquakes might be considered as a potential mainshock.In this section we use Epidemic Type Aftershock Sequence (ETAS) model (Ogata 1988) applying the algorithm developed by Marsan et al. (2013).Time variation of the background seismicity rate will be: 2"3$ & 4"3$ 5 6"3$ & 4"3$ 5 ∑   8 ( 9 :;< = ><?"3*3 = !?$ " 3 = #$ (2) Where t & and M & are the occurrence times and magnitudes of earthquakes.This formulation separates the time-dependent forcing (background) rate µ(t) and the earthquake rate υ(t) related to earthquakeearthquake triggering, where parameters c and p come from the Omori law (Utsu et al. 1995) and K ( and α are related to the magnitude-dependent aftershock productivity.The applied algorithm iteratively estimates the four parametersK ( , α , c, and p of the ETAS model by maximizing the log-likelihood value inside a time interval and then estimating the time-dependent background rate using the ± n-nearest neighbors.Figure 4 illustrate the temporal changes in the background seismicity rate for five obtained ETAS parameters in Table 1.From the descending shape of the decay it is evident that apart from the jump in the rate after the earthquake on November 7th, the background activity in late 2012 and early 2013 experience a smooth growth and finally decreases to ~0.05 (per day) or one event with M>2.5 per 100 days.

Assessing the source model for the doublet
First Shock: Concluding from the active faulting map that Ghods et al. (2015) provide together with the rupture alignments the faulting strike is divided into five vertical segments (Fig. 6).Such segmentation along almost 19 km might overestimate the change in the strike however it does not complicate the model but just build a more detailed geometry.For continental tectonic setting and a strike slip fault, the selected fault with 19.08 km surface rupture is empirically relevant with Mw=6.41 (Stirling et al. 2013) Figure 5: The cumulative number of events before and after ETAS modelling.In the previous studies and calculation by national and international agencies distinct depth for the first shock is determined.Donner et al. (2015) assume that the start of rupture is located at 14 km (Ghods et al. 2015) and the centroid they estimate using waveform inversion is located at shallower depth (6±1) km.Adopting their estimations, the constructed fault plane in this study have 14 km width and all main five segments start from surface and slop down 14 km (dip =90) as shown in Figures 7.
For weighting the net-slip we introduce a cushion shape model which implies assignation of the maximum weight for net-slip to the central area of the fault.The weighting function ., we use to make this cushion model, is the probability density function of a circle-arc with angle of θ (θ ≤180°).This function basically gives a simple curve for constructing weights along both length ."/$,and wide ."0$ of the fault.We did modifications on the circle-arc angle and offset of its center from depth of 7 km in order to simulate the measurements on the surface rupture.The resultant weighs for net-slip along depth ."0$ is shown in Figure 8 as gray belts for each vertical segments.In this figure, ."/$ is also shown for the net-slip along fault's strike; which is flattened in the middle section regarding to lake of evidence for centralized slip along specific longitude.
Considering the measured oblique slip at the west part, a gradual change in rake (increments of 5 ) is assumed; λ changes from 180 in eastern segment to 160 in western one.Finally, multiplying all weights by a constant value we simulate the measured displacements at the surface with maximum strike-slip of 73 cm and complete the construction of our source model (Fig. 9).Second mainshock; In contrary to the first mainshock, there is no common agreement about the orientation of the ruptured plane for the second mainshock.Previous studies argue both nodal planes.Donner et al. (2015) argue that the earthquake waveform data together with the the rock mechanism suggest that the slip on eastwest plane is less probable and the north-south plane with eastward dipping is more presumable as the causative fault for the second mainshock.Nevertheless, one possibility is that the second mainshock is an elastic rebound which responds to the abrupt change of Coulomb stress arrangements close to the first mainshock's fault.
Across each imaginary plane inside a block of rock the Coulomb failure stress ?@ or CFS, is defined as the difference between absolute value of tectonic shear stress on it and frictional stress that resist against rupture; The failure happens when ? @ exceeds a special positive value which depends on internal characteristics of the rock (King et al. 1994).Here in this calculate the resultant ∆CFS after the first mainshock for the second mainshock's CMT focal mechanism solutions.The calculations of ∆CFS for calibrated epicentral location (Ghods et al. 2015) at 38.425, 46.777 for both nodal planes is illustrated in Figure 10.Regarding to the calculated centroid and hypocentral depth for the second mainshock, we considered a depth range of 12±2 km for ∆CFS estimations.As it is evident the mean of ∆CFS for depth range 10 to 14 km shows how the north west of the first mainshock is significantly charged for left-lateral striking along north-south direction and uncharged for right-lateral striking along east-west direction.This result fortifies previous arguments that consider the northsouth oriented plane for the second mainshock.Hereafter, in order to introduce a source model for the second mainshock on the north-south CMT solution we consider and the hypocenter distribution relocated events (Fig. 3) we propose a 12x10 km ruptured area in accordance with Mw=6.2 (Wells and Coppersmith 1994).Such plane dips down with an angle of 50 to 19.66 km depth.In order to avoid any unknown complexity about the slip model and make a simple uniform model, six concentric planes are used to smooth the slip at fault's edges.This way maximum slip at the central area reaches to 105 cm which result a total seismic moment of 2.2 H 10 => Dyne.cm and Mw=6.2 (Aki, 1966, Stirling et al. 2013)

∆CFS on optimum failure directions
A stablished tectonic stress that address a block of rock at depth, produces different shear stresses along different orientations in the block.For a rock subjected to principal tectonic compressive stresses ?I JJJK, ?= JJJJK and ?L JJJJK where ?I M ?= M ?L , it can be shown that the maximum of CFS lays on a plane parallel the intermediate stress ?= JJJJK direction.This plane makes an angle of NO with the greatest principal stress ?I JJJK (King et al. 1994) where angle O directly depends of the effective coefficient of friction μ * .
Optimum faulting orientations expose the closest plane for further seismic ruptures.Although, in most cases we are dealing with stablished faults not future failures.These existing faults are planes of weakness and in another words their existence is changing the isotropic condition of the overburden blocks.However, Ghods et al. (2015) remarking the consistency of independently determined stress field, suggest that all the quaternary faults in the Ahar-Varzeghan area are potentially active and may remain aseismic for long time.Their results present the modern stress field of area with a maximum horizontal compressive stress ?I which on average is oriented along ∼ N128.5E.Let's consider that the smallest principal compressive stress ?L JJJJK is horizontal and ?= JJJJK vertical (?I JJJK and ?L JJJJK have no plunge).Therefor CFS gets larger on a vertical plane with strikes on ∼ 128.5°N 34.1° applying E * & 0.4 for pure strike-slip.
On the other hand the stress state calculated using fault kinematics data reveal a compressional stress arrangement in the area with ?L JJJJK vertical (Ghods et al. 2015).Thus the planes with maximum CFS might have strike along ?= JJJJK or ∼ 128.5°N 90° with dip of 29.5° (E * & 0.6) representing oblique-slip.
Figure 11 shows mean of ∆CFS along mentioned optimum orientations in depth ranges; 4-9 km and 16-21 km and because of the introduced source models for the AV-doublet in section 5.This figure is overlayed with a selection of relocalized events in both depth ranges.These selection is done based on a probability value (between 0 and 1) that ETAS modelling assigned to each events and indicates how likely that event belongs to the back ground population.Events with probabilities higher than 0.6 and less than 0.4 are slected for overlaying obtained ∆CFS maps.

Conclusions
The 2012 Ahar-Varzeghan double earthquakes standing among the largest events observed in this area destroyed many villages and provoked serious damages and losses (more than 300 deaths).Such catastrophic aspect together with the unknown ruptured faults and lake of seismological knowledge about the area, lately encouraged many semiologist and geologist to investigate this zone.In this paper the temporal evolution of the earthquake-triggering rate is estimated applying ETAS modeling.The result shows more than 55% contribution of aftershocks among over 820 events with M≥2.5.
The α-value=1.61which determine the magnitude dependence of the triggering mechanism is less than 2.3 which is expected for static stress triggering as the main earthquake generation process (Hainzl et al. 2010).This medium α-value together with b-value less than 1 are slightly indicating a higher level of confining stress that might be influencing the background seismicity in the area.
In this study we also tried to characterize a new source model mainly based on field measurements and geomorphological observations by Ghods et al. (2015).
The constructed source model with five 14 km vertically dipped segments makes a fault with 19 km length along east-west orientation (Fig. 7).We run ∆CFS calculation due to the first shock with an introduced cushion slip model.The result coincides with positive ∆CFS for northsouth nodal plane (CMT solution) for the second large shock at its hypocenter location given by Ghods et al. (2015).In contrary, we observe negative ∆CFS for eastwest strike-slip focal mechanism at this place.Accordingly, a perpendicular faulting geometry is suggested as the AV-doublet source.This double plane consists of a vertically dipped east-west plane along 19 km at shallow depth (0-14 km) and an eastward dipped north-south oriented plane along 12 km at higher depth (~12-20 km).
The associated ∆CFS because of the AV-doublet was calculated along optimum planes of ruptures respecting to the modern stress regime in the area.The positive ∆CFS areas are expected to be where the aftershock seismicity have significant contribution.Mainly because the internal elastic loading is addressing these weeker orientations (Stein et al. 1994).As it is illustrated in Figure 11 at higher depth aftershocks are pretty well placed on the charged zones wheras in the shallower depth the portion of aftershocks in discharged areas is not neglectable.Nevertheless, not all the aftershocks occur on optimally oriented ruptures and the probability of being from the background seismicity is taken from only one sample.Indeed, bootstapping is needed for having a more precise selection of events with aftershock attribution.
Mentioned facts together with limitation of our database in terms of the amount of provided phases and missed depth errors are directly affecting any spatial analysis.

Figure 2 :
Figure 2: Temporal variation of daily rate of seimicity during 26 months and magnitude dirstribution of the IRSC online catalogwith more than 5000 events with M≥0.7.

Figure 3 :
Figure 3: Result of hypoDD for over 300 relocated events projected to horizontal and vertical planes.

Figure 4 :
Figure 4: Time varying background seismicity rate obtained by ETAS modelling.

Figure 6 :
Figure 6: The ruptured area and location of mesernment.Yellow line indicates the strike of five segments we considered in this study.The red line shows the surface rupture trend.

Figure 7 :
Figure 7: The designed vertical fault plane and pdf of netslip on it for source model of the first shock.

Figure 8 :
Figure 8: The slip distribution along width and strike of the plane in figure 7.

Figure 9 :
Figure 9: The achived net-slip model for the first shock fault.The empirical formula byStirling et al. (2013)for stable continental tectonic regime results M 1 =6.41 for total

Figure 10 :
Figure 10: Changes in Coulomb failure stress due to the first shock for both east-west (up) and north-south (down) planes of the second large shock CMT focal mechanism solution.

Figure 11 :
Figure 11: The result of ∆CFS overalyed with events with probability of being of background population greater and equal to 0.6 (black squars) and less than and equal to 0.4 (yellow

Table 1 :
The ETAS parameter for which the Akaike.Information Criterion (AIC) yeilds n=21 as optimal smooting window