Uncovering a Seismogenic Fault in Southern Iran through Co-Seismic Deformation of the Mw 6.1 Doublet Earthquake of 14 November 2021

Abstract

On 14 November 2021, a doublet earthquake, each event of which had an Mw of 6.1, struck near Fin in the Simply Folded Belt (SFB) in southern Iran. The first quake occurred at 12:07:04 UTC, followed by a second one just a minute and a half later. The SFB is known for its blind thrust faults, typically not associated with surface ruptures. These earthquakes are usually linked to the middle and lower layers of the sedimentary cover. Identifying the faults that trigger earthquakes in the region remains a significant challenge and is subject to high uncertainty. This study aims to identify and determine the fault(s) that may have caused the doublet earthquake. To achieve this goal, we utilized the DInSAR method using Sentinel-1 to detect deformation, followed by finite-fault inversion and magnetic interpretation to determine the location, geometry, and slip distribution of the fault(s). Bayesian probabilistic joint inversion was used to model the earthquake sources and derive the geometric parameters of potential fault planes. The study presents two potential fault solutions—one dipping to the north and the other to the south. Both solutions showed no significant difference in strike and fault location, suggesting a single fault. Based on the results of the seismic inversion, it appears that a north-dipping fault with a strike, dip, and rake of 257°, 74°, and 77°, respectively, is more consistent with the geological setting of the area. The fault plane has a width of roughly 3.6 km, a length of 13.4 km, and a depth of 5.6 km. Our results revealed maximum displacements along the radar line of sight reaching values of up to −360 mm in the ascending orbit, indicating an unknown fault with horizontal displacements at the surface ranging from −144 to 170 mm and maximum vertical displacements between −204 and 415 mm. Aeromagnetic data for Iran were utilized with an average flight-line spacing of 7.5 km. The middle of the data observation period was considered to apply the RTP filter, and the DRTP method was used. We calculated the gradient of the residual anomaly in the N-S direction due to the direction of the existing faults and folds. The gradient map identified the fault and potential extension of the observed anomalies related to a fault with an ENE-WSW strike, which could extend to the ~ E-W. We suggest that earthquakes occur in the sedimentary cover of the SFB where subsurface faulting is involved, with Hormuz salt acting as an important barrier to rupture. The multidisciplinary approach used in this study, including InSAR and magnetic data, underscores the importance of accurate fault characterization. These findings provide valuable insights into the seismic hazard of the area.

1. Introduction

The Zagros orogenic belt in southwestern Iran is a prime example of a fold-and-thrust belt at the leading edge of the Arabian–Eurasian continental collision zone. Seismic activity in the Zagros Fold-and-Thrust Belt (ZFTB) is interpreted from two distinct perspectives. One viewpoint posits that the majority of moderate to large mainshocks are due to basement faulting (e.g., [1,2,3]). Conversely, an alternative viewpoint primarily attributes these occurrences to faults within the middle and lower layers of the sedimentary cover (e.g., [4,5]). Most earthquakes occur along major thrust faults, such as the High Zagros Fault (HZF), the Mountain Front Fault (MFF), or transverse strike–slip faults. The recent perspective reinforces the notion that the basement is rarely affected by moderate to large earthquakes [4,5,6].

The ZFTB is generally divided into two parallel zones: the High Zagros (HZ) and the Simply Folded Belt (SFB) (Figure 1). This division is based on differences in topography, geomorphology, exposed stratigraphy, and seismicity [7]. The seismic activity in the SFB is primarily caused by blind thrust faults. These are not often associated with surface ruptures but are known for causing destructive earthquakes [3,7,8,9,10,11]. On 14 November 2021, at 12:07:04 UTC, a strong earthquake of Mw 6.1 struck southern Iran within the SFB. One and a half minutes later, another earthquake of similar magnitude occurred. The events were recorded by the Global Centroid Moment Tensor (GCMT), formerly known as the Harvard CMT project [12,13] (Figure 1).

