Seismogenic Source Model of the 2019, M_w 5.9, East-Azerbaijan Earthquake (NW Iran) through the Inversion of Sentinel-1 DInSAR Measurements

Abstract

In this work, we investigate the M_w 5.9 earthquake occurred on 7 November 2019 in the East-Azerbaijan region, in northwestern Iran, which is inserted in the tectonic framework of the East-Azerbaijan Plateau, a complex mountain belt that contains internal major fold-and-thrust belts. We first analyze the Differential Synthetic Aperture Radar Interferometry (DInSAR) measurements obtained by processing the data collected by the Sentinel-1 constellation along ascending and descending orbits; then, we invert the achieved results through analytical modelling, in order to better constrain the geometry and characteristics of the seismogenic source. The retrieved fault model shows a rather shallow seismic structure, with a center depth at about 3 km, approximately NE–SW-striking and southeast-dipping, characterized by a left-lateral strike-slip fault mechanism (strike = 29.17°, dip = 79.29°, rake = −4.94°) and by a maximum slip of 0.80 m. By comparing the inferred fault with the already published geological structures, the retrieved solution reveals a minor fault not reported in the geological maps available in the open literature, whose kinematics is compatible with that of the surrounding structures, with the local and regional stress states and with the performed field observations. Moreover, by taking into account the surrounding geological structures reported in literature, we also use the retrieved fault model to calculate the Coulomb Failure Function at the nearby receiver faults. We show that this event may have encouraged, with a positive loading, the activation of the considered receiver faults. This is also confirmed by the distribution of the aftershocks that occurred near the considered surrounding structures. The analysis of the seismic events nucleated along the left-lateral strike-slip minor faults of the East-Azerbaijan Plateau, such as the one analyzed in this work, is essential to improve our knowledge on the seismic hazard estimation in northwestern Iran.

1. Introduction

On 7 November 2019 (22:47 UTC), a M_w 5.9 earthquake took place in the East-Azerbaijan (hereinafter referred to as E-Azerbaijan earthquake), in northwestern Iran, about 100 km east of Tabriz, the fourth largest city of Iran, with over two million citizens (Figure 1a). The event was clearly felt in most parts of the E-Azerbaijan province, killing at least five people, injuring hundreds, and causing widespread damage to the surrounding villages. During the next days, the mainshock was followed by more than 170 events with 2.5 ≤ M_w ≤ 4.8, recorded by the Broadband Iranian Network (BIN) [1]. The mainshock was recorded by different institutions: the Iranian Seismological Center (IRSC) [2], the United States Geological Survey (USGS) [3], the GFZ-Potsdam (GEOFON) [4] and the Harvard Global CMT catalog [5,6], that furnished four different hypocentral locations and focal mechanisms parameters. Despite such different solutions, the analysis of the available focal mechanisms and hypocenters suggests that the E-Azerbaijan earthquake nucleated along a strike-slip fault, located within the complex structural setting of the Bozgush Range [7,8] (Figure 1b,c). Although no medium-high intensity earthquakes have affected this seismogenic zone during the last years, an extensive active deformation and the consequent seismicity have interested the considered Iranian region and have represented a serious risk to the local cities and villages since ancient times. The historical seismic activity rate of this region is high, as testified by the many documented instrumental or historical earthquakes with magnitudes above M_w 6.0 [9]. In particular, one of the strongest historical seismic events took place on 22 March 1879, with M_w 6.7 along the South Bozgush Fault (black star in Figure 1b), causing 2000 casualties and heavy damage to the surrounding villages.

The seismicity of the E-Azerbaijan Plateau can be related to the oblique Arabia-Eurasia convergence, which is accommodated in northwestern Iran, partly as seismic activity along strike-slip and thrust faults [7,10,11,12,13,14,15,16,17] (Figure 1a). Accordingly, the E-Azerbaijan Plateau represents a complex mountain belt, which contains three internal major fold-and-thrust belts and, in particular, the Bozgush Range, bounded by two major faults, the North- and South-Bozgush Faults [7,8] (Figure 1b).

A wealth of data exists from historical seismicity and the related focal mechanisms [9,10,18,19,20], stress state reconstruction [7], field observations [21] and geodetic measurements around the involved seismogenic zone [14,16]. Thanks to the integration of these multidisciplinary data, it is possible to provide a detailed picture of the kinematics and the main geological structures of the region. Moreover, the fault geometry and slip distribution analysis for medium-large earthquakes are very significant for improving our knowledge of both the active faults distribution and the associated seismic hazard of a given region.

