Rheology of the Zagros Lithosphere from Post-Seismic Deformation of the 2017 Mw7.3 Kermanshah, Iraq, Earthquake

We use 2018–2020 Sentinel-1 InSAR time series data to study post-seismic deformation processes following the 2017 Mw 7.3 Kermanshah, Iraq earthquake. We remove displacements caused by two large aftershock sequences from the displacement field. We find that for a six month period the response is dominated by afterslip along the up-dip extension of the coseismic rupture zone, producing up to 6 cm of radar line-of-sight displacements. The moment magnitude of afterslip is Mw 5.9 or 12% of the mainshock moment. After that period, the displacement field is best explained by viscoelastic relaxation and a lower crustal viscosity of ηlc = 1 +0.8 −0.4 × 10 19 Pas. The viscosity of the uppermost mantle is not constrained by the data, except that it is larger than 0.6 × 1019 Pas. The relatively high lower crustal and uppermost mantle viscosities are consistent with a cold and dry lithosphere of the Zagros region.


Introduction
The Mw 7.3 Kermanshah earthquake of 12 November 2017 near the Iran/Iraq boundary occurred along a thrust fault of the Zagros Mountains. It killed more than 400 people and destroyed 12,000 buildings. There were serious economic losses caused by this event [1]. It was one of the largest earthquakes in the northwestern Zagros during the instrumental era (others were the 1909 Mw 7.3 earthquake and 1957 Mw 7.1 earthquake) [2,3]. Observations of post-seismic deformation provide an opportunity to investigate the fault and rheological properties of the Zagros lithosphere.
Post-seismic displacements can be caused by (i) afterslip along the fault, either along the patch that ruptured during the earthquake or next to it [4][5][6][7]; (ii) viscoelastic relaxation of the lower crust and/or upper mantle [8][9][10]; (iii) poroelastic rebound following earthquake-induced pore-fluid pressure changes [11][12][13]; or (iv) large aftershocks. These processes have different spatial-temporal characteristics [14]. Afterslip and poroelastic rebound cause ground displacements in the near field and last for a few months to years. In contrast, viscoelastic relaxation causes displacements in the far field and can last several decades. Understanding the origin of post-seismic deformation is important for seismic hazard because it changes the crustal stress field, loading or unloading nearby faults [15][16][17].
In this study we use 2.3 years (until February 2020) of Sentinel-1 data to investigate post-seismic displacements after the Kermanshah earthquake. After removing the coseismic deformation of the largest aftershocks, we investigate whether surface displacements reflect afterslip and/or viscoelastic relaxation and explore the afterslip characteristics and the lithospheric rheology of this region.

