**1. 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–7]; (ii) viscoelastic relaxation of the lower crust and/or upper mantle [8–10]; (iii) poroelastic rebound following earthquake-induced pore-fluid pressure changes [11–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–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.

#### **2. 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–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.

**Figure 1.** Map of the study area. The red star represents the location of the mainshock (from USGS earthquake Catalog). The black dashed rectangle marks the coseismic fault [3]. The black dot indicates the upper edge of fault. The red lines are nearby faults (MRF: the Zagros Main Recent Fault, HZF: the High Zagros Fault, MFF: the Zagros Mountain Front Fault). Yellow circles mark aftershocks (Mw > 3.0) up to 839 days after the mainshock (from USGS earthquake catalog: https: //earthquake.usgs.gov/earthquakes/).


**Table 1.** Aftershock sequences a.

 a earthquake information came from USGS earthquake catalog: https://earthquake.usgs.gov/earthquakes/.

#### **3. 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.

#### **4. 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–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:

$$\chi^2 = (\mathbf{d}\_{\rm obs} - \mathbf{d}\_{\rm sim})\mathbf{C}^{-1}(\mathbf{d}\_{\rm obs} - \mathbf{d}\_{\rm sim})^{-1}$$

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.