In this work, we first exploit the Differential Synthetic Aperture Radar Interferometry (DInSAR) measurements obtained by processing the data collected by the Sentinel-1 (S1) satellite of the Copernicus European Program along ascending and descending orbits. Then, we apply an analytical modelling approach to the computed coseismic DInSAR displacements, with the aim of better constraining the kinematics of the main seismic source. In addition, starting from the retrieved fault model characteristics and by considering the already known surrounding geological structures, we perform an analysis of the Coulomb stress transfer on the nearby faults, in order to investigate possible fault interaction processes.

2. Tectonic Setting

E-Azerbaijan, in northwestern Iran, is characterized by a complex active tectonics, linked to the main geodynamic regime controlled by the oblique convergence between the Arabian and Eurasian plates, with a present-day deformation rate of about 13–15 mm/year [7,10,11,14,16] (Figure 1a). In particular, this region represents a segment of the Alpine-Himalayan orogenic belt, whose origin is ascribable to several intense geodynamic processes (i.e., multiple accretion and continental collision consequent to the subduction and closure of the Tethys Ocean), which occurred between the late Eocene and the Oligocene epochs [7,18]. The local expression of this large-scale tectonic framework is represented by intense deformations and seismic activity, mainly located along strike-slip faults and thrusts, and analyzed through GPS measurements and earthquakes focal mechanisms (GCMT) catalog [12,16,23,24], respectively. Indeed, the recorded historical [25] and instrumental seismicity [19] in Iran suggests an intracontinental deformation concentrated in several mountain belts surrounding relatively aseismic blocks (i.e., central Iran, Lut and South Caspian blocks [14]).

Our study area is inserted in the tectonic context of the E-Azerbaijan Plateau, an articulated mountain belt containing three internal major fold-and-thrust belts (the Arasbaran, Ghoshe Dagh and Bozgush Ranges) [7], and in particular of the Bozgush Range, bounded by two major faults, the North- and South-Bozgush Faults (Figure 1b). Several interpretations about their kinematics have been proposed in literature: Zamani and Masson [7] suggest that these faults can be identified as major thrusts associated with the original paleogeographic configuration and build upon inherited, pre-existing structures in the late Cenozoic; on the other hand, Faridi et al. [8] propose a right-lateral component according to the stream deflections and the present day kinematic measurements [23]. However, all the authors suggest that the Bozgush Range can be interpreted as a positive flower structure, also known as a pop-up zone [26]. Moreover, recent structural investigations in the mountain discovered local range-parallel dextral faults and range-transverse conjugate sinistral faults, with a strong normal component in Neogene and Quaternary sediments. The eastern Bozgush Range terminates on the roughly N–S-striking Garmachay sinistral fault that, with the Shalgun–Yelimsi Fault, defines the eastern Bozgush Range as left-stepping step-over breached and fragmented by multiple faults. These many minor faults accommodate about a 30° counterclockwise rotation of the East Bozgush Mountain with respect to the West Bozgush Range and accommodated uplift of the Precambrian and Paleozoic rocks on the southeastern flanks of the range. In addition to the sub-parallel NW–SE right-lateral strike-slip major faults, which attracted much attention because they are responsible for historical and destructive earthquakes, the study area also includes sinistral NNE–SSW-striking faults, which, in spite of their recent activity, contribute to the accommodation of the Holocene tectonic regime [8] (Figure 1b,c).

3. DInSAR Measurements