Geologic Setting and Aftershock Sequence
The Kermanshah earthquake is located on the southwestern edge of the northwestern Zagros fold-and-thrust belt (ZFTB), one of the most seismically active intra-continental belts and relatively young orogens in the world [8,18]. ZFTB is formed due to continental collision between the Arabian and Eurasian plates in western Asia. The convergence rate between these two plates is approximatelỹ 2 cm/year [19] and ZFTB accommodates one third of the total collision rate, according to the NUVEL-1A plate motion model [20]. GPS data show that the collision rate decreases steadily from 9 mm/year in the southeastern section to 7 mm/year and 4 mm/year in the central and northwestern Zagros, respectively [20]. The Zagros Main Recent Fault (MRF), the High Zagros Fault (HZF) and the Zagros Mountain Front Fault (MFF) are major faults in the northwestern Zagros region ( Figure 1).
Several studies using seismic and geodetic data concluded that this earthquake occurred on a shallowly east-dipping blind fault with 15 • dip angle and between 13 and 20 km depth [3,18,[20][21][22][23][24]. For the viscoelastic relaxation simulation below we use the source parameters of Vajedian [3].
The Kermanshah earthquake was followed by strong aftershock sequences. There were two sequences of aftershocks with several events with magnitude 5.0 or larger that occurred in close vicinity within a short time period. The first one was a one-hour-long sequence on 11 January 2018 with magnitude ranging from Mw 5.1 to Mw 5.5; the second one was a seven-hour-long sequence on 25 November 2018 with magnitude ranging from Mw 5.0 to Mw 6.3. These sequences ( Figure 1, Table 1) are referred to as aftershock sequence 1 and 2, respectively. The third largest aftershock was an Mw 6.0 on 25 August 2018. In this study we use 2.3 years (until February 2020) of Sentinel-1 data to investigate post-seismic displacements after the Kermanshah earthquake. After removing the coseismic deformation of the largest aftershocks, we investigate whether surface displacements reflect afterslip and/or viscoelastic relaxation and explore the afterslip characteristics and the lithospheric rheology of this region.

Geologic Setting and Aftershock Sequence
The Kermanshah earthquake is located on the southwestern edge of the northwestern Zagros fold-and-thrust belt (ZFTB), one of the most seismically active intra-continental belts and relatively young orogens in the world [8,18]. ZFTB is formed due to continental collision between the Arabian and Eurasian plates in western Asia. The convergence rate between these two plates is approximately ~2 cm/year [19] and ZFTB accommodates one third of the total collision rate, according to the NUVEL-1A plate motion model [20]. GPS data show that the collision rate decreases steadily from 9 mm/year in the southeastern section to 7 mm/year and 4 mm/year in the central and northwestern Zagros, respectively [20]. The Zagros Main Recent Fault (MRF), the High Zagros Fault (HZF) and the Zagros Mountain Front Fault (MFF) are major faults in the northwestern Zagros region ( Figure 1).
Several studies using seismic and geodetic data concluded that this earthquake occurred on a shallowly east-dipping blind fault with 15° dip angle and between 13 and 20 km depth [3,18,[20][21][22][23][24]. For the viscoelastic relaxation simulation below we use the source parameters of Vajedian [3].
The Kermanshah earthquake was followed by strong aftershock sequences. There were two sequences of aftershocks with several events with magnitude 5.0 or larger that occurred in close vicinity within a short time period. The first one was a one-hour-long sequence on 11 January 2018 with magnitude ranging from Mw 5.1 to Mw 5.5; the second one was a seven-hour-long sequence on 25 November 2018 with magnitude ranging from Mw 5.0 to Mw 6.3. These sequences ( Figure 1, Table  1) are referred to as aftershock sequence 1 and 2, respectively. The third largest aftershock was an Mw 6.0 on 25 Aug 2018.

InSAR Data and Processing Methodology
We used a total of 134 ascending Sentinel-1 A/B data from Path 72 (from 17 November 2017 to 5 February 2020; the first image 5 days after the mainshock) and 54 descending images from Path 79 (from 18 November 2017 to 31 January 2020; the first image 6 days after the mainshock).
We used ISCE stack processor [25,26] to generate the interferograms with 21 looks in the azimuth and 7 looks in the range direction. Precise orbit data and the 3 arc-second DEM from the Shuttle Radar Topography Mission (SRTM) [27] were used to simulate and remove the phase error caused by topography and earth curvature from each interferogram. Multi-looked and filtered interferograms were registered to a master SAR image by finding offsets between single-look complex (SLC) images and the master SLC using DEM and orbit vectors. Finally, the phases of the coregistered interferograms were unwrapped using the statistical-cost network-flow algorithm (SNAPHU) [28].
For time series processing we used the routine workflow of the Miami InSAR time series software in Python (MintPy) [29]. In this workflow the network of interferograms was inverted for the raw phase time series and then corrected for the tropospheric delay (we used the ERA5 global atmospheric model) and topographic residuals to obtain the noise-reduced displacement time series from which the average line of sight (LOS) velocity was estimated pixel by pixel.

Modelling Approach
To estimate the aftershock and afterslip parameters, we consider uniform rectangular dislocations in a homogeneous elastic half-space [30]. The model parameters are fault location, length, width, dip angle, strike angle, depth, strike slip and dip slip. To estimate the rheological structure we consider Maxwell rheology for the lower crust and upper mantle, each characterized by two parameters, the steady state shear modulus and steady state viscosity. We use the RELAX software [31][32][33], which uses the Fourier-domain elastic Green's functions and the equivalent body-force representation of coseismic and post-seismic deformation processes to calculate the viscoelastic relaxation displacement. The model parameters are the lower crust and upper mantle viscosities.
We sample from the data using a downsampling method combining uniform and quadtree sampling. We use quadtree sampling in the region with significant deformation and in the far field the uniform sampling approach (see Supplemental Material for downsampling results). We calculate for each model a misfit function χ 2 defined as: where d obs and d sim are the observed and simulated displacements, respectively, and C is the covariance matrix [34]. We solve the non-linear inversion problem for the elastic dislocations using a Bayesian approach with the Geodetic Bayesian Inversion Software (GBIS) [34]. We also use GBIS to estimate C. We invert for the rheological structure using a grid search approach. Figure 2 shows the line-of-sight (LOS) displacements from ascending and descending orbits after the earthquake (17 and 18 November 2017, respectively) to February 2020, where the positive value means movement towards the sensors and negative values mean movement away from the sensor. The ascending data shows a lobe of up to 8 cm of LOS decrease (yellow/red color) and the descending data shows lobes of both LOS decrease and LOS increase (up to 10 cm, blue color). The ascending and descending displacement time series for points in the epicentral area show initially rapidly and then slowly growing and decaying signals, respectively ( Figure 2c). We therefore consider two time periods characterized by different displacement patterns, from 17 November 2017 to 10 June 2018 (time period 1) and from 10 June 2018 to February 2020 (time period 2). We consider four different areas (Figure 2a,b). In addition to the area used for studying viscoelastic relaxation (area 1), we consider two aftershock areas (areas 2 and 3) and an area to investigate afterslip (area 4).

Post-Seismic Displacement Field
Remote Sens. 2019, 11, x FOR PEER REVIEW 4 of 14 Figure 2 shows the line-of-sight (LOS) displacements from ascending and descending orbits after the earthquake (17 and 18 November 2017, respectively) to February 2020, where the positive value means movement towards the sensors and negative values mean movement away from the sensor. The ascending data shows a lobe of up to 8 cm of LOS decrease (yellow/red color) and the descending data shows lobes of both LOS decrease and LOS increase (up to 10 cm, blue color). The ascending and descending displacement time series for points in the epicentral area show initially rapidly and then slowly growing and decaying signals, respectively ( Figure 2c). We therefore consider two time periods characterized by different displacement patterns, from 17 November 2017 to 10 June 2018 (time period 1) and from 10 June 2018 to February 2020 (time period 2). We consider four different areas (Figure 2a,b). In addition to the area used for studying viscoelastic relaxation (area 1), we consider two aftershock areas (areas 2 and 3) and an area to investigate afterslip (area 4).

Aftershocks
We consider both aftershock sequences as composite events and we obtain elastic dislocation models for aftershock sequence 1 and 2 so that their contributions can be removed from the post-seismic displacement fields. We do not consider the 25 August 2018 aftershock located east of the epicenter because there is no clear signal associated with it in ascending and descending data ( Figure S1) and no geodetic focal mechanism can be estimated. For both ascending and descending orbits, we use interfergrams reconstructed from the raw time series. They are less noisy than the

Aftershocks
We consider both aftershock sequences as composite events and we obtain elastic dislocation models for aftershock sequence 1 and 2 so that their contributions can be removed from the post-seismic displacement fields. We do not consider the 25 August 2018 aftershock located east of the epicenter because there is no clear signal associated with it in ascending and descending data ( Figure S1) and no geodetic focal mechanism can be estimated. For both ascending and descending orbits, we use interfergrams reconstructed from the raw time series. They are less noisy than the observed interferograms because the network inversion of interferograms filters out temporal decorrelation noise.
For aftershock sequence 1, the best-fitting dislocation has a length and width of 12.9 km and 5.4 km and 0.4 m slip. For aftershock sequence 2, the best-fitting dislocation has a length and width of 12.7 km and 8.8 km and 1.8 m slip ( Table 2). The observations, modelled displacement and difference between observations and model are shown in Figure 3a,b. The red lines are nearby faults. The geodetic focal mechanisms (shown in Figure 3) are very similar to the Centroid-Moment-Tensor (CMT) solutions (shown in Figure 1).

Afterslip for Time Period 1
To test whether the surface displacements are consistent with afterslip, we use the ascending and descending data until June 2018 (Figure 4a) after removal of the contribution from the January 2018 aftershock sequence. We estimate the best-fitting uniform dislocation using a homogeneous elastic half space model. The best-fitting solution ( Table 2,  observed interferograms because the network inversion of interferograms filters out temporal decorrelation noise. For aftershock sequence 1, the best-fitting dislocation has a length and width of 12.9 km and 5.4 km and 0.4 m slip. For aftershock sequence 2, the best-fitting dislocation has a length and width of 12.7 km and 8.8 km and 1.8 m slip ( Table 2). The observations, modelled displacement and difference between observations and model are shown in Figure 3a,b. The red lines are nearby faults. The geodetic focal mechanisms (shown in Figure 3) are very similar to the Centroid-Moment-Tensor (CMT) solutions (shown in Figure 1).

Afterslip for Time Period 1
To test whether the surface displacements are consistent with afterslip, we use the ascending and descending data until June 2018 (Figure 4a) after removal of the contribution from the January 2018 aftershock sequence. We estimate the best-fitting uniform dislocation using a homogeneous elastic half space model. The best-fitting solution ( Table 2, black solid rectangle in Figure 4a

Viscoelastic Relaxation for Time Period 2
We assume that the displacements in time period 2 are caused by viscoelastic relaxation. In order to constrain the rheological structure we use a 25 km-thick elastic upper crust and a 20 km-thick viscoelastic Maxwell lower crust over a viscoelastic Maxwell upper mantle (Figure 5a; the average depth of the regional Moho is 45 km according to the CRUST 1.0 model) [35,36]. In the grid search the viscosities of the lower crust and upper mantle range from 1 × 10 17 Pas to 1 × 10 21 Pas; we use a step size for the logarithm of viscosity of 0.25.
Our result (Figure 5b) yields the best estimates of the parameters η lc = 1 +0.8 −0.4 × 10 19 Pas, where the confidence interval is calculated using F-Test [37][38][39]. Our data resolve the viscosity of the lower crust but not the viscosity of the uppermost mantle, except that it is larger than 0.6 × 10 19 Pas.
The comparison between the observed and simulated displacement (Figure 4b) shows that the model explains the descending data well but there is a discrepancy for the ascending data. This discrepancy can be attributed to the larger tropospheric noise for the ascending data which are acquired during the day when the tropospheric variability is high (local acquisition time of 15:00 and 3:00 for ascending and descending data, respectively). What is more, the color discontinuity in the ascending data (around center-bottom of Figure 4b, top) is caused by the subtraction of coseismic displacement of aftershock sequence 2.
Furthermore, the viscoelastic relaxation signal during this 1.75 year period is relatively small (up to 4 cm). The small signal is a consequence of the fault geometry. A horizontal fault causes largely horizontal flow.
We also explored afterslip models and concluded that they can be ruled out. The fault plane of the best-fitting model locates a few kilometers above the coseismic rupture in the phanerozoic cover ( Figures S2 and S3), which is geologically not plausible. Besides, the best-fitting viscoelastic relaxation model is characterized by a significantly better rms (rms = n i=1 residual 2 i n ) than the best-fitting afterslip model (0.015 m versus 0.025 m). We can also rule out afterslip on the coseismic rupture or its downdip extension because of even higher rms. We did not consider combined models with both afterslip and viscoelastic relaxation because the data lack the resolution to constrain all model parameters.

Viscoelastic Relaxation for Time Period 2
We assume that the displacements in time period 2 are caused by viscoelastic relaxation. In order to constrain the rheological structure we use a 25 km-thick elastic upper crust and a 20 km-thick viscoelastic Maxwell lower crust over a viscoelastic Maxwell upper mantle (Figure 5a; the average depth of the regional Moho is 45 km according to the CRUST 1.0 model) [35,36]. In the grid search the viscosities of the lower crust and upper mantle range from 1 × 10 to 1 × 10 ; we use a step size for the logarithm of viscosity of 0.25.
Our result (Figure 5b) yields the best estimates of the parameters = 1 . . × 10 , where the confidence interval is calculated using F-Test [37][38][39]. Our data resolve the viscosity of the lower crust but not the viscosity of the uppermost mantle, except that it is larger than 0.6 × 10 . The comparison between the observed and simulated displacement (Figure 4b) shows that the model explains the descending data well but there is a discrepancy for the ascending data. This discrepancy can be attributed to the larger tropospheric noise for the ascending data which are acquired during the day when the tropospheric variability is high (local acquisition time of 15:00 and 3:00 for ascending and descending data, respectively). What is more, the color discontinuity in the ascending data (around center-bottom of Figure 4b, top) is caused by the subtraction of coseismic displacement of aftershock sequence 2.
Furthermore, the viscoelastic relaxation signal during this 1.75 year period is relatively small (up to 4 cm). The small signal is a consequence of the fault geometry. A horizontal fault causes largely horizontal flow.
We also explored afterslip models and concluded that they can be ruled out. The fault plane of the best-fitting model locates a few kilometers above the coseismic rupture in the phanerozoic cover ( Figures S2, S3), which is geologically not plausible. Besides, the best-fitting viscoelastic relaxation model is characterized by a significantly better rms ( = ∑ ) than the best-fitting afterslip model (0.015 m versus 0.025 m). We can also rule out afterslip on the coseismic rupture or its downdip extension because of even higher rms. We did not consider combined models with both afterslip and viscoelastic relaxation because the data lack the resolution to constrain all model parameters. The yellow star is model 1 (best-fitting model) in Figure 6; the yellow square is model 2 in Figure 6; the yellow triangle is model 3 in Figure 6. The black dash line in (a) is the coseismic fault. The black solid line in (a) is the afterslip fault. The gray-shaded area in (b) denotes 95% confidence interval.

Summary of Post-Seismic Models
To visualize the fits of the post-seismic models we combine the ascending and descending displacement fields for the two time periods to the quasi-vertical and quasi-horizontal displacements The yellow star is model 1 (best-fitting model) in Figure 6; the yellow square is model 2 in Figure 6; the yellow triangle is model 3 in Figure 6. The black dash line in (a) is the coseismic fault. The black solid line in (a) is the afterslip fault. The gray-shaded area in (b) denotes 95% confidence interval.

Summary of Post-Seismic Models
To visualize the fits of the post-seismic models we combine the ascending and descending displacement fields for the two time periods to the quasi-vertical and quasi-horizontal displacements in east-west direction [40] (Figure 6a,b). The north-south displacement cannot be calculated because only right-looking InSAR data are available and because polar-orbiting sensors are not sensitive to north-south displacements. We generate simulations for models with lower and higher viscosities in both the lower crust and uppermost mantle. We use η lc = η um =10 18 Pas for the low viscosity model (model 2) and η lc = η um =10 20 Pas for the high viscosity model (model 3), respectively.
For time period 1, the low viscosity model produces vertical displacements comparable to the observations, but the east-west displacements are different than observed (model 2, Figure 6a). This reinforces that afterslip was the dominating process during this period. For time period 2, the low and high viscosity models predict significantly more and less displacements than observed, respectively (models 2, 3 in Figure 6f,g), reinforcing that the viscosities are of the order of 10 19 Pas (model 1, Figure 6d). in east-west direction [40] (Figure 6a,b). The north-south displacement cannot be calculated because only right-looking InSAR data are available and because polar-orbiting sensors are not sensitive to north-south displacements. We generate simulations for models with lower and higher viscosities in both the lower crust and uppermost mantle. We use = =10 for the low viscosity model (model 2) and = =10 for the high viscosity model (model 3), respectively. For time period 1, the low viscosity model produces vertical displacements comparable to the observations, but the east-west displacements are different than observed (model 2, Figure 6a). This reinforces that afterslip was the dominating process during this period. For time period 2, the low and high viscosity models predict significantly more and less displacements than observed, respectively (models 2, 3 in Figure 6f,g), reinforcing that the viscosities are of the order of 10 (model 1, Figure 6d).

Afterslip for Time Period 1
Our analysis shows that the post-seismic response through about June 2018 was dominated by afterslip along the up-dip continuation of the coseismic rupture, consistent with Barnhart [18] and Yang [24] who studied shorter time periods. The moment released by afterslip was 1.13 × 10 19 Nm, which was about 12% of that by the coseismic slip. To put these results into perspective, we compiled the location of afterslip and the ratios between afterslip and coseismic moment release found in previous studies of earthquakes, which are listed in the supplementary Table S1.
The afterslip on the up-dip continuation of the coseismic rupture indicates a zone along the fault with velocity-strengthening frictional behavior [41] above a zone with velocity weakening behavior that ruptured during the earthquake. Up-dip afterslip was also observed for the 2003 Mw 6.9 Boumerdes [42] and the 2003 Mw6.8 Zemmouri, Algeria [43] earthquakes. For many other earthquakes, afterslip occurs on the down-dip extension of the coseismic rupture zone. Examples include the 2015 Mw 7.8 Gorkha earthquake [44], the 2011 Mw 9.0 Tohoku earthquake [45][46][47] and the 2008 Mw 7.9 Wenchuan earthquake [48].
The afterslip duration time of six months is comparable to that of the 2014 Napa earthquake [10] but much smaller than for example for the 2004 Mw6.0 Parkfield, California, earthquake [49] after which it continued for about 12 years. The relative amount of afterslip moment release is at the lower end compared to other earthquakes (10%, 28% and 56% for the Gorkha, Tohoku and Pakistan earthquakes, respectively).
Afterslip is promoted by clay-rich sediments/fault with clay-rich gouge, high temperature and elevated pore-fluid pressure [7]. Therefore, the low heat flow (average of 74 mW/m 2 [50]) and relatively low pore-fluid pressure (indicated by the relative low Vp/Vs ratio about 1.73 [51]) might contribute to the relative short duration time and small moment of Kermanshah afterslip compared to the mainshock moment.
Laboratory experiments and previous studies show that lithospheric viscosity is mainly dependent on the composition, temperature and water content [69][70][71]. Lower temperature and less water lead to higher viscosity. The relatively high upper mantle viscosity of the Kermanshah region compared to the western US (Landers, Hector Mine, Central Nevada Seismic Belt, EL Mayor Cucapah, Heban Lake earthquakes), South central Alsaka (Denali earthquake) and Iceland is consistent with the relatively lower average heat flow (Iran: 74 mW/m 2 [50], Western US: 91 mW/m 2 [50], South central Alsaka [72]: 89 mW/m 2 ; Iceland: 175 mW/m 2 [73]). It is also consistent with the presence of the high-velocity anomalies in the upper mantle beneath the Zagros [74] which mark the presence of dry and cold lithospheric mantle [48]. What is more, the relatively high viscosity of the lower crust is also consistent with the relatively low average Vp/Vs value (Iran: 1.73 [51]; Chi-Chi: 1.9 [75]), which can be an indication of low fluid content [76]. Lower temperatures might also contribute to a higher viscosity in the lower crust.
There was no significant relaxation deformation in the first six months, which can be explained by the relatively high lower crust and upper mantle viscosities. Furthermore, laboratory experiments also suggest that the power-law viscous flow often occurs on the hot lithosphere rocks [77], thus the relatively cold lithosphere rocks in our study region may mean the low possibility of the power-law viscous flow occurring in the first few months, which may also contribute to this no significant relaxation deformation in the first six months phenomenon.

Conclusions
We obtained the post-seismic displacement field of the Kermanshah earthquake for the November 2017 to February 2020 time period, longer than previous studies (4 months for Barnhart's study, 7 months for Yang's study), from Sentinel-1 ascending and descending data. We came to the following conclusions: For the first six months (until June 2018) the post-seismic displacement was due to the afterslip along the up-dip extension of the rupture zone. The moment released by afterslip was 12% of that released by coseismic slip. This was at the lower end of the observed range of afterslip durations and relative moment release. Possible explanations for little afterslip were the relatively low heat flow and pore-fluid pressures in the Zagros region.
For the subsequent period (we considered data until February 2020), the post-seismic displacement field was not consistent with afterslip but more consistent with viscoelastic relaxation in the lower crust. Assuming that this was the only process causing surface displacements, we found a best-fitting viscosity of η lc = 1 +0.8 −0.4 × 10 19 Pas for lower crust. The data do not have the resolution to constrain the viscosity of the uppermost mantle but we could infer a lower bound of 0.6 × 10 19 Pas. A relatively strong lower crust and upper mantle was consistent with the relatively lower heat flow, low average Vp/Vs value and the existence of the upper mantle high-velocity anomalies in the Zagros region.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/12/2032/s1: Figure S1: coseismic displacement field caused by the aftershock of 25 August 2018; Figure S2: uniform slip inversion result for afterslip for time period 2; Figure S3: conceptual model of the spatial relationship between coseismic fault slip and afterslip for time period 2; Figure S4: sampled data used for modelling of aftershock sequence 1 and aftershock sequence 2; Figure S5: same as Figure S4 but for time period 1 and time period 2; Table  S1: reported ratio of the afterslip moment release relative to the coseismic moment release and the position of afterslip relative to the coseismic slip; Table S2: rheologic structures inferred in previous studies; Table S3: range, optimal and uncertainties of InSAR-inversion fault slip parameters for the afterslip for time period 2.