The earthquakes caused 1 death and 99 injuries. The affected area comprises 1 city and 137 villages, with a population of approximately 9121 households (36,484 people). The events destroyed buildings in the affected area, known as Fin, and caused damage to various facilities and infrastructure. The earthquakes were also felt in neighboring countries, including the cities of Abu Dhabi, Dubai, and Ras al-Khaimah in the UAE.

Several earthquakes were reported in the nearby study area, which were attributed to the main faults of the ZFTB [9,14]. Roustaei et al. [14] used Interferometric Synthetic Aperture Radar (InSAR) to analyze the co-seismic displacements of the 25 March 2006 Fin earthquakes with Mws of 5.7, 5.5, 5.2, 5.0, and 4.9. Their findings revealed significant differences between the ~ENE–WSW strike of the fault identified by InSAR and the ~E–W strike indicated by body-wave solutions. The study also provided detailed insights into the well-resolved upper and lower segments of the fault, located at depths of 5–6 km and 9–10 km, respectively.

In reviewing the literature specifically concerned with the 14 November 2021 earthquakes, Fathian et al. [15] used radar data from the Sentinel-1 and ALOS-2 satellites and calculated co-seismic interferograms using the Differential Interferometric Synthetic Aperture Radar (DInSAR) method, which were then resampled and subjected to non-linear inversion to estimate the rupture mechanism and geometric parameters of the uniform slip source. On this basis, they identified two possible fault planes striking in a west–east direction but dipping oppositely towards the north and the south and with a maximum slip of 2 m.

Golshadi et al. [16] conducted an InSAR study combined with an analysis of seismicity data, including foreshocks and aftershocks. These data were recorded by sparse, locally established seismographs of the Iranian Seismological Center (IRSC), the International Institute of Earthquake Engineering and Seismology of Iran (IIEES), and a UAE national seismic network station.

Toker et al. [17] used teleseismically recorded broadband P-velocity waveforms in their analysis of finite-fault source inversion using data from the IRSC and identified two faults as the cause of the events.

Golshadi et al. [16] and Toker et al. [17] attributed the events to two different faults. However, these faults are approximately 50 km apart and have different orientations. Rezapour and Jamalreyhani [18] found that the distribution of relocated aftershocks from the IRSC data indicates that the fault(s) associated with the Fin doublet earthquake extend along a nearly west–east striking, SSE-dipping thrust fault.

Significant strides have been made in identifying earthquake precursors, with notable examples including advances in land surface temperature (LST) and ionospheric scintillation. In a recent investigation, Boudriki Semlali et al. [19] demonstrated that anomalies in LST occurring a few days before earthquakes could effectively be precursors for earthquake events. Molina et al. [20,21] established a positive correlation between ionospheric scintillation and the period preceding earthquakes.

Ommi and Smirnov [22] used the region–time–length (RTL) parameter and the time-to-failure method to analyze the earthquakes that occurred on 14 November 2021. They identified a seismic quiescence anomaly ten months before the main shock. This anomaly developed gradually and reached its minimum value only one month before the main quake.

The ZFTB is known for an abundance of diapirs that formed from the Ediacaran to early Cambrian Hormuz salt. These salt diapirs are particularly widespread in the eastern part of the ZFTB and the Persian Gulf. However, the central region of the ZFTB displays salt diapirs distributed along tear fault systems trending nearly N-S and NE-SW and crossing the mountain belt obliquely [6]. Many researchers have documented the diapirs in detail (e.g., [6,23,24,25]).

Nissen et al. [26] investigated a cluster of three large earthquakes (which occurred on 27 November 2005, 28 June 2006, and 10 September 2008) and numerous smaller shocks at Qeshm Island within the SFB, located 50 km south of our research area. Their findings suggest that the Hormuz salt acts as an important regional barrier to rupture and that subsurface faulting contributes to the occurrence of the largest earthquakes in the SFB.

Nissen et al. [7] reinterpreted the seismicity of the ZFTB based on InSAR, teleseismic data, and geomorphology. Their analysis revealed that the seismogenic zone extends to a depth of about 20 km, with most earthquakes of moderate magnitude (Mw 5–6) occurring in the lower sedimentary cover at a 5–10 km depth and not in the basement. This agrees with the estimations of the sedimentary cover thickness in the Persian Gulf and the eastern part of the ZFTB, ranging between about 12 and 17 km [6,24,25]. Importantly, these faults often have no expression at the surface because they are decoupled from the surface folding by incompetent strata at different levels of the sedimentary cover.