In order to investigate the ground displacements associated with the considered seismic event, we exploited the Differential SAR Interferometry (DInSAR) technique [27], which allows the analysis of surface displacement phenomena, by providing a measurement of the ground deformation projection along the radar line-of-sight (LOS). The SAR data considered to retrieve the ground deformation associated with the occurred seismic event were acquired by the Sentinel-1 (S1) constellation of the Copernicus European Program. Benefiting from the short revisit time and the small spatial baseline separation of the S1 constellation, we generated several coseismic differential interferograms (Figure 2), with a spatial resolution of about 80 m. Among them we selected, for the seismic source modelling discussed in the following, those less affected by undesired phase artifacts (atmospheric phase delays, decorrelation noise, etc.), thus preserving good spatial coverage and interferometric coherence, which is very relevant in order to correctly retrieve the deformation pattern through the exploited phase unwrapping procedures [28,29]. In particular, the employed S1 data were acquired on 15 October and 20 November 2019, (Figure 2c,g), and on 16 October and 9 November 2019 (Figure 2e,h) along the ascending (ASC) and the descending (DESC) orbits, respectively (

Table 1. Main characteristics of the interferometric pairs considered for the seismic source modelling of the E-Azerbaijan earthquake.

Sensor	Dinsar Pair	Orbit	Perpendicular Baseline (m)	Track
Sentinel-1	15102019–20112019	ASC	149	174
Sentinel-1	16102019–09112019	DESC	−12.19	6

). On both interferograms, several fringes located near the epicentral area are clearly visible (Figure 2g,h); note that each fringe corresponds to a LOS-displacement of about 2.8 cm (i.e., half of the employed S1 C-band wavelength λ = 5.56 cm). Subsequently, starting from these interferograms, we generated their corresponding LOS displacement maps (Figure 3a and b) through an appropriate phase unwrapping operation [29]. In particular, the ascending track (Figure 3a) presents a deformation pattern characterized by negative values down to about −4 cm and positive values up to about 7 cm, indicating a sensor-to-target distance increase and decrease, respectively; moreover, the descending track (Figure 3b) shows a clearer deformation pattern, characterized by both negative LOS displacement values down to about −5 cm and positive LOS displacement values up to about 5 cm.

4. Modelling

4.1. Analytical Modelling

In order to retrieve the fault parameters, we jointly inverted the S1 DInSAR displacements acquired from ascending and descending orbits (see previous section), by performing a consolidated two-step approach that consists of a non-linear optimization to constrain the fault geometry assuming a uniform slip, followed by a linear inversion to retrieve the slip distribution on the fault plane (details about the adopted algorithms can be found in [30]).

The LOS displacements retrieved from the DInSAR interferograms were modelled with a finite dislocation fault in an elastic and homogeneous half-space [31], also applying a compensation for the topography (in order to take into account the real depth from the ground surface) and assessing possible offsets and ramps affecting the DInSAR measurements. Moreover, they were preliminarily sampled over a regular grid (240 m in all the considered area) to reduce the computation load.

Starting from a non-linear inversion algorithm based on the Levenberg–Marquardt (LM) least-squares approach [32], we searched for the source parameters and, thanks to multiple random restarts implemented within the LM approach, it was possible to catch the global minimum during the optimization process.

The results of the non-linear inversion are reported in Figure 3, where we show the LOS-projected original (Figure 3a,b) and modelled (Figure 3c,d) DInSAR measurements, as well as their residuals (Figure 3e,f), obtained as the difference between the original and the modelled data [30]. The retrieved results reveal that the modelled seismic source allowed to effectively replicate the DInSAR measurements and the best-fit solution consists of a left-lateral strike-slip fault. In particular, the retrieved fault plane is located approximately at 3000 (± 633) m depth b.s.l. (i.e., the depth of the center of the fault plane) and presents a strike of 29.17° (± 5), a dip of 79.29° (± 7), a rake of −4.94° (± 7), and a uniform slip of about 0.73 m (± 0.23), distributed over a source of about 6220 (± 1540) m × 5630 (± 2530) m (see

Table 2. Best-fit parameters of the seismic source retrieved from non-linear inversion modelling. The 1-σ uncertainty is also reported.

Parameter	Best-fit
Length (m)	6219.80 (±1540)
Width (m)	5629.90 (±2530)
Center Depth (m)	3096.60 (±633)
Dip (deg)	79.29 (±7)
Strike (deg)	29.17 (±5)
Rake (deg)	−4.94 (±7)
East (m)	720,969.80 (±667)
North (m)	4,176,230.30 (±737)
Slip (m)	0.73 (±0.23)

and Figure S2).

The second step of our modeling is represented by the linear inversion process with the computation of the non-uniform slip distribution, in order to have a more accurate estimate of the slip along the fault plane. In particular, the linear inversion was performed by using as starting model that obtained from the previously realized non-linear inversion (Table 2) and inverting the following system:

$$\left[\begin{array}{c}{\mathbf{d}}_{\mathbf{DInSAR}}\\ 0\end{array}\right]=\left[\begin{array}{c}\mathbf{G}\\ \mathbf{k}·{\nabla }^2\end{array}\right]·\mathbf{m}$$

where d_DInSAR is the DInSAR displacements vector, G is the Green’s matrix with the point-source functions, m is the vector of slip values, and ∇² is a smoothing Laplacian operator weighted by an empirical coefficient k to guarantee a reliable slip varying across the fault; further constraints were introduced by imposing non-negative slip values obtained via non-linear inversion [33]. The fault length and width were extended to reduce the border effects as much as possible; the fault plane was subdivided into patches of 0.4 × 0.4 km.

