The Source Rupture Models and Seismogenic Structures of the Iran 2017 M_W 7.3 Earthquake

Abstract

The 12 November 2017 M_W 7.3 Iran earthquake was further studied. By analyzing Rayleigh-wave dispersion data, crustal models in the surrounding of the epicenter were obtained. It was found that there are high-velocity layers over a low-velocity zone. Using the obtained crustal models and a grid search procedure, the initial rupture depth of about 16.4 km and the rupture propagation velocity of about 1.62 km/s were retrieved. The source rupture models were established using the obtained rupture initial depth and the rupture velocity. The key features are as follows: The earthquake occurred on a shallow dip-angle fault, with ruptures spanning high-velocity layers in a depth range of approximately 7–25 km. A noteworthy observation from comparing crustal and rupture models is the existence of a low-velocity zone (layers) beneath the major rupture region (below about 25 km). It was realized that the seismogenic structure of this earthquake showed that high-velocity layers lay a low-velocity zone in the Zagros mountain seismic belt.

1. Introduction

On 12 November 2017, a significant earthquake with a moment magnitude (M_W) of 7.3 struck the Iran-Iraq border region (Figure 1), causing widespread devastation and loss. The seismic impact was felt as far away as the United Arab Emirates. Tragically, the earthquake resulted in the loss of over 600 lives and injured more than 7000 individuals. Moreover, approximately 70,000 people were rendered homeless, with around 12,000 homes destroyed and an additional 15,000 homes damaged.

Figure 1.