Nissen et al. [4] and Elliott et al. [5] suggested that the crystalline basement undergoes predominantly aseismic shortening, contributing to a more sophisticated understanding of the seismic behavior and structural dynamics of the ZFTB.

Most of our understanding of active faults in the ZFTB is based on seismology, involving solutions for focal mechanisms and earthquake source parameters [11]. Nevertheless, there are still considerable uncertainties in detecting the seismogenic faults responsible for triggering earthquakes in the region. The lack of a dense local seismological network suggests that the infrastructure for data collection in the area was suboptimal for capturing the full extent and details of the foreshock and aftershock sequences.

The inconsistencies in the reports regarding the causes of the doublet earthquake led us to reconsider the search for its origins. To fill in the gaps and verify the results of the identified causative faults, the study performed an integrated analysis that included both InSAR and magnetic interpretation, which have not been used before.

This study used the DInSAR method with radar data from Sentinel-1. The next step involved the decomposition of the signals into vertical and horizontal displacements in the east–west direction using data from ascending and descending paths. This approach provided a better insight into the observed effects of earthquakes on the ground surface. Subsequently, seismic inversion was performed, which plays an important role in our study, as the inversion allowed the orientation of the fault to be determined. Both the amount of the fault slip and the affected area are related to the magnitude of the occurrence. They can be modeled using finite-fault inversion based on the displacements obtained from InSAR observations, Global Navigation Satellite System (GNSS) data, and/or recorded shock waveforms to reconstruct the slip history. In this study, we used inversion theory to infer seismic source parameters based on the indirect measurements mentioned above (InSAR).

The overarching scientific investigations that guided this research raised the following key questions: Which major fault(s) contributed to the occurrence of the doublet earthquake on 14 November 2021? Can the seismic activity be attributed to a single fault, or was it a complex interaction of several faults?

We hypothesized that a multidisciplinary investigation combining InSAR and magnetic surveys could help us identify the fault(s) responsible for the earthquake, which originated in the sedimentary cover closely linked to the main fault trends in the SFB. It is assumed that the earthquake did not originate from a point source but from a fault of a finite zone. Based on the elastic half-space model [27], we performed a non-linear inversion and estimated the most probable geometric parameters of the fault plane where the rupture caused the Fin doublet earthquake.

The results were compared to subsurface data using magnetic surveying as an effective tool to map tectonic structures, including faults.