The final slip distribution over the modelled fault plane is shown in Figure 4a and reaches a maximum value of 0.8 m. The retrieved seismic source has a geodetic moment of 8.63 × 10¹⁷ Nm, corresponding to a moment magnitude of 5.92, coinciding with the value estimated by the IRSC (M_w 5.9) (Figure 4). Aside from a slight correlation between slip and width and slip and rake, the variance analysis indicates that the implemented inversion technique is able to properly solve all the parameters. The best model parameters (best-fit) only slightly differ from the mean values and this is evident from the 1D distributions, which resemble Gaussian distributions, and from the computed standard deviations (see Figure S2).

4.2. Coulomb Failure Function

Static stress changes play an important role in the occurrence of earthquakes, which tend to nucleate on faults that have experienced an increase in Coulomb stress by previous events, and typically reach values included in a range of 0.1–10 bars [34,35,36]. Since the early 1990s, several studies have shown that the space-time earthquake clustering can be explained by faults interaction [37,38]. Indeed, when an earthquake occurs, the state of the stress in the surrounding volume is altered as a consequence of the seismic dislocation and, if the parameters describing the seismogenic fault and the receiver faults are known, the change in the stress field can be computed through the standard elasticity theory. In particular, under the Coulomb failure hypothesis, it is possible to verify if failure on a receiver fault is promoted or not by the slip on the seismogenic source fault by using the Coulomb failure function (CFF) [38,39], which is defined as:

$$\mathrm{CFF}=\Delta \tau +{\mu }^{\prime }\Delta {\sigma }_n$$

where Δτ is the shear stress change, Δσ_n is the normal stress change and μ′ is the effective fault friction coefficient on the receive fault.

In this work, we use the fault model retrieved from the analytical modelling of the DInSAR measurements, discussed in Section 3, as source fault, and assume a uniform slip distribution to compute the Coulomb stress changes at a reference depth of 5 km that corresponds to the mean depth value of the aftershocks recorded by the BIN during the 7 November–24 December 2019 time interval. For this purpose, we use the software Coulomb 3.3 [40,41].

Starting from the geological structures reported in Faridi et al. [8] and according to the Coulomb 3.3 software convention [40,41], we set the receiver fault mechanism to strike = 95°, dip = 80° and rake 180° for the Fault 1 (F1); 190°, 80°, 0° for the Fault 2 (F2), known as the Shalgun–Yelimsi Fault; 275°, 80°, 180° for the Fault 3 (F3), known as the South Bozgush Fault; and 190°, 80°, 0° for the Fault 4 (F4). For each selected receiver fault mechanism, we tested two effective friction coefficients (i.e., 0.4 and 0.6) (Figure 5). We verified that this latter coefficient (i.e., μ′ = 0.6) only slightly modifies the Coulomb stress change map. Given the mechanism of the source fault (strike = 29°, dip = 80°, rake = −4.9°), the receiver faults having the larger portion falling in a positive CFF area are F1 and F3 (Figure 5a,b,e,f and Figure S3). However, this does not preclude the activation of the other two faults F2 and F4 (Figure 5c,d,g,h and Figure S3). Overall, our results indicate that the main event may have encouraged (i.e., positively stressed), with a positive loading, the activation of all the considered receiver faults. This is confirmed also by the distribution of the aftershocks (black dots in Figure 5) that occurred in proximity or exactly on the considered faults.

5. Discussion

Starting from the retrieved modelling results, we suggest that the seismic source responsible for the 7 November 2019, M_w 5.9, E-Azerbaijan earthquake consists of a predominantly left-lateral strike-slip fault (strike = 29.17°, dip = 79.29°, rake = −4.94°), with a center depth at about 3 km, located within the complex structural setting of the Bozgush Range [7,8]. By comparing our solution with the different focal mechanisms and hypocentral locations provided by the IRSC, the USGS, the GEOFON and the CMT, we remark that the geometric parameters of the seismogenic source presented in this study mostly agree with the fault plane solutions proposed by the IRSC. Our model indicates a slip distribution (retrieved maximum slip of 0.80 m), mainly located near the ground surface, which can explain the displacement field retrieved by the S1 DInSAR measurements; therefore, we can suggest that our solution can be considered compatible with a shallow source of the earthquake occurring in the complex tectonic setting of the Bozgush Range [7,8]. Furthermore, by analyzing the hypocentral distribution of the earthquakes nucleated after the mainshock and until 24 December 2019, we remark that the seismicity shows an evident southeast-dipping high angle alignment, which is in good agreement with our modelled fault plane (Figure 4). Moreover, the aftershocks distribution shows a good relationship with the CFF values computed along the surrounding structures [35] (Figure 5); indeed, also without accounting for the specific fault mechanism, about 35% of the aftershocks are located in positive CFF areas and, therefore, we can suggest that their origin can be attributable to the static stress transfer.

According to the retrieved results, and by comparing our modelled fault with the detailed structural framework reported in Faridi et al. [8], we suggest that our solution reveals a minor fault located west of the Shalgun–Yelimsi Fault and not mapped in the geological maps available in the open literature [7,8,21,26], whose kinematics is compatible with that of the surrounding structures and with the local and regional stress states (Figure 4). Moreover, Zamani and Masson [7] and Faridi et al. [8] have furnished a detailed reconstruction of the subsurface geology of the considered seismogenic area and have produced some geological sections of the examined region. Starting from these sections, we can suggest that the geometry and characteristics of the retrieved source are in good agreement with those of the represented adjacent faults.

The occurrence of the analyzed seismic event along a left-lateral structure is also confirmed by the fieldwork, the photogrammetry and the drone study, performed by a team of the Northwest branch and the Seismotectonics and Seismology Department of Geological Survey of Iran [21]. In addition, as presented in [21], the reported failures caused by the mainshock nucleation are located in proximity of our modelled structure and correspond mostly to some dynamic surface and slope instabilities (e.g., centimetric cracks, rockfalls, landslides, stone jumping features, and change in spring water-colors), which have been especially concentrated within the hangingwall of the causative fault zone. Moreover, we remark that the fractures generated by this event, as reported in [21], are secondary failures compatible with the kinematics of our modelled fault (Figure 4b). We further remark that, if we should consider as seismogenic source a structure consistent with the reported location and orientation of the Shalgun–Yelimsi Fault [8], the geodetic inversion of the exploited DInSAR measurements would result in a best-fit solution whose residuals are significantly worse than those achieved for our model (Figure S4). This is an additional confirmation of the validity of our findings.

The retrieved fault solution can be related with the reconstruction of the stress states performed by Zamani and Masson [7]; in particular, the complex active tectonics of the studied E-Azerbaijan region is characterized by several compressive structures, such as the Arasbaran, the Ghoshe Dagh and the Bozgush fold-and-thrust-belts, which have been generated by two distinct compressional stress systems, with NE–SW and NW–SE directions. The evolution of these compressive structures have also been associated with considerable Neogene to Quaternary volcanic activity and NE–SW and NW–SE strike-slip faults. The origin and kinematics of the local left-lateral strike-slip faults, such as the causative fault of the considered E-Azerbaijan earthquake, and the related seismicity, can be linked to the accommodation of the nearly N–S shortening between the Arabian and Eurasian plates, with the subsequent eastward extrusion of regional crustal scale blocks [42,43].

6. Conclusions

In this work, we have investigated the M_w 5.9 E-Azerbaijan (NW Iran) earthquake, in order to retrieve the causative fault of this seismic event. To this aim, we have exploited the available DInSAR measurements obtained by processing the SAR data collected by the Sentinel-1 constellation along ascending and descending orbits. Moreover, we have applied an analytical modelling two-step approach to the computed coseismic DInSAR displacements, in order to better constrain the kinematics of the main source.

Our main findings can be summarized as follows:

• The source model reveals a rather shallow seismic structure approximately NE–SW-striking and characterized by a left-lateral strike-slip, southeast-dipping faulting mechanism. The retrieved source reveals a minor fault not mapped in the geological maps available in the open literature, but it is characterized by a kinematics compatible with that of the surrounding structures, the local and regional stress states and with some of the field observations.
• Starting from the retrieved fault model characteristics and by considering the known surrounding geological structures, we have performed an analysis of the Coulomb stress transfer on the nearby faults, in order to investigate possible fault interaction processes. Our results indicate that the considered receiver faults may have been positively stressed by the main event and this is confirmed by the aftershocks distribution.
• The analysis of the seismic events nucleated along the left-lateral strike-slip minor faults of the East-Azerbaijan Plateau, such as the one analyzed in this work, is essential to improve our knowledge of the seismic hazard estimation in northwestern Iran.