Abstract
We have examined the temporal and spatial slip distribution of the 2017 Mw 7.3 Iran-Iraq border region earthquake, utilizing 49 broadband teleseismic P-wave records. Based on the nonnegative least square method and multi-time window, a finite fault model was used to parameterize the rupture process. According to the L-curve, the optimal inversion result was detected. The inversion results showed that the earthquake was a shallow-dip thrusting event. Rupture duration was 20 s, and the total seismic moment was 0.9 × 1020 N·m. There was only one asperity in the fault plane, which indicated the rupture process was simple. The slip was mainly distributed around the initial rupture point, dominated by thrust motion with a small amount of right-lateral strike slip, and the maximum slip was 5.2 m, located on a subfault of the initial rupture. The entire rupture lasted 20 s, and 75% of the energy was released in the first 10 s. The rupture area was 825 km2, and the estimated static stress drop was 6.1 MPa.
1. Introduction
The 12th November 2017 Mw 7.3 earthquake occurred about 200 km northeast of Baghdad, the capital of Iraq, at 18:18:17 UTC [1]. The earthquake, located near the area of Iran-Iraq border, caused serious damage over a large region, especially in the town of Sarpol-e Zahab in the Kermanshah province of Iran. In total, about 530 people were killed and more than 7,200 injured (https://baike.baidu.com/item). The seismogenic district was in a mountainous area, causing landslides, which were a serious impediment to the relief work after the earthquake.
In the seismogenic area, the Arabian plate is moving towards the Eurasian plate at a rate of about 20∼25 mm/yr according to the observed global positioning system (GPS) data [2, 3]. The motion of these two plates produced the Zagros fold and thrust belt. The thrust movement in the northwestern part of the belt caused this earthquake [4]. The earthquakes in the area occur mainly along the junction of the Arabian Plate and Eurasian Plate and around the Zagros Mountains (Figure 1). In the USGS catalogue, information about historical earthquakes began in 1924. There were 225 intracrustal earthquakes in this area with magnitudes larger than 5.0 (USUS), and there were about 26 earthquakes greater than Mw 6.0. For these earthquakes, the depths have a large range of 0∼50 km with magnitudes less than 6.5. The depths are concentrated around 10∼25 km for magnitudes larger than 6.5 (Figure 1). The region has experienced four earthquakes greater than M 7.0. An Mw 7.1 earthquake on May 6, 1930 was located at 44.78°E, 38.09°N. On November 24, 1976, an Mw 7.0 earthquake occurred at 43.98°E, 38.88°N. An Mw 7.4 earthquake on June 20, 1990 (49.52°E, 36.95°N), caused significant damage [1]. The most recent of these was an Mw 7.1 earthquake on October 23, 2011 (43.51°E, 38.72°N).
After an earthquake larger than 6.0, the Incorporated Research Institutions for Seismology (IRIS) can provide a teleseismic record in about 2∼3 hours. Some research institutions (e.g., USGS; Institute of Geophysics, China Earthquake Administration) and researchers will give the rupture process of the earthquake in about 3 hours. This is useful for earthquake relief work after the event. Although the existing research results for the rupture process inversion show that slip distribution for one earthquake is largely due to the kind of data used [7], the results based on teleseismic field data are usually relatively stable, neglecting the resolution of this data. For example, in different studies on the 2008 Mw 7.9 Wenchuan earthquake, despite multiple fault planes, fault model variety, data numbers, source functions, and frequency bands, the slip in the fault was mainly located near the Hongkou and Beichuan area [8–13]. In this study, we use broadband teleseismic records to establish the spatial and temporal distribution of slip for the earthquake using a nonnegative least square method and multi-time window. The major purpose of this paper was to quickly reveal the rupture history of this earthquake.
2. Data
Forty-nine broadband P-waves with epicentral distances ranging from 30 to 90° were used for the inversion. The records were from the Global Seismographic Network (GSN) network of IRIS. The data are displacement by deconvolving of the instrument response with a zero-phase bandpass filter between 0.01 and 0.5 Hz with a time step of 0.1 s. The total length of each record was 50 s. The locations of these stations are shown in Figure 2, and the maximum amplitude and epicentral distance are shown in Figure 3. If necessary, the initial time of some records was manually adjusted when the P-wave arrival times derived from velocity model AK135 [14] have timing errors. Most of the P-waveforms were very similar, but with different first motion directions. Meanwhile, the P-waveforms had only a significant large packet with a duration of about 30 s, which indicated that the fault plane had a large slip area.
3. Inversion Method
The rupture history inversion is based on the representation theorem. Here, the nonnegative least-squares inversion method was used to invert the rupture history. Since Hartzell and Heaton [15] proposed this method, large earthquakes throughout several decades from over the world have been successfully explained [16].
Following this method, the fault plane is divided into a series of subfaults with the same size. Based on the ray theory method, the synthetic waveforms of these subfaults are calculated by the computer programs in seismology [17] using velocity structure model Ak135. To reveal the slip direction of each subfault, the direction is divided into two mutual vertical directions. Then, total slip can be composed, and slip can change within the 90° range. Slip amplitude for the two directions of every subfault is solution vector x. Meanwhile, the synthetic waveforms of each subfault for a unit amount of slip in two directions are also calculated (Green’s functions). For the model with a shallow dip angle, the directions are both a right-lateral strike slip dislocation (rake angle 180°) and a thrust dislocation (rake angle 90°). For a steep dip angle, the two directions are rake angles of 35° and 125°. The matrix of Green’s functions is A, and the filter band and sampling interval are the same as recorded. Including observed data matrix b, the over determined system is as follows:
Each column of A represents the synthetic records in all stations for a particular subfault. Slip of this particular subfault is the mutual vertical direction. x is the slip value of every subfault. The observed record is b. The size of A is n × m, x is m × 1, and b is n × 1. The size of n depends on the amount of stations, time step of each record, and record components. The size of m is controlled by the number of subfaults and time windows.
The linear least-squares approach can be used to solve equation (1), but A is an ill-conditioned matrix, making the solution unstable. This means that with a small change in matrix A or b, greater variation happens in matrix x. Therefore, we appended linear stability constraints for equation (1). The smoothing constraint (S) was utilized to minimize the slip between adjacent subfaults (xi − xi+1 = 0) along both strike and downdip, as well as the adjacent windows (multi-time window) for most of the studies on the rupture process inversion. Considering that there was no large slip around the fault, the constraint (B) was added to constrain the peripheral area of the fault beside the surface area [12]. The solution was constrained to be nonnegative (xi ≥ 0). This was because slip cannot regresses. Then function (1) was rewritten aswhere is a priori data covariance matrix to set each record in the inversion accounts to the same weight. To account for the heterogeneity of the rupture duration, slip amplitude, and direction, a multi-time window was adopted. Any time window can slip, and a complicated rise time form can be composed. A constant rupture front velocity was used, and if slip occurred in subsequent time windows, it was interpreted as a lower rupture velocity, otherwise equal to the set value.
Linear weights λ1 and λ2 control the trade-off between satisfying these constraints and fitting the data. For inversion problems, the important thing is to determine the best smoothing weight. Here, following Hansen [18] and Mendoza and Hartzell [19], we utilized L-curve analysis to distinguish optimal value. For the L-curve, the residual norm was the label on the x-axis and the smoothing norm was the y-axis label. L represented the constraints matrix. The L-curve represented the trade-off between the residual norm and constraint norm for different smoothing weights. Mendoza [19] used a smoothing matrix and minimizing total seismic moment matrix. Here, we used a smoothing matrix.
The value near corner of the L-curve, which can both meet the smoothing constraint and fit the recorded data, was an optimal one. Like the judgment criterion used by Nakamura et al. [10], the best estimated result was obtained by the minimum value of the waveform misfit (equation (4)). Here, xi and yi were the observed and synthetic waveform data:
4. Rupture Fault Model
In order to distinguish the seismogenic fault, two models with different strike and dip angles were utilized according to the USGS Centroid Moment Tensor (CMT) solution. For models I and II, the strike and dip angles were 351°/16° and 122°/79°, respectively. The fault sizes were 80 km and 70 km along the strike and dip directions. The hypocenter was located at a depth of 25 km. The fault plane was composed of 224 subfaults, each with a fixed dimension of 5 km by 5 km. In order to simulate the possible rise time, eight isosceles triangle time windows were used, 0.2 s in duration and separated from one another by 0.2 s. The total rupture time was 1.6 s. A complicated rise time function can be constructed on each subfault, if required by the data. The rupture velocity was 2.6 km/s, which was about 0.8 times the shear-wave velocities [20]. For model I with a shallow dip angle, the shallowest fault was at a depth of 15 km to reduce the number of subfaults and contain the area of the initial rupture. The depth range of the fault was 15.6∼33.6 km. For model II with a steep dip angle, the fault started from the surface area and had a maximum depth of 66 km. For various smoothing weights, the L-curves of models I and II are shown in Figure 4. Cross stars represent the corner points (λ = 80). According to Mendoza [19], the values were 2∼3 times higher or lower than the corner value, yielding a similar residual as the 2010 Mw 7.0 Baja California earthquake. Here, the residual and moment of λ from 40 to 120 are listed in Table 1.
(a)
(b)
First, the residual increased about 5% for models I and II when λ ranges from 40 to 120. Slip distributions of the smoothing weight (λ = 40, 80, and 120) for models I and II are shown in Figure 5 and Table 1. The results showed similar slip distribution patterns for different smoothing constraints. Meanwhile, the slip area was more continuous and uniform with the increase of the constraint factor. For λ = 40, the dislocation distribution was discrete, and the maximum slip was 7.6 and 5.7 m for models I and II. For λ = 80 and 120, the slip areas had small differences. Secondly, for the same constraint weight, the residual of model II was larger than model I (Figure 6). For λ = 80, the residual was about 5% larger, but the resolution of the teleseismic data tended to be lower. Such a degree of difference for the residual could not recognize the real rupture model without additional constraints. Projection of the spatial-slip distribution and ten day aftershock distributions for models I and II are shown in Figure 7. The locations of epicenter of these aftershocks have some uncertainties. These aftershocks mainly occurred along the strike for model I. On the contrary, the aftershocks were nearly perpendicular to the strike for model II. Additionally, Kobayashi [6] obtained the real rupture model with a northeast (NE) dipping plane according to interferometric synthetic aperture radar (InSAR) images for this earthquake. Therefore, through residuals of the records, aftershock distribution, and the InSAR image result, we evaluated the slip model with NE-dipping as the real fault plane. This agreed with the results of the USGS and Wang et al. [21], but this was different from the results of Institute of Geophysics, China Earthquake Administration. The strike of the fault model was 123° (http://www.cea-igp.ac.cn/tpxw/276001.html).
(a)
(b)
5. Rise Time of Source Function
Time window quantities depend on the maximum dislocation duration. For kinematic source inversion, the number of windows is an important influential factor. The more windows there are for inversion, the smaller the values of residual will be, but the result will be unstable with more unknown quantities. In order to detect the rise time of this earthquake, 1, 2, 4, 6, 8, 10, and 12 windows with duration of 0.2 s each lagged from the previous one by 0.2 s are utilized. With these situations, the total rise times were 0.2, 0.4, 0.8, 1.2, 1.6, 2.0, and 2.4 s, respectively. The results are shown in Figure 8. Figure 8(a) is change of residual, Figure 8(b) is L-curve, Figure 8(c) is slip distribution, and Figure 8(d) represented comparison of observed and synthetic waveforms.
(a)
(b)
(c)
(d)
The results showed that the multiple time windows reduced the residual. For windows from 1 to 12, the residual was reduced by 8%. When the window was more than 6, the residual had little change. Therefore, in the next analysis, the maximum number of windows was 6. For different windows, the trend of the L-curve was similar. The corner value was about λ = 80. Meanwhile, the slip distribution was also similar, but with different maximum slips. The synthetic records were in good agreement with the observed ones for various windows. Here, we needed to focus on the situation where the synthetic records had very small differences for different numbers of windows. It meant that major slip of the subfaults in asperity happened in the first time window, so for only one window, the synthetic records showed good agreement with the observed ones.
For the rupture process inversion, when deciding on the duration of the source function, it was important to consider the first node of its Fourier spectrum. For the source time function, the first node of the Fourier spectrum should be larger than the maximum frequency [22, 23]. Here, we took source functions with different durations of 0.2, 0.4, 0.8, 1.2, 1.6, 2.0, 2.4, 2.8, and 3.2 s. Six time windows were used, and each was separated by 0.2 s. The first node of the Fourier spectrum for a triangle source time function was 2/T, where T was the duration. Therefore, the first nodes were 10.0, 5, 2.5, 1.7, 1.3, 1.0, 0.8, 0.7, and 0.6 Hz, respectively (Figure 9). For the inversion, the maximum frequency was 0.5 Hz.
(a)
(b)
(c)
With the increase of the source time duration, the residual also gradually became larger. The duration was less than 2.0 s with the first node of the Fourier spectrum larger than 1.0 Hz, and the difference of the residual was within 8%. However, the residual increased significantly when the duration was larger than 2.0 s. When the duration was 3.2 s, the residual was 24% larger than when the duration is 0.2 s. From the frequency domain, for the inversion with a maximum frequency of 0.5 Hz, when the first node of the Fourier spectrum was less than 1 Hz (larger than 2.0 s duration), the degree of conformity for the synthetic and observed records was reduced. Therefore, for the source time function, the first node of the Fourier spectrum should be larger than the maximum frequency. Using near field records, Hartzell et al. [23] also concluded that if the spectrum of source time function was relatively large below the maximum frequency for the inversion; the result had little change.
6. Spatial and Temporal Distribution of Slip
Through the above analysis, we adopted a source rise time of 0.2 s and 6 time windows, each separated by 0.2 s. The results are exhibited in Figures 3 and 10–12. An asperity can be defined as enclosing fault elements that have a slip 1.5 times larger than the average slip over the fault [20].
(a)
(b)
(c)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
Figure 10(a) is the slip distribution, Figure 10(b) is the moment rate function, and Figure 10(c) represents the slip amplitude of every window of the subfaults in the blue rectangle in Figure 10(a). In the fault plane, the spatial slip distribution was concentrated around the hypocenter. The maximum slip near the initial rupture point was 5.2 m with only one asperity. The asperity was mainly a thrust with a small amount of right-lateral strike slip. Suppose that the fault was a circle and the asperity had an average slip (ΔD) of 1.5 m and a radius (R) of 16.2 km with total area 825 km2, the estimated static stress drop was 6.1 MPa using the expression Δσ = 7πμΔD/(16R) [12, 24, 25]. This was consistent with the observed stress drop range of 1∼10 MPa [26]. The depth of the asperity ranged from 21 to 28 km. The total seismic moment was 0.90 × 1020 N·m, which was smaller than the results found by the USGS (1.1 × 1020 N·m) and the Harvard CMT solution (1.7 × 1020 N·m). This was probably due to the maximum frequency of the inversion and the different records used. Meanwhile, slip did not rupture to the surface. On the whole, the slip distribution was consistent with the results of the USGS and Wang et al. [21]. The slip was both focused on the initial rupture point and dominated by thrust as well as a small amount of right-lateral strike slip. The moment rate function showed that 75% of the energy was released in the first 10 s, corresponding to the slip near the initial rupture point (Figure 11). After 10 s, the remaining 25% of the energy was generated, corresponding to the slip near the fault edge (Figure 11). Meanwhile, the slip value at the periphery of the fault was small. The total rupture duration was about 20 s. The initial rupture and the upper side subfaults had the maximum slip amplitude. The slip happened in the first three time windows with a duration 0.6 s, and the first window had the maximum value. Slip of the windows along the other subfaults had a smaller value. The slip durations in some subfaults were 1.2 s.
Overall, the synthetic record fitted the observed ones fairly well, both in amplitude and phase (Figure 3). There were 25 aftershocks with magnitudes of 4.0∼5.0 around the epicenter within 10 days of the main shock and 2 aftershocks with magnitudes of 5.0∼6.0. The locations of these aftershocks were mainly to the southwestern side of the epicenter and in the periphery of the large slip (Figure 12).
7. Conclusion
After an earthquake, in order to recover the rupture process in a very short time, the teleseismic data are the best choice. The ray path is simple for teleseismic records with distances of 30∼90° and Green’s functions are easy to calculate. The teleseismic record can distinguish shallow from deep slip [9]. In this paper, the spatiotemporal variation of the rupture process of the 2017 Mw 7.3 Iran-Iraq border region earthquake was revealed by 49 broadband P-wave teleseismic records. Common features of our inversion result showed that this earthquake was a shallow-dip thrusting event with a small amount of right-lateral strike slip. The rupture history was relatively simple, and the slip concentrated around the initial rupture point with depths ranging from 21 to 28 km. The total rupture duration was about 20 s, and the seismic moment was 0.9 × 1020 N·m.
Data Availability
The teleseismic records used for the inversion are from the GSN of Incorporated Research Institutions for Seismology (IRIS) at http://ds.iris.edu/wilber3/find_event.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Authors’ Contributions
Yin Deyu collected and processed the inversion data, performed inversion, and analyzed the inversion result of the 2017 Mw 7.3 Iran-Iraq border region earthquake. Yin Deyu wrote the manuscript. Qifang Liu participated in writing the computer program for the inversion method and analyzed the inversion result. Yuexin She helped to draft the manuscript and participated in the design of this inversion. Yun Dong helped to draft the manuscript.
Acknowledgments
The National Natural Science Foundation of China (51804129) supported the study. The National Natural Science Foundation of China (51708245) and the Natural Science Foundation for Colleges and Universities in Jiangsu Province (17KJB620002) also supported and provided the data for this inversion. The project with Science and Technology of Huai’an City (HAG201606), the Natural Science Foundation of Jiangsu Province (BK20160426), and Huaian Natural Science Project (HAB201708) helped to write the manuscript. Science and Technology Guidance Project of Jiangsu Housing and Urban-Rural Construction Department (2018ZD034 and 2018ZD154) also provided fund to the study. The authors thank Wu Jingke for his contribution to this article.