Figure 1. The base map shows the three geological zones of Iran (Zagros, Central Iran, and Makran). The yellow rectangle shows the study area. The active faults are marked as HZF, MFF, and ZFF in the ZFTB [2]. The ZMP marks the transition between the Zagros collision and the Makran subduction zone [28]. The beachballs show the mechanisms of the faults that caused the earthquakes (the data are available at https://www.globalcmt.org/CMTsearch.html, accessed on 19 October 2023). The blue arrow shows the GPS velocities relative to the Eurasian plate near the Strait of Hormuz (KHAS) [28]. The blue square is Bandar Abbas city. The suture between the Arabian margin and central Iran is labeled MZRF and MRF. The identification of faults on the map is based on data from the Iranian faults map [29].

Figure 1. The base map shows the three geological zones of Iran (Zagros, Central Iran, and Makran). The yellow rectangle shows the study area. The active faults are marked as HZF, MFF, and ZFF in the ZFTB [2]. The ZMP marks the transition between the Zagros collision and the Makran subduction zone [28]. The beachballs show the mechanisms of the faults that caused the earthquakes (the data are available at https://www.globalcmt.org/CMTsearch.html, accessed on 19 October 2023). The blue arrow shows the GPS velocities relative to the Eurasian plate near the Strait of Hormuz (KHAS) [28]. The blue square is Bandar Abbas city. The suture between the Arabian margin and central Iran is labeled MZRF and MRF. The identification of faults on the map is based on data from the Iranian faults map [29].

2. Tectonic Setting

The study area lies within the ZFTB, which spans approximately 1800 km from northern Iraq to the Strait of Hormuz in Iran (Figure 1). The ZFTB is a significant part of the Alpine–Himalayan orogenic belt and the most seismically active fold-and-thrust belt in the world.

This region is an active continental collision zone that has promoted convergence between the Arabian and Eurasian plates for approximately 35 million years [30]. The ZFTB consists of sedimentary layers covering the entire Phanerozoic and has a thickness of 10–15 km [7,31,32,33,34], reaching up to 17 km in some parts of the Persian Gulf and the eastern part of the belt [6,24,25]. The stratigraphy, characterized by a combination of platform carbonates, evaporites, and marine shales, exerts a significant influence on the style of deformation [7]. The collision led to a shortening of the belt, and its prominent feature is the NW-SE double plunge folds at the surface, which are supported by steep reverse and thrust faults in the basement and NW-SE to N-S strike–slip faults [2,3]. This general trend deviates significantly in the easternmost part of the ZFTB, including the study area, where W-E and SW-NE folds with frequent fold-axis rotations are present.

The northeastern (and topographically highest) part of the ZFTB is referred to as the HZ, while the southwestern (and topographically lower) part is referred to as the SFB [14] (Figure 1). The HZ is bounded on its northeastern side by the Main Zagros Reverse Fault (MZRF) and on its southwestern side by the High Zagros Fault (HZF). The MZRF is commonly interpreted as a suture between the former rocks of the Arabian continental margin and the igneous and metamorphic rocks of central Iran [32]. While the MZRF is no longer active in reverse motion, the northern Zagros region of the same fault zone exhibits a right-lateral component of Arabian–Iranian convergence, known as the Main Recent Fault (MRF) [3,35], confirming a slip rate of 3 ± 2 mm per year [28]. The HZF is a large, NE-dipping thrust and marks the southern boundary of Paleozoic rocks in the Zagros at most locations [2]. The SFB extends from the HZF to the Persian Gulf. Most of the active deformation in the Zagros is concentrated in this zone and likely migrated southwestward from the HZF at an earlier stage [3,36,37,38]. A notable change in stratigraphy, often associated with elevation changes, occurs in certain folds within the SFB. These folds are closely related to prominent north-dipping basement faults, known as master blind thrusts, which play a crucial role in shaping structures in the region [2]. The most significant of these faults is the MFF (Figure 1), which marks the southwestern limit of surface exposures. The second major basement fault is the Zagros Foredeep Fault (ZFF) (Figure 1), which delineates the boundary of the Quaternary alluvium in the northern coastal plain of the Persian Gulf [7].

The study area, referred to as Bandar Abbas, is located in the SFB and serves as a transition zone between the Zagros and the Makran accretionary prism. The Zendan–Minab–Palami (ZMP) faults play a crucial role in adjusting the velocity gradient between the eastern ZFTB and western Makran areas (Figure 1) and delineate this transition zone [28].

3. Seismicity and Structural Geology

This section aims to enhance our understanding of seismic events and surface structures in the study area. The SFB is characterized by parallel double-plunge anticlines and synclines that dominate the short-wavelength topography and surface structure.

These features were initially categorized as detachment folds formed by buckling of the sedimentary cover above the Hormuz salt [39].

The study area displays parallel folds trending approximately east–west (Figure 2). The Hormuz salt has been exposed in the Handun, Namak, Guniz, and Anguru anticlines [40]. The original thickness of the Hormuz salt may have been several kilometers, although much of it was likely removed by diapirism and erosion [5].

The strongest earthquakes documented in the Bandar Abbas area of the SFB did not exceed a magnitude of Mw 6.7. These events are attributed to detachment levels comprising weak evaporitic and/or shale horizons [7]. The salt layers act as vertical separations within the overburden’s seismogenic layer, hindering seismic rupture propagation [7,41]. The Hormuz salt is identifiable at the surface by numerous diapirs, primarily located in the southeastern part of the SFB [42,43,44]. Subsurface salt diapirs and pillows have been discovered through the analysis of seismic profiles in the Persian Gulf [24,25,45].

According to GCMT data [12,13], the earthquakes occurred at depths of 17.2 and 13.8 km at coordinates 27.57°N, 55.98°E and 27.54°N, 56.03°E. Both earthquakes exhibited a reverse fault mechanism with a minor strike–slip component.

The surface faults in the study area were digitized based on the geologic map of Bandar Abbas at a scale of 1:250,000 [40] and the Iranian faults map [29] (Figure 3). An examination of geological maps and relevant reports on the faults responsible for the earthquakes [46] revealed no identifiable faults in proximity to the earthquake epicenters.

The MFF is considered the main fault near the epicenter of these earthquakes and is associated with numerous folds along its segments [2]. The origin and tectonic significance of the MFF system and its sinuosity remain unclear. Most hypotheses suggest strong structural control from geological inheritances [47]. Several earthquakes were reported near our study area, attributed to the main faults of the ZFTB [9,14].

To analyze and compare historical earthquakes regardless of any adjustments or relocation to the epicenters, we utilized the seismic catalog for the study area obtained from the United States Geological Survey (USGS). The dataset contains all earthquakes with magnitudes greater than four that occurred from 1900 until 13 November 2021, one day before the earthquakes of 14 November 2021 (Figure 3).

4. Materials and Methods

This section details the materials and methods used to identify the causative fault(s) of the events. Our approach integrates InSAR, seismic inversion, and magnetic studies.

4.1. InSAR

The InSAR technique played a crucial role in this study, offering a valuable alternative for visualizing fault surface projections, and is essential for understanding earthquake behavior. The analysis of interferograms primarily reveals surface deformation over time. We utilized Sentinel-1 SAR images acquired before and after the earthquakes (

Table 1. Sentinel-1 pairs were used to detect the deformations caused by the doublet earthquake on 14 November 2021.

No.	Satellite	Flight Direction	Path	Master’s Date	Slave’s Date
1	Sentinel-1A	Ascending	57	13 November 2021	25 November 2021
2	Sentinel-1A	Descending	166	9 November 2021	21 November 2021

) and calculated interferograms using the DInSAR method.

Interferogram calculations were performed using GMTSAR software (V6.2) [48], while phase unwrapping was executed with Snaphu software (v2.0.5) [48]. The Shuttle Radar Topography Mission (SRTM) version 1 [48] was employed to correct the wave phase relative to the ground surface (Figure 4).

The ground surface displacement determined using the DInSAR method reflects the displacement along the satellite’s line of sight (LoS).

Decomposition was used [49,50] to calculate two displacement components: the vertical component (d_V) and the horizontal component (d_EW) in the east–west direction.

$$\left[\begin{array}{c}{d}_{asc}\\ d_{dsc}\end{array}\right]=\left[\begin{array}{cc}-sin\left({\theta }_{asc}\right)cos\left({\theta }_{asc}\right)& cos\left({\theta }_{asc}\right)\\ -sin\left({\theta }_{dsc}\right)cos\left({\theta }_{dsc}\right)& cos\left({\theta }_{dsc}\right)\end{array}\right]\left[\begin{array}{c}{d}_{EW}\\ d_V\end{array}\right]$$

(1)

where $d_{asc}$ and $d_{dsc}$ are the LoS displacement vectors ($d_{LOS}$) from the descending and ascending paths, respectively, and $\theta$ is the incidence angle.

4.2. Seismic Inversion and Fault Modeling

We characterized the earthquake source by inferring the geometric parameters of possible fault planes through measurements and indirect observations made via InSAR, using the theory of seismic inversion. The probabilistic approach consisted of a combination of observations, theoretical data, and a priori information derived from the probability distributions of the earthquake. It was supported by the results of an a posteriori probability distribution. Building upon this foundation, we determined the most probable parameters of the earthquake source model, encompassing the offset, depth, length, and width of the fault plane, along with the slip, rake, dip, and direction of movement, which were consistent with the assumptions and observations.

Inversion of earthquake sources for a finite rectangular fault plane based on InSAR data was performed using Pyrocko scripts and Grond software (V1.6.1) employing the Bayesian bootstrap-based probabilistic joint inversion method. This technique estimates the distribution of estimation errors using multiple sampling with out-of-sample replacement [51]. The Bayesian method allowed us to determine the set of all likely source model parameters consistent with the observations and assumptions about the earthquake source and wave propagation.

The input data we used were unwrapped co-seismic interferograms obtained from the ascending and descending paths to determine the true course of the deformation field on the surface. Before the calculations, each of the interferograms was corrected for the influence of the troposphere with the use of the Generic Atmospheric Correction Model (GACOS) [52], which is based on the Iterative Tropospheric Decomposition Model (ITD). The data were divided into quadtrees to reduce the time required for calculations. In this process, four key parameters were determined: the variance threshold in each quadtree’s tile, the percentage of allowed NaN pixels per tile, and the minimum and maximum size of the tiles. In the case of the ascending path interferogram, the following values were finally used for the mentioned parameters: 0.007, 0.90, 0.01°, and 0.35°, while for the data from the descending path, the values were equal to 0.010, 0.90, 0.01°, and 0.32°. Based on this analysis, a variance–covariance matrix, subsequently utilized in the optimization process, was generated using the Kite software (V1.3.0).

For the optimization process, in addition to the InSAR data obtained for both the ascending and descending paths, along with additional parameters such as timing (date, centroid time, and assumed duration), location (depth, longitude, and latitude), and seismic parameters (magnitude, scalar moment, moment tensor, and focal mechanism) derived from the GCMT catalog, Green’s functions (precomputed for a regional CRUST2 Earth model) were used. Parameter estimation was conducted following data processing and the setting of search limits and the number of iterations.

4.3. Magnetic Study

To unveil the hidden fault(s) that caused the seismic events in the study area, we used aeromagnetic data for Iran within the longitude range of 55.60°–56.30°E and the latitude range of 27.20°–27.90°N. The data were surveyed by the Aeroservice Company of the USA under commission from the Geological Survey of Iran during 1974–1977, with an average flight-line spacing of 7.5 km [34,53].

To mitigate spurious frequencies that might arise during the grid generation process, we applied a low-pass filter with a cutoff wavelength of 15 km to the data. The survey was performed mostly at a constant barometric height of 2 km above mean sea level for the study area [53], which provides a suitable resolution for tectonic studies. The effects of the diurnal changes in the Earth’s magnetic field and the effects of the Earth’s core magnetic field had already been removed from the data [53]. Therefore, the data only needed to be corrected by reduction to the magnetic pole (RTP) filter.

RTP correction centers the asymmetrical magnetic anomalies over their sources, like the measured anomalies at the Earth’s magnetic poles, thus making magnetic interpretation easier [54,55]. To apply the RTP filter, we considered the middle of the data observation period. Additionally, we used the Differential Reduction to Magnetic Pole (DRTP) method [56,57] to assess the variation in the inclination and declination of Earth’s magnetic field. The reduced-to-pole magnetic anomaly of the study area is shown in Figure 5a.

There are several methods that can be used to enhance potential field data to identify the edges of buried structural contacts and lineaments. One of the simplest and most efficient methods uses directional derivatives, such as horizontal derivatives, which are well suited for locating the lateral contrasts of density and magnetic susceptibility of structures placed next to each other [56]. In this method, the gradient of a field anomaly is calculated considering the orientation of the structures. The gradient of the potential field is highest in the direction perpendicular to the structures, facilitating better detection of the edges of anomalous bodies. Therefore, this perpendicular direction is preferred for calculating the gradient of the field anomaly. Before performing any process, such as calculating the derivative, it is better to separate the regional and residual anomalies originating from deeper and shallower depths, respectively (Figure 5c,d). This process facilitates the detection of structures located at different depths. There are many methods that can be used to separate the anomalies, such as removing the polynomial trends [58], upward continuation [59], and matched filtering [60,61,62].

We used the matched filtering method, which involves spectral analysis of field data. The power spectrum of reduced-to-pole magnetic anomalies shows the average depths of 5 and 25 km for structures that create the anomalies (Figure 5b). We designed a matched filter to extract the residual anomalies corresponding to shallower structures.

5. Results

This section presents the results obtained via the methods described in Section 4.

5.1. InSAR

The decomposition of the data into horizontal and vertical displacement components offered a significant advantage: it enabled precise determination of the movement along the primary earthquake fault zone [63].

We calculated the vertical and horizontal displacements in the east–west direction by analyzing data from two paths. The primary displacements along the LoS (Figure 6) revealed values of up to −360 mm for the ascending orbit. This indicates an unidentified fault with horizontal surface displacements ranging from −144 to 170 mm and maximum vertical displacements between −204 and 415 mm (Figure 7).

5.2. Seismic Inversion and Fault Modeling

We modeled the earthquake source by inferring the geometric parameters of a potential fault plane through measurements and indirect observations, specifically using InSAR and the theory of seismic inversion.

Based on the calculations, we derived two geodetic solutions for the fault plane that best fitted the displacements obtained from the InSAR data. The first solution features a north-dipping, ENE–WSW-oriented fault plane with the following parameters: strike: 257°, dip: 74°, and rake: 77°. The fault plane dimensions are approximately 3.6 km in width, 13.4 km in length, and 5.6 km in depth (Figure 8).

The second solution represents a south-dipping, ENE–WSW-oriented fault plane with the following parameters: strike: 80°, dip: 13°, and rake: 89°. The fault plane’s dimensions are a depth of 4.0 km, a length of 13.4 km, and a width of 7.3 km (Figure 9).

A graph with a probability density function (PDF) is presented for each parameter in both of the obtained solutions, north-dipping (Figure 10) and south-dipping (Figure 11).

5.3. Magnetic Study

We calculated the gradient of the residual anomaly in the N-S direction, considering the orientation of the existing faults and folds. This analysis enabled us to identify the fault, as recognized by seismic inversion, through an examination of the horizontal gradient of the residual magnetic anomaly. The gradient map reveals a potential extension of the observed anomalies, attributed to a fault striking in an ENE-WSW direction, likely extending towards the ~ E-W(Figure 12). We determined the fault extension by simultaneously investigating the gradient map, the approximate locations of the main shocks, and the area’s seismicity map. In this analysis, we focused on interpreting the major fault related to the seismic events and avoided evaluating all existing anomalies in the gradient map. The black dashed line depicts the proposed fault extension (Figure 12).

6. Discussion

The study aimed to investigate the underlying causes of the doublet earthquake that occurred in Fin on 14 November 2021. Seismic inversion results indicate the presence of two potential fault solutions: one dipping towards the north and the other towards the south. Both suggest a single fault and show no significant variation in strike and location.

The examination of the horizontal gradient of the residual magnetic anomaly confirms the fault identified in the seismic inversion results. The analysis suggests the potential extension of observed anomalies resulting from a fault striking in an ENE-WSW direction, likely extending towards the ~E-W. However, it remains unclear whether the fault with the north or south dip aligns with the geological setting of the area.

To address this, we thoroughly analyzed the history of structural sections in the region. Notably, our study area was limited to a few published regional structural cross-sections (Figure 13e,f) explicitly highlighted by Molinaro et al. [64] and Jahani et al. [45].

Molinaro et al. [64] identified a low-dip thrust fault on the southern limb of the Handun anticline (Figure 13e). This fault poses challenges and uncertainties because the observation is based solely on surface data, which may not provide a comprehensive understanding.

The focal mechanisms of the events provide evidence for a thrust fault (Figure 13a). Furthermore, displacement distribution observed through the LoS indicates upward movement in the northern part of the proposed fault (Figure 13b,d), suggesting that the fault dips to the north rather than the south.

The strikes and dips for potential fault solutions in this study align closely with the results of Fathian et al. [15], who used radar data from the Sentinel-1 and ALOS-2 satellites to identify two fault planes. Both planes are oriented east–west and do not emphasize the source geometry. The first plane had a northward dip (71° dip, 263° strike), while the second plane had a southward dip (15° dip, 84° strike).

Previous studies, such as those conducted by Golshadi et al. [16], Toker et al. [17], and Rezapour and Jamalreyhani [41], predominantly relied on local seismic networks.

Golshadi et al. [16] reported that the first modeled plane corresponded to a NNE-dipping fault with parameter values of 284°, 53°, and 92° for the strike, dip, and rake, respectively. The second earthquake was associated with a SSE-dipping fault with parameter values of 85°, 6°, and 81°, with depths ranging between 7 and 14 km for the first event and between 1 and 7 km for the second. Golshadi et al. [16] endeavored to distinguish between a thrust and a back thrust in the Khurgu (Namak) and Handun anticlines. Toker et al. [17] presented the following results: the first fault plane was characterized by parameter values of 114°, 32°, and 91° for the strike, dip, and rake, with a SW dip, while the second fault plane had parameter values of 49°, 24°, and 27°, with a SE dip. They stated that the co-seismic displacement was confined to a depth range of 3 to 16 km. Rezapour and Jamalreyhani [41] suggested a more likely orientation of the fault plane with a west–east strike and a S-SE dip.

The investigations by Golshadi et al. [16] and Toker et al. [17] revealed significant discrepancies in fault locations attributable to differing interpretations and strike orientations of the causative faults. From our perspective, the identified faults, particularly the back thrust, are only approximations due to insufficiently reliable data. Further studies with robust references are needed to validate the authors’ structural interpretations. Additionally, Toker et al. [17] and Rezapour and Jamalreyhani [41] relied on local seismic stations, resulting in limited data coverage and inconsistencies.

Discrepancies in the attribution of faults, differing interpretations by researchers, and challenges in defining seismic source characteristics underscore the need for further research. To verify and enhance our understanding of the causes of earthquakes, we conducted the present study using a novel method that combines InSAR and magnetic data.

We propose a single causative fault with a northward dip, suggesting that two earthquakes occurred along it. This aligns with typical seismic events observed in the SFB, as highlighted by Nissen et al. [65], who indicated that most regional earthquakes with magnitudes around 6 (Mw) occur along reverse faults with dips of approximately 50°. Furthermore, these faults do not reach the surface but are buried at depths of approximately 4 km. Our seismic inversion confirmed that the fault plane depth is 5.6 km, consistent with results from magnetic interpretation (Figure 5b). This is particularly relevant, as earthquakes in the region are triggered within the sedimentary cover [4,5,6,65].

We refute the possibility of a south-dipping fault or back-thrust mechanism due to the lack of concrete evidence supporting upward movement in the southern segment. Confirmation of upward movement in the southern part of the block is necessary to consider such mechanisms.

The research findings indicate that the earthquake resulted from a single fault dipping northward, which extends approximately E-W, through the interpretation of magnetic anomalies.

Our study emphasizes the importance of a multidisciplinary approach to understanding the complexities of seismic hazard assessment in tectonically active regions like the ZFTB. Integrating different methodologies is crucial for developing effective risk mitigation strategies and fully understanding seismic hazards.

7. Conclusions

The current study employed a multidisciplinary approach to identify the causative fault and seismic sources of the doublet earthquake that occurred on 14 November 2021 in the Fin area. The study concluded that a single fault with a northward dip was responsible for both earthquakes.

The primary finding is that the fault has a strike, dip, and rake of 257°, 74°, and 77°, respectively. The fault plane measures approximately 3.6 km in width and 13.4 km in length. This previously unidentified fault showed horizontal surface displacements ranging from −144 to 170 mm and maximum vertical displacements between −204 and 415 mm during the earthquakes.

Complementing these findings, our magnetic study detected a concealed fault using aeromagnetic data and directional derivatives. The possible extension of observed anomalies was caused by a fault striking in an ENE-WSW direction, which likely extends towards the ~E-W, suggesting the presence of a major fault not previously identified.

The seismic inversion confirmed a fault plane depth of 5.6 km, consistent with the magnetic interpretation results. This is particularly relevant, as earthquakes in this region are triggered within the sedimentary cover.

While this study addressed some data gaps and uncertainties in previous research, future studies could incorporate structural fieldwork, structural cross-section analysis, non-seismic and microseismic investigations, and physical modeling.