Rupture History of the 13th November 2017 M w 7.3 Iran-Iraq Border Region Earthquake Based on Teleseismic Data

We have examined the temporal and spatial slip distribution of the 2017 M w 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 ﬁnite 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 20s, and the total seismic moment was 0.9 × 10 20 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.2m, located on a subfault of the initial rupture. The entire rupture lasted 20s, and 75% of the energy was released in the ﬁrst 10s. The rupture area was 825km 2 , and the estimated static stress drop was 6.1MPa.


Introduction
e 12th November 2017 M w 7.3 earthquake occurred about 200 km northeast of Baghdad, the capital of Iraq, at 18:18:17 UTC [1]. e 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). e 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]. e motion of these two plates produced the Zagros fold and thrust belt. e thrust movement in the northwestern part of the belt caused this earthquake [4]. e 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. ere were 225 intracrustal earthquakes in this area with magnitudes larger than 5.0 (USUS), and there were about 26 earthquakes greater than M w 6.0. For these earthquakes, the depths have a large range of 0∼50 km with magnitudes less than 6.5. e depths are concentrated around 10∼25 km for magnitudes larger than 6.5 ( Figure 1). e region has experienced four earthquakes greater than M 7.0. An M w 7.1 earthquake on May 6, 1930 was located at 44.78°E, 38.09°N. On November 24, 1976, an M w 7.0 earthquake occurred at 43.98°E, 38.88°N. An M w 7.4 earthquake on June 20, 1990 (49.52°E, 36.95°N), caused significant damage [1]. e most recent of these was an M w 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. is 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 M w 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][9][10][11][12][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. e major purpose of this paper was to quickly reveal the rupture history of this earthquake.

Data
Forty-nine broadband P-waves with epicentral distances ranging from 30 to 90°were used for the inversion. e records were from the Global Seismographic Network (GSN) network of IRIS. e 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. e total length of each record was 50 s. e 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.

Inversion Method
e 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. en, 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°).    shock and the centroid moment tensor solution (CMT), ten days aftershocks and historical earthquakes since 1924. e red star is the epicenter of the main shock. e black and pink circles are ten days aftershocks larger than 4.0. e brown, yellow, and red circles are historical earthquakes greater than 5.0. e CMT solution of some historical earthquake larger than 6.5 is shown in this figure. (b) e relationship between magnitude and depth for historical earthquake. e red circle is the 2017 Iran-Iraq border region earthquake. e information of these earthquakes is from USGS. e red line is the plate boundary [5,6].
For a steep dip angle, the two directions are rake angles of 35°a nd 125°. e 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: (1) 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. e observed record is b. e size of A is n × m, x is m × 1, and b is n × 1. e size of n depends on the amount of stations, time step of each record, and record components. e size of m is controlled by the number of subfaults and time windows. e linear least-squares approach can be used to solve equation (1), but A is an ill-conditioned matrix, making the solution unstable. is means that with a small change in matrix A or b, greater variation happens in matrix x. erefore, we appended linear stability constraints for equation (1). e smoothing constraint (S) was utilized to minimize the slip between adjacent subfaults (x i − x i+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]. e solution was constrained to be nonnegative is was because slip cannot regresses. en function (1) was rewritten as where C −1 d 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 ‖Ax − b‖ was the label on the x-axis and the smoothing norm ‖Lx‖ was the y-axis label. L represented the constraints matrix. e 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.
e 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, x i and y i were the observed and synthetic waveform data:

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. e fault sizes were 80 km and 70 km along the strike and dip directions. e hypocenter was located at a depth of 25 km. e 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. e total rupture time was 1.6 s. A complicated rise time function can be constructed on each subfault, if required by the data. e 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. e 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  Table 1.
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    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. e locations of epicenter of these aftershocks have some uncertainties. ese 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. erefore, through residuals of the records, aftershock distribution, and the InSAR image result, we evaluated the slip model with NEdipping as the real fault plane. is 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.

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. e 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. e results are shown in Figure 8. e 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. erefore, in the next analysis, the maximum number of windows was 6. For different windows, the trend of the L-curve was similar. e corner value was about λ � 80. Meanwhile, the slip distribution was also similar, but with different maximum slips.
e 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. e first node of the Fourier spectrum for a triangle source time function was 2/T, where T was the duration. erefore, 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.
With the increase of the source time duration, the residual also gradually became larger. e 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. erefore, 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.

Spatial and Temporal Distribution of Slip
rough the above analysis, we adopted a source rise time of 0.2 s and 6 time windows, each separated by 0.2 s. e 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]. 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. e maximum slip near the initial rupture point was 5.2 m with only one asperity. e 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 km 2 , the estimated static stress drop was 6.1 MPa using the expression Δσ � 7πμΔD/(16R) [12,24,25]. is was consistent with the observed stress drop range of 1∼10 MPa [26]. e depth of the asperity ranged from 21 to 28 km. e 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). is 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]. e slip was both focused on the initial rupture point and dominated by thrust as well as a small amount of right-lateral strike slip. e 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. e total rupture duration was about 20 s. e initial rupture and the upper side subfaults had the maximum slip amplitude. e 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. e 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). ere 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. e locations of these aftershocks were mainly to the southwestern side of the epicenter and in the periphery of the large slip ( Figure 12).

Conclusion
After an earthquake, in order to recover the rupture process in a very short time, the teleseismic data are the best choice. e ray path is simple for teleseismic records with distances of 30∼90°and Green's functions are easy to calculate. e teleseismic record can distinguish shallow from deep slip [9]. In this paper, the spatiotemporal variation of the rupture process of the 2017 M w 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. e rupture history was relatively simple, and the slip concentrated around the        initial rupture point with depths ranging from 21 to 28 km. e total rupture duration was about 20 s, and the seismic moment was 0.9 × 10 20 N·m.

Data Availability
e 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
e authors declare that they have no conflicts of interest.