Background map depicting various elements related to the 2017 M_W 7.3 earthquake study. The beach balls show the projections of the focal mechanisms for four large earthquakes: The 10 October 2011 Turkey M_W 7.1; the 12 November 2017 Iran M_W 7.3; the 6 February 2023 Turkey M_W 7.8; and the 6 February 2023 Turkey M_W 7.5 earthquakes (the used data were retrieved from http://ds.iris.edu). The arrows with letter P indicate the direction of compressive force. The yellow dashed line running through the beach ball indicated with M7.3 illustrates the strike of the shallow-dipping nodal plane. The solid circles show the epicenters of earthquakes with magnitude ≥5, occurred between the 1 January 2000 and the 12 July 2023 (retrieved from http://ds.iris.edu). The long dashed lines show the 15 Rayleigh wave paths, used to obtain the crustal models in this article. The four triangles connected with the paths are the stations used to measure Rayleigh-wave dispersion data. The epicentral trend near and on the M7.3 beach ball shows the distribution of aftershocks that occurred within 1 month following the mainshock. The trend closely aligns with the strike direction. The 2017 M_W 7.3 earthquake occurred in the northwest of the Zagros mountain range [2]. This figure was created using the generic mapping tools [3].

The epicenter of this M_W 7.3 earthquake was located in the northwest part of the Zagros mountain range. The Zagros Mountains span southwestern Iran, extending in a northwest-southeast direction from the border areas of eastern Turkey and northern Iraq to the Strait of Hormuz. These mountains were formed through the ongoing collision between the Arabian plate and Central Iran, which began in the Miocene period [1].

This M_W 7.3 earthquake is one of the most significant instrumentally recorded earthquakes in the Zagros region, generating substantial interest among researchers. For instance, Barnhart et al. [4] employed interferometric synthetic aperture radar (InSAR) to analyze co- and post-seismic surface displacements, employing a fault inversion approach to determine permissible geometries, locations, and slips. Their findings indicated that the earthquake ruptured the shallow-dipping (14°–19°) Mountain Front Fault (MFF) within the Arabian crystalline basement at depths of approximately 12–22 km. Other researchers have also utilized InSAR techniques to study this M_W 7.3 earthquake (e.g., [2]).

The United States Geological Survey (USGS) has provided preliminary finite fault results for this earthquake. More than 10 papers have been published on or related to this earthquake (e.g., [2, 5, 6, 7, 8]).

However, given the varied strengths and weaknesses of existing findings, we were motivated to employ a waveform inversion method to retrieve more plausible source rupture models and surface wave dispersion data to retrieve crustal models in the epicentral region. Given the tectonic activity in the region, it is crucial to study this earthquake and the crustal structures in the epicentral region to enhance tectonic research and improve seismic hazard assessment.

In this article, we present the following: (1) introductions to the crustal modeling method, source rupture inversion method, and the grid search method; (2) the average crustal model obtained through the analysis of Rayleigh wave dispersion data in the epicentral region; (3) the initial rupture depths and rupture propagation velocities derived using a grid search procedure; (4) the source rupture models obtained by using shallow-dipping nodal planes and the identification of the active fault associated with the earthquake; and (5) the exploration of the seismogenic structure, which potentially reveals high-velocity crustal layers over a low-velocity zone in the Zagros mountain seismic zone.

2. Method introduction

In this article, three methods were used: crustal model retrieval using Rayleigh-wave dispersion data, source rupture model inversions, and a grip search procedure.

2.1 Crustal velocity modeling method

The crustal velocity model utilized in this study was obtained using a Rayleigh wave dispersion method. The method consists of two steps: (1) measuring the Rayleigh-wave group velocities and (2) modeling Rayleigh-wave dispersion data to extract S-wave velocities.

2.1.1 Measurement of Rayleigh-wave group velocities

Seismograms often exhibit prominent surface wave trains. These surface waves have different frequency components, and waves with different frequencies travel at different speeds, a phenomenon known as dispersion. As these surface waves propagate through the Earth, their amplitudes decay exponentially with depth. As a result, seismograms recorded at surface stations contain valuable information about the Earth’s structure. Rayleigh waves are generated by P and S waves, thus, S-wave velocity structures can be retrieved by analyzing Rayleigh-wave dispersion.

We measured the group velocities for the selected records using the multiple filter technique (MFT) developed by Dziewonski et al. [9]. A computer program developed by Herrmann and Ammon [10], was employed to implement the MFT method. In this technique, the group time for a specific frequency ${\omega }_n$ was determined as the time when the envelope of the filtered seismic record, as defined by Eq. (1), reaches its maximum. This time, regarded as the center of the Gaussian filter, is used to calculate the group velocity by dividing the epicentral distance by the group time.

$$h_n\left({\omega }_n,t,r\right)=\underset{-\propto }{\overset{\propto }{\int }}\left|F\left(\omega \right)\right|e^{-\alpha \left(\frac{\omega -{\omega }_n}{{\omega }_n}\right)2}cos\left[k\left(\omega \right)r-\mathit{\omega t}\right]\mathit{d\omega }$$E1

where $F\left(\omega \right)$ is the Fourier transform of the analyzed record, $e^{-\alpha \left(\frac{\omega -{\omega }_n}{{\omega }_n}\right)2}$ is a Gaussian window function, k(ω) is the wave number, and r is the epicentral distance.

The envelope of the filtering Rayleigh wave record can be computed using the formula proposed by Båth [11]:

$$g_n\left(t\right)=\sqrt{h_n^2\left({\omega }_n,t\right)+{\stackrel{-}{h}}_n^2\left({\omega }_n,t\right)}$$E2

where $\stackrel{-}{h}\left({\omega }_n,t\right)$ is the Hilbert transform of $h_n\left({\omega }_n,t\right)$ [9]. Herrmann and Ammon [10] introduced the formulas they used in their programs’ manual.

2.1.2 Modeling Rayleigh wave dispersion data for S-wave velocities

After obtaining Rayleigh-wave dispersion data at a particular seismic station, we can utilize these data to model the S-wave velocities along the path traversed by the Rayleigh waves.

The initial step involves establishing an initial crustal model and then the model is revised by comparing the observed dispersion with the predicted dispersion generated using the crustal model. The inversion is in fact to identify a revised model that best fits the observed data.

In our inversion process, we assumed that $\mathit{\delta U}\left({\omega }_k\right)$ was the group velocity difference between the observed and synthetic waveforms at frequency ${\omega }_k$, the crustal model had N horizontal layers on the top of a half-space, and the S-wave velocity is ${\beta }_n$ in the nth layer (n is the index for the layer). The difference at frequency ${\omega }_k$ (k = 1, …, M) could then be approximated using the following formula:

$$\mathit{\delta U}\left({\omega }_k\right)=\frac{\partial U\left({\omega }_k\right)}{\partial {\beta }_1}{\mathit{\delta \beta }}_1+\frac{\partial U\left({\omega }_k\right)}{\partial {\beta }_2}{\mathit{\delta \beta }}_2+\dots +\frac{\partial U\left({\omega }_k\right)}{\partial {\beta }_N}{\mathit{\delta \beta }}_N$$E3

The partial derivatives $\frac{\partial U\left({\omega }_k\right)}{\partial {\beta }_n}$ (n = 1, 2, …, N) were calculated using the formula proposed by Rodi et al. [12]:

$$\frac{\partial U\left({\omega }_k\right)}{\partial {\beta }_n}=\frac{U\left({\omega }_k\right)}{C\left({\omega }_k\right)}\left(2-\frac{U\left({\omega }_k\right)}{C\left({\omega }_k\right)}\right)\frac{\partial C\left({\omega }_k\right)}{\partial {\beta }_n}+{\omega }_k\frac{U^2\left({\omega }_k\right)}{C^2\left({\omega }_k\right)}\frac{\partial }{\partial \omega }\left(\frac{\partial C\left(\omega \right)}{\partial {\beta }_n}\right)$$E4

where $C\left({\omega }_k\right)$ is the phase velocity. $C\left({\omega }_k\right)$, $U\left({\omega }_k\right)$, and $\frac{\partial C\left({\omega }_k\right)}{\partial {\beta }_n}$ were obtained using standard Thomson-Haskell matrix calculations, and $\frac{\partial }{\partial \omega }\left(\frac{\partial C\left(\omega \right)}{\partial {\beta }_n}\right)$ were obtained by numerically differentiating $\frac{\partial C\left({\omega }_k\right)}{\partial {\beta }_n}$ (Rodi et al., [12]).

If we let m = (${\mathit{\delta \beta }}_1$, ${\mathit{\delta \beta }}_2$, …, ${\mathit{\delta \beta }}_N$), d = ($\mathit{\delta U}\left({\omega }_1\right)$, $\mathit{\delta U}\left({\omega }_2\right)$, …,$\mathit{\delta U}\left({\omega }_M\right)$, and $G_{k,n}$ = $\frac{\partial U\left({\omega }_k\right)}{\partial {\beta }_n}$ (k = 1, 2, …, M; n = 1, 2, …, N), we obtained a linear equation:

$$\mathit{D}{=}\mathit{Gm}$$E5

After solving this equation, we obtained a set of corrections to the S-wave velocities in the N layers (${\mathit{\delta \beta }}_1$,${\mathit{\delta \beta }}_2$, …, ${\mathit{\delta \beta }}_N$). We could form a new model by assigning the S-wave velocity in the nth layer, ${\beta }_{n\_\mathit{new}}={\beta }_n+{\mathit{\delta \beta }}_n$. The above steps were repeated if this new model could not generate a satisfied synthetic dispersion curve. The methods to solve Eq. (5) were extensively studied and used (e.g., [13, 14, 15]). We use the sur96 program in the package [10] to obtain our solutions.

As for the theoretical background and the technique of Rayleigh wave dispersion inversion, except for the above introduction, the pioneering work by Haskell [16], Dorman [17], and Dorman and Ewing [18] provide valuable references.

2.2 The source rupture modeling method

The procedure used to establish an earthquake rupture model is described below. One of the two nodal planes is used as the earthquake rupture plane, and a Cartesian coordinate system is set up on this plane. Usually, the x-axis is along the strike direction, and the y-axis is along the dipping direction. The selected rupture plane is divided into many small rectangles where the length of each small rectangle along the x-axis is dx and along the y-axis is dy. Each small rectangle is called a sub-fault and is treated as a point source. At a specific station, a synthetic seismogram is the summation of the synthetic seismograms generated by all rectangle sources. The source time function of all the sub-faults is usually depicted as overlapped triangles. The layout of the source rupture model can be found in Hartzell and Heaton [19]. A unit constant rupture slip vector for each sub-fault can be divided into two orthogonal vector components (one aligned along the strike and the other aligned along the dip direction). At a specific station, each component can be used to calculate a synthetic seismogram called Green’s function. Any slip vector on the sub-fault can be obtained by multiplying the two constant vector components with appropriate coefficients. An inversion aims to obtain the coefficients of all the sub-faults from the observed seismograms to establish the rupture model for the earthquake.

Once a nodal plane is selected as the rupture plane, it is divided into M × N sub-faults. The slip function (source time function) of each sub-fault is depicted by overlapped L triangles with a rise time τ, which is half the length of the bottom side of the triangle. The initial constant unit slip direction (slip0) of each sub-fault is obtained through the moment tensor solution. The initial slip is separated into two components in the directions of (slip0 + 45°) and (slip0 − 45°). This breakdown is convenient for Green’s function calculations.

Assuming that on a sub-fault mn (m = 1, ···, M; n = 1, ···, N) at the k^th component direction (k = 1, 2), the slip corresponding to the l^th triangle is X_mnlk, and the vertical components of the Green’s functions generated at Station j are g_mnkj, and then the vertical component synthetic seismogram W_j at the station generated by all sub-faults at the time point t_i can be expressed as follows:

$$W_j\left(t_i\right)=\sum _{\mathit{\text{mnlk}}}{X}_{\mathit{\text{mnlk}}}{g}_{\mathit{\text{mnkj}}}\left(t_i-\left(l-1\right)\tau -T_{\mathit{mn}}-{\mathit{dl}}_{\mathit{mn}}\right)$$E6

where T_mn is the rupture start time at the mn^th sub-fault, and dl_mn is a time delay generated by the different travel path lengths between the P waves generated at the mn^th sub-fault and the rupture start sub-fault m₀n₀ (hypocenter). The set of X_mnlk that can generate a W_j, which is most similar to the observed seismogram at Station j, is the best fitting rupture model for the earthquake.

The Eq. (6) was solved using the nonnegative least squares (NNLS) method [20]. To ensure stability in the slip solution, a smoothness constraint was imposed on the spatial distribution of the total slips using a Laplacian differential operator [21]). To calculate the time delay dl_mn in Eq. (6), a rupture velocity is required. Additionally, to compute Green’s functions, a reasonable initial focal depth and a crustal model are necessary.

2.3 Grid search method

A Grid Search is an optimization algorithm commonly used in various research fields. It is beneficial for selecting the best parameter values to optimize a model or solve a problem. The Grid Search helps identify the optimal parameter values that yield the best results by systematically evaluating the model’s performance across all possible parameter combinations.

Grid searches are particularly effective when the parameter space is not too large or when the relationship between the parameters and the model’s performance is poorly understood. However, they can be computationally expensive when dealing with many parameters or when the parameter space is extensive.

Let a model parameter vector V = (ν₁, ν₂, …, ν_m), and an objective function f = f(V). A simple way to set up a grid search consists in defining a vector of lower bounds a = (a₁, a₂, …, a_m) and a vector of upper bounds b = (b₁, b₂, …, b_m) for the vector V. The grid search involves taking n equally spaced points for each component. All the objective function values at each grid point are calculated and compared. Pick up the vector that corresponds to the minimum (or maximum, depending on how to define the objective function) of the objective function as the solution.

The problem with this method is that the number of evaluations increases exponentially as n and m increase. Because we cannot reduce m, decreasing n is the only possible way of assuring that the evaluation process stops in a reasonable time, but this decreases the validity of the solution. For specific situations, some vector components may need dense spacing, some not; as such, the spacing number n may differ for different vector components. In our case, V = (ν₁, ν₂); ν₁ is the source rupture initial depth, whereas ν₂ is the rupture travel velocity.

3. Data for the source rupture inversions and the Rayleigh-wave dispersion inversion

The data used for the source rupture inversions are P-wave segments in the teleseismic waveform vertical records. These records contain valuable information about the earthquake source. The waveform records were downloaded from the Incorporated Research Institutions for Seismology (IRIS) database. Specifically, records within an epicentral distance of 30°–90° surrounding the epicenter were selected for analysis. Only the high-quality P-wave segments were considered in the analysis to ensure high signal-to-noise ratios. The selected records underwent an instrument response correction to account for the recording instrument’s characteristics. Additionally, a band-pass filter with frequencies ranging from 0.01 to 0.1 Hz (equivalent to periods of 10–100 s) was applied to the records. The resulting dataset used for the source rupture inversion comprised 52 records.

The digital waveform records generated by moderate earthquakes around the M_W 7.3 main shock, registered at regional seismic stations, were first retrieved from IRIS. Then, the instrument responses were removed, and they converted the velocity records into displacement records using the SAC2000 software [22]. After conducting spectral analysis on the waveform records, it was found that the strong amplitudes were within periods of approximately 2–40 s. Therefore, a band-pass filter (1–50 s) was applied to the records. Once the displacement records were obtained, those that clearly displayed Rayleigh waves within the period range of interest were selected for the dispersion data measurements.

4. Crustal velocity models in the M_W 7.3 epicentral region

The P-wave segments generated by a crustal earthquake, recorded at a teleseismic station, contain the direct P wave radiated from the earthquake source, the reflected wave pP, and the S-wave converted P wave sP. To generate Green’s functions at teleseismic stations, an earth model is required. In our study, we formed an earth model by replacing the crustal part in the preliminary reference Earth model (PREM; [23]) with the crustal model obtained from Rayleigh-wave dispersion data.

At some stations, the Rayleigh waves generated by moderate earthquakes were strong enough to measure surface wave dispersion data. We retrieved waveform records from IRIS and selected 15 clean Rayleigh wave vertical records. Figure 1 shows the 15 Rayleigh wave travel paths, connecting the four stations and the nine epicenters of the selected moderate earthquakes. These travel paths pass through the epicenter region of the M_W 7.3 earthquake.

Once a group dispersion curve of a waveform record was measured at a given station, it is used to determine an average 1-D S-wave velocity model along the source station path. The P-wave velocity and density are obtained by converting the S-wave velocity model using a Poisson ratio (Vp/Vs = 1.732) and the Nafe-Drake relation [24].

Fifteen crustal velocity models were retrieved around the epicenter of the M_W 7.3 mainshock, and an averaged crustal model was obtained. The average model was simplified from 20 layers to 7 to expedite the calculation time for Green’s function. Figure 2 depicts the simplified crustal model. The model parameters are listed in Table 1. This model was utilized for the calculation of Green’s functions. For more details on the utilization of Rayleigh-wave dispersion data for retrieving crustal models, refer to the following studies: Dziewonski et al. [9], Herrmann and Ammon [10], and Motazedian and Ma [25]).

Figure 2.

This figure displays the average crustal velocity model obtained using Rayleigh-wave dispersion data from 15 Rayleigh wave travel paths. The two thick lines represent the final average crustal velocity model. The two thin lines represent the initial model. The Rayleigh wave travel paths used for the dispersion data are those shown in Figure 1.

∆H	Vp	Vs	density
8.0	5.3075	3.0656	2.5615
6.0	5.4695	3.1557	2.5941
3.0	5.6423	3.2551	2.6290
6.0	5.9186	3.4167	2.6861
3.0	5.7640	3.3283	2.6526
10.	5.4377	3.1405	2.5877
4.0	6.4516	3.5985	2.8356
0.0	6.8496	3.8205	2.9480

Table 1.

The simplified average crustal model.

In the table, column ∆H represents the layer thickness (km), Vp represents the P-wave velocity (km/s), Vs represents the S-wave velocity (km/s), and density (g/cm³). They are listed from top to bottom in Figure 2.

5. Source rupture models

In this section, we introduce the source rupture models obtained.

5.1 Selection of the rupture plane from the two nodal planes

Selecting a rupture plane from the two nodal planes is necessary when establishing a source rupture model. Typically, the distribution of aftershocks can be used to determine the rupture plane. In Figure 1, the smaller solid circles within and nearby the beach ball indicated by M7.3 represent the distribution of moderate aftershocks that occurred within 1 month after the main shock. The distribution trend aligns closely to the strike direction (351°) of nodal Plane 1 (the shallow-dipping plane). Thus, nodal Plane 1 was selected as the rupture plane.

5.2 The determination of the initial rupture depth and the rupture velocity

An initial rupture depth and a rupture propagation velocity are required to conduct the source rupture process inversion. A grid search scheme was employed to obtain reasonable values for these parameters. The depth range from 14.0 km to 21.0 km with an increment of 0.5 km was set, whereas the range for rupture velocity was from 1.40 km/s to 1.80 km/s with an increment of 0.05 km/s. At each grid point, the source rupture inversion was performed using the inversion code developed by Kikuchi and Kanamori, provided by Lingling Ye (personal communication). To speed up the calculations of the Green’s functions, we made some revisions to the code. The average variance from all the utilized records in each grid point (each inversion) was recorded. This variance represents the fit between the synthetic seismograms generated by the retrieved rupture model and the observed seismograms. A smaller variance value indicates a better fit. A total of 135 variance values were obtained from 135 grid points.

The contour map (Figure 3) illustrates the variance value changes with respect to the initial rupture depth and the rupture velocity. In the search procedure, the shallow-dipping nodal plane obtained by USGS (strike/351°, dip/16°, and rake/137°) was used. Two minima were found. One is at the rupture velocity of 1.47 km/s and an initial depth of 18.9 km (1.47, 18.9); its variance value is 0.1264. The second minimum is at (1.61, 16.4); its value is 0.1269.

Figure 3.

The contour map of variance values. At a specific point on the map, the color of a pixel shows the variance value at that point. A deeper color shows smaller variance value. The color bar at the right side shows the variance values represented by colors. A variance value shows the fit between the observed and the synthetic seismograms. The smaller the value, the better the fit. This map depicted the changes in variance values with the rupture’s initial depth and rupture propagation velocity. When the shallow nodal plane provided by the United States geological survey (USGS) was used as the rupture plane (strike/351°, dip/16°, rake/137°) to perform the grid search procedure, two minima were found (two small circles). One is at the rupture velocity of 1.47 km/s and the initial depth of 18.9 km (1.47, 18.9); its variance value is 0.1264. The second minimum is at (1.61, 16.4); its value is 0.1269. The second minimum is our preferred solution.

5.3 Retrieving rupture models using the obtained minima

With the initial rupture depths and rupture velocities at the two minima, we retrieved two rupture models for the M_W 7.3 earthquake. Figure 4 displays the comparison between the rupture distributions obtained using the two minima. The upper panel shows that there are four rupture patches. The phenomenon that one earthquake has four separated patches may not easily occur naturally. The bottom panel shows that two patches combined to form a larger prolate patch, indicated by letter A. This shape is similar to the shapes obtained by Nissen et al. ([2]; their Figure 5). For the same earthquake, the rupture models obtained using different methods should be similar. From this common sense, we prefer to take the bottom panel of Figure 4 as our solution.

Figure 4.

The comparison between the rupture distributions obtained using different minima. (a) the rupture distribution obtained using the USGS shallow-dipping nodal plane (strike/351°, dip/16°, rake/137°) at the minima (rupture velocity of 1.47 km/s, initial depth of 18.9 km). The star (*) marked with SP denotes the rupture start point. Four red areas show four rupture patches. The dashed circles represent the rupture fronts at specific time intervals, indicated as 10 (s), 20 (s), and so on. (b) the rupture distribution using the same nodal plane at the minima (rupture velocity of 1.61 km/s, initial depth of 16.4 km). One larger rupture area, indicated by letter A, was obtained.

Figure 5.

The comparison between the first 24 observed and synthetic seismograms. For each pair the upper trace (solid line) represents the recorded waveform, whereas the lower trace (dashed line) represents the synthetic waveform generated using the slip distribution depicted in Figure 4(b). Both the observed and synthetic waveforms were filtered with a band-pass filter, ranging from 0.01 to 0.1 Hz. The symbols and numbers on the left side of each pair indicate the station name, the P-wave vertical component, the station azimuth in degrees, and the ratio between the observed and synthetic maximum amplitudes. Among all the 52 waveform records used in the inversion, the poorest ratios were observed at stations SACV (0.64) and BBSR (1.32). They show the most discrepancies between the observed and synthetic waveforms. To save page space, the remaining 28 records were not shown up. The fit quality is similar to those displayed.

Figure 5 compares seismograms at 24 stations, generated using the rupture model in Figure 4(b). The fit between the observed (upper) and synthetic (bottom) traces is generally good. For brevity, the comparison for the remaining 28 stations is not provided; the waveform fit at those stations is similar to those shown in Figure 5. This figure allows for a visual comparison between the observed and synthetic waveforms, highlighting the agreement or discrepancies between the recorded and the simulated waveforms based on the obtained slip distribution. It provides, to some extent, insights into the modeling approach’s accuracy and the simulated waveforms’ fidelity in capturing the actual earthquake’s characteristics.

5.4 Rupture models retrieved using nodal planes obtained by other authors

We also utilized shallow-dipping nodal planes obtained by other authors. First, we show the results in detail using the nodal plane provided by Nissen et al. [2] as the rupture plane; then, list the key results using the nodal planes provided by other authors. The contour map in Figure 6 shows variance value changes obtained using the nodal plane (strike/353.7°, dip/14.3°, rake/136.8°) provided by Nissen et al. [2]. Three minima were found. In the deeper colored region, the rupture velocity range is about 1.45–1.75; the initial depth range is about 15–19 km.

Figure 6.

The contour map of variance values. When the nodal plane (strike/353.7°, dip/14.3°, rake/136.8°) provided by Nissen et al. [2] was used as the rupture plane, three minima were found (the small circles). One minimum is at the rupture velocity 1.62 km/s and the initial depth 18.2 km (1.62, 18.2); its variance value is 0.1199. The second is at (1.61, 16.1); its value is 0.1194. And the third is at (1.72, 16.0); its value is 0.1211. The second minimum is our preferred solution.

Three rupture models were generated using the above three minima (Figure 7). The upper panel shows the rupture distribution obtained using the nodal plane (strike/353.7°, dip/14.3°, rake/136.8°; Nissen et al. [2]) at the minimum (rupture velocity 1.62 km/s, initial depth 18.2 km). Several rupture patches were obtained. The middle panel shows the rupture distribution using the same nodal plane at the minimum (1.61 km/s, 16.1 km). One larger prolate rupture area, indicated by letter A, was obtained. The bottom panel shows the rupture distribution using the same nodal plane at the minimum (1.60 km/s, 17.2 km). One larger, oval rupture area, indicated by letter A, was obtained. The rupture distributions retrieved using different methods for the same earthquake should be similar. Because the rupture distribution in Panel (b) was similar to those obtained by Nissen et al. [2], so, we preferred to select Figure 7(b).

Figure 7.

The comparison between the rupture distributions obtained using different minima. (a) the rupture distribution obtained using the nodal plane (strike/353.7°, dip/14.3°, rake/136.8°) provided by Nissen et al. [2] at the minimum (rupture velocity of 1.62 km/s, initial depth of 18.2 km). The star (*) marked with SP denotes the rupture start point. Several rupture patches were obtained. (b) the rupture distribution using the same nodal plane at the minimum (1.61 km/s, 16.1 km). One larger prolate rupture area, indicated by letter A, was obtained. (c) the rupture distribution using the same nodal plane at the minimum (1.60 km/s, 17.2 km). One larger, oval rupture area, indicated by letter A, was obtained. The rupture distribution (b) was our preferred solution.

Other authors have also retrieved rupture planes for the 2017 M_W 7.3 Iran earthquake. The key results obtained using their nodal planes for preferred minima are listed in Table 2.

Authors	data	Strike; (°)	Dip; (°)	Rake; (°)	i.d.	Vr	va	m. s.; depth
USGS	L. P.	351.0	16.0	137.0	16.4	1.61	0.1269	14
GCMT	L. P.	351.0	11.0	140.0	16.4	1.59	0.1153	16
Nissen et al. [2]	B. W.	354.0	17.0	142.0	16.3	1.65	0.1268	13
Nissen et al. [2]	InSAR	353.7	14.3	136.8	16.1	1.61	0.1194	14
Barhart et al. [4]	InSAR	351.0	15.0	128.0	16.6	1.65	0.1258	14
Vajedian et al. [5]	InSAR	354.4	17.5	141.5	16.5	1.65	0.1278	14
Feng et al. [6]	InSAR	353.5	14.5	135.6	16.4	1.60	0.1210	14
Ding et al. [7]	InSAR	354.7	16.3	137.3	16.1	1.63	0.1249	14
Chen et al. [8]	InSAR	351.0	15.0	135.0	16.5	1.62	0.1235	14
Average		352.7	15.2	137.0	16.4	1.62	0.1235	14
GCMT	St. NP	121.	83.	82.	15.4	1.62	0.1535	18

Table 2.

The rupture plane parameters by different authors and the key rupture parameters.

In the table, the initial rupture depths (km, i.d.) and the rupture velocities Vr (km/s) were measured on a computer screen. The major slip (m. s.) depth (km) means the centre depth of the major rupture patch, visually estimated. L.P. means long period; B. W. means body waves; St. NP means Steep nodal plane. The other parameters are understandable.

6. A rupture model retrieved using a steep nodal plane as the rupture plane

All the previous models were obtained using shallow-dipping nodal planes as the rupture plane. This section presents the modeling results obtained using the same procedure but with a steep-dipping nodal plane.

Figure 8 displays the contour map depicting the variance values. In the search procedure, the steep-dipping nodal plane by The Global Centroid Moment Tensor (GCMT) (strike/121°, dip/83°, and rake/82°) was used. The small circle represents the obtained minimum variance value (0.1535), occurred at a depth of approximately 15.4 km and a rupture velocity of about 1.62 km/s. Figure 9 illustrates the slip distribution on the steep-dipping nodal plane. Several rupture patches were scattered within a 5–70 km depth range. Figure 10 compares the first 20 observed and synthetic seismograms generated using the rupture model shown in Figure 9. The poorest ratio occurred at station CMLA, with a ratio of 0.52.

Figure 8.

The contour map of variance values. When the steep-dipping nodal plane (strike/121°, dip/83°, rake/82°) provided by GCMT was used as the rupture plane, one minimum variance at the rupture velocity was 1.62 km/s and the initial depth of 15.4 km was found; its variance value is 0.1535.

Figure 9.

The distribution of the slip on the GCMT steep dip nodal plane (strike/121°, dip/83°, rake/82°). The star (*) marked with SP indicates the rupture start point at depth of 15.4 km. Multiple patches along the rupture plane were scattered at depths ranging from approximately 5–70 km. The rupture velocity of 1.62 km/s was used in the inversion.

Figure 10.

The comparison between the first 20 observed and synthetic seismograms. For each pair the upper trace represents the recorded waveform, whereas the lower trace represents the synthetic waveform generated using the slip distribution depicted in Figure 9. Among all the 52 waveform records used in the inversion, the poorest ratios were observed at stations CMLA (0.52) and ADK (1.23). To save page space, the remaining 32 records were not shown up; the fit quality is similar to those displayed. The variance value is 0.1535.

This test demonstrates that the minimum variance value (0.1535) obtained using a steep-dipping nodal plane was notably larger than those obtained using shallow-dipping planes (Table 2). One major rupture patch was as deep as more than 60 km (Figure 9), which is beneath the crust. It is impossible for this M_W 7.3 crustal event to have a major rupture patch beneath the crust. Therefore, the steep-dipping nodal plane could be ruled out to be the rupture plane.

7. The seismogenic structure of the 2017 M_W 7.3 Iran earthquake

The 2017 M_W 7.3 Iran earthquake occurred along the boundary between the Arabian and Eurasian plates in the Zagros mountain seismic zone (e.g., Nissen et al., [2]). The compressive force in the epicentral region is oriented northeastward (Figure 1), approximately perpendicular to the trend of the Zagros Mountains. In the epicentral region, a low-velocity zone exists from a depth of about 25 km to approximately 36 km beneath a high-velocity zone (Figure 2). The materials in the low-velocity zone are considered to be ‘soft,’ whereas stress accumulates in the ‘hard’ layers (high-velocity zone) at depths of about 7–25 km. When the stress strength in the ‘hard’ layers becomes strong enough, the crust in that region undergoes breakage, resulting in the M_W 7.3 earthquake. The source rupture modeling presented in Figure 4(b) and Figure 7(b) indicates that the ruptures occurred within a depth range of about 7–25 km, consistent with the depth range of the ‘hard’ layers. This suggests that the seismogenic structure of the 2017 M_W 7.3 earthquake possibly involves the high-velocity layers overlying the low-velocity layers along a significant seismic activity belt.

Overall, the modeling results introduced above indicate that the seismogenic structure of the 2017 M_W 7.3 earthquake likely involves the interaction between high-velocity layers and low-velocity layers along a seismic activity belt. The rupture models obtained using different shallow-dipping nodal planes consistently showed rupture patches distributed within the depth range of the high-velocity layers, reinforcing the correlation between the seismic activity and the properties of the subsurface layers.

8. Conclusions

The source rupture models of the 2017 M_W 7.3 Iran earthquake were established, and the crustal model in the source region was constructed. The ruptures obtained using different shallow-dipping nodal planes provided by different authors showed similar trends. The optimal initial rupture depths were about 16.4 km, and the rupture velocities were around 1.62 km/s. The rupture models retrieved using those shallow-dipping nodal planes exhibited rupture patches distributed within a depth range of about 7–25 km.

Unlike the common practice that the rupture’s initial depth and the rupture velocity are assumed or by trade-and-errors when the rupture modeling is performed, we used the grid search method to determine the optimal initial depth and rupture velocity. For each nodal plane used as the rupture plane, two or three minima formed by the initial depth and the rupture velocity were found. One minimum was selected by analyzing the shapes of ruptures and comparing the solutions obtained by other authors for the same earthquake. If solutions from other authors for the same earthquake are not available at the moment, we do not have an effective way to select the best minima.

To validate the selection of the shallow-dipping nodal plane as the rupture plane, we tested a steep-dipping nodal plane as the rupture plane. The test revealed that the minimum variance (0.1535) was obviously larger than that (∼ 0.1235) achieved using the shallow-dipping nodal planes. The worst maximum amplitude ratio between observed and synthetic data was 0.52 at station CMLA, a significant deviation from the ideal ratio of 1.0. In contrast, when the shallow nodal plane was employed as the rupture plane, the worst maximum amplitude ratio was 0.64. The rupture distribution associated with the steep-dipping nodal plane also exhibited a scattered pattern ranging from approximately 5–60 km in depth. As a crustal event, the ruptures less likely occurred beneath the crust. Based on the above discussion, it can be concluded that the selection of the shallow nodal plane as the rupture plane was justified. This test also indicates that modeling two nodal planes for the 2017 M_W 7.3 earthquake facilitated the identification of the most likely rupture plane. However, it is important to note that for other earthquakes, the rupture plane may not be discernible through modeling two nodal planes, as demonstrated by Motazedian and Ma [26].

Barnhart et al. [4] used InSAR observations of both the co-and post-seismic displacement to image the M_W 7.3 earthquake sequence. They implemented an iterative fault-slip inversion approach designed to elucidate the geometry and location of faults associated with this sequence. They found that InSAR observations are best described by co-seismic slip on a shallow-dipping (15–19°) plane that dips eastward at depths of 12–22 km. The rupture plane we identified was the shallow-dipping plane. The depth range of the major rupture we obtained (Figure 4b) was consistent with that obtained by Barnhart et al., [4].

The compressive force in the epicentral region is in the northeast direction. The epicentral region has a low-velocity zone under a high-velocity zone. This feature was displayed by the crustal models retrieved using the Rayleigh-wave dispersion data (Figure 2). The stress is accumulated in the ‘hard’ layers (high-velocity layers) above the ‘soft’ zone (low-velocity layers). When the accumulated stress is strong enough in the ‘hard’ layers, the crust there had to be broken, leading to the occurrence of this M_W 7.3 earthquake. Based on results obtained using several shallow-dipping nodal planes as the rupture planes, we conclude that the 2017 M_W 7.3 Iran earthquake occurred on a shallow-dipping fault at a depth range from about 7–25 km in the high-velocity layers over a low-velocity zone; the initial rupture depth was about 16.4 km, and the rupture propagation velocities were about 1.62 km/s. The central depths of the major ruptures were about 14 km.

The features exhibited by Figure 2 were also supported by the Love-wave dispersion inversion. The related article will be published. The factors leading to the uncertainty in the crustal model are multiple. One of them is the uncertainty in the measured dispersion data. It turned out that this uncertainty did not change the crustal features retrieved using the Rayleigh-wave dispersion data [27].

The epicenter of the 2017 M_W 7.3 earthquake is in the Lurestan arc region (e.g. Nissen et al., [2]). In the Lurestan arc region, the depth to the magnetic basement from spectral analysis of aeromagnetic data is about 16 km [28]. The initial rupture depth we retrieved is about 16.4. This may imply that the earthquake nucleated at the bottom of the basement.

Viewing this M_W 7.3 earthquake in a larger background, it can be seen that it is located in the boundary region between the Arabian Plate and the Eurasian Plate (Figure 1). Other large earthquakes also occurred in this corner region: 2011 M_W 7.1; 2023 M_W 7.8, and M_W 7.5 Turkey earthquakes. The time intervals between the 2011 Turkey M_W 7.1 and 2017 Iran M_W 7.3, 2017 Iran M_W 7.3 and 2023 Turkey M_W 7.8 and M_W 7.5 are about 6 years. The compressive forces were at a similar direction for the M_W 7.1 and M_W 7.8; they were at another similar direction for the M_W 7.5 and M_W 7.3. The locations, intervals, and forces of large earthquakes are another topic. Here, we only mention the above phenomena.

Nissen et al. [2] estimated that the rupture velocity of this 2017 M_W 7.3 earthquake is about 1.5–2.0 km/s, whereas the rupture velocity we obtained is about 1.62 km/s. These velocity values are at the lower end of the range of rupture velocities observed in large earthquakes globally (e.g., Chounet et al., [29]).

The grid search method was used to look for the optimal values of the rupture velocity and the rupture’s initial depth. We found that for each nodal plane used as a rupture plane, two or three minima appeared. It is crucial to consider the unique characteristics of each seismic event. Further investigations into the phenomenon of multiple minima and improved methodologies for rupture plane determination will advance our understanding of earthquake source characterization and hazard assessment.

By source rupture modeling and crustal velocity modeling, the 2017 M_W 7.3 Iran earthquake occurred in the high-velocity layers over a low-velocity zone along the boundary between the Arabian plate and the Eurasian plate within the seismic zone of the Zagros. This finding highlights the significance of the seismogenic structure in earthquake occurrence and contributes to our understanding of the earthquake processes in the region, providing valuable insights for future seismic hazard assessments.

Acknowledgments

This research was made possible through the support of the Natural Sciences and Engineering Research Council of Canada under the Discovery Grant program. We would like to express our gratitude to editor Walter Salazar for the improvement of this chapter. We would also like to acknowledge the following programs and tools in our data processing and figure preparation: SAC2000, Rdseed, geotool, MATLAB, and Generic Mapping Tools (GMT). Special thanks go to Lingling Ye at the Department of Earth and Planetary Sciences, University of California, Santa Cruz, California, for providing us with a version of the source rupture modeling code. We are sincerely grateful for her assistance.

Data sources

In this study, the seismograms, GCMT, and USGS nodal plane solutions utilized were obtained from the Incorporated Research Institutions for Seismology (IRIS) database, accessed at http://www.iris.edu (last accessed May 29, 2023). Certain information presented in the introduction section was revised based on a web page of the U.S. Geological Survey (https://earthquake.usgs.gov/earthquakes/eventpage/us2000bmcg#executive; last accessed the 31 July 2018).