Logarithmic Model Joint Inversion Method for Coseismic and Postseismic Slip: Application to the 2017 Mw 7.3 Sarpol Zahāb Earthquake, Iran

Interferometric synthetic aperture radar (InSAR) has become an important technique for studying earthquake cycle deformation. However, due to the limited satellite revisit time, it is often difficult to fully separate coseismic and postseismic slip from InSAR data. Nevertheless, accurately estimating spatiotemporal coseismic and postseismic fault slip distribution models is important for quantifying earthquake slip, understanding the kinematics of seismogenic faulting, and evaluating seismic hazards. Here, we develop a logarithmic model‐based method (LogSIM) for the joint inversion of coseismic and postseismic fault slip using InSAR data from multiple platforms in different orbits. This method considers the nature of early postseismic slip following logarithmic decay. The coseismic slip, the decay time constant, and the decay amplitude of the logarithmic model can be jointly estimated from unwrapped interferograms without the need for InSAR time series calculations, thereby reducing the number of unknown parameters and stabilizing the inversion. The robustness of LogSIM is first validated with synthetic experiments. We then apply LogSIM to invert for the coseismic and postseismic slip models of the 2017 Mw 7.3 Sarpol Zahāb earthquake. We find that the maximum afterslip value and postseismic moment release during the first four months reported in previous studies are underestimated by ~26% and ~31%, respectively, due to the unavailability of the first five days of postseismic deformation. The estimated afterslip distribution spatially overlaps with the aftershocks in the updip region of the fault plane. LogSIM could potentially be extended to integrate different space geodetic data in earthquake cycle deformation modeling.


Introduction
Interferometric synthetic aperture radar (InSAR) data have been successfully used to study earthquake cycle deformation, especially when ground-based data are not available (Elliott et al., 2013;Fialko et al., 2005;Fielding et al., 2009;Li et al., 2009;Peltzer, 1999;Raucoules et al., 2010;Xu et al., 2016;Zhou et al., 2018). InSAR data can be inverted to derive finite fault slip models, which are important for understanding of fault geometry, fault slip, and fault frictional properties, as well as for assessing potential regional seismic hazards (Elliott et al., 2016;Marone, 1998;. However, due to the limited revisit time of radar satellites, coseismic deformation maps generally contain a few days or months of early postseismic deformation, which makes it difficult to accurately estimate a coseismic slip model (Feng et al., 2014;Heki & Tamura, 1997;Hudnut et al., 1996;Langbein et al., 2006;Yoshioka, 2001). Heki and Tamura (1997) showed that the first 24 hrs of postseismic slip after the 1994 Mw 7.8 Sanriku-Haruka-Oki earthquake was comparable to~30% of the corresponding coseismic slip. This suggests that coseismic slip values can be significantly overestimated if early postseismic signals are included in coseismic deformation maps. The estimated coseismic moment release of the 2004 Mw 6.0 Parkfield earthquake estimated from seismic data was about 28% less than that estimated using InSAR data, which contain the first two days of the postseismic signal (Kim & Dreger, 2008). The coseismic slip model of the 2011 Mw 7.1 Van earthquake was significantly overestimated due to the existence of the first three days of the early postseismic signal (Feng et al., 2014). Similar discrepancies between coseismic slip and moment estimates are observed for other events: the 1994 Mw 6.7 Northridge thrust earthquake Wald et al., 1996), the 1999 Mw 7.6 Chi-Chi thrust earthquake (Ma et al., 2001;Wang et al., 2001;Yoshioka, 2001), and the 1996 Mw 5.9 normal faulting event in the Pumqu-Xainza Rift . In terms of postseismic slip inversions, Yano et al. (2014) found that the first day of postseismic moment release of~6.0 × 10 17 Nm after the 2009 Mw 6.3 L'Aquila earthquake is nearly equal to the moment of 6.5 × 10 17 Nm released from the second day to 60 days (Cheloni et al., 2010). Ragon et al. (2019) found that the peak early afterslip amplitude of the 2009 L'Aquila earthquake was up to three times greater than previously reported (Cheloni et al., 2014;D'Agostino et al., 2012). These studies highlight that the afterslip can be seriously underestimated, while the coseismic slip can be overestimated if direct measurements of the early postseismic deformation are absent.
To jointly estimate coseismic and postseismic slip models with geodetic data, Johanson et al. (2006) developed a linear joint inversion model to combine GPS and InSAR data with the aim of estimating the afterslip amplitudes and coseismic slip of the 2004 Mw 6.0 Parkfield earthquake by establishing a prior uniform afterslip decay time derived from the GPS data. However, the coseismic slip derived with this method still includes the afterslip on the first day due to the daily sampling of GPS data (Johanson et al., 2006). In addition, the prior decay time parameter is required in advance for the least squares solution. Subsequently, Johnson et al. (2009) proposed a semilinear joint model to invert the coseismic and postseismic slip and other parameters, such as the asthenospheric viscosity, lithospheric thickness, and frictional parameters, associated with the 2002 Mw 7.9 Denali earthquake in Alaska. To avoid estimating the coseismic slip contaminated by the early postseismic signals, Yano et al. (2014) proposed a strategy to jointly invert coseismic and postseismic slip models simultaneously with seismic and geodetic data. They found substantial afterslip (~20% of coseismic slip) between 10 s and one day after the 2009 L'Aquila earthquake. However, this method relies on both seismic and geodetic data, and the number of parameters in this method will increase significantly with extended temporal afterslip (Yano et al., 2014). Ragon et al. (2019) modified this strategy by considering model uncertainties with the Bayesian method. They found an inverted peak coseismic slip was 30% higher than that reported in previous studies (Cheloni et al., 2014;D'Agostino et al., 2012). The cumulative early afterslip was significantly underestimated when the postseismic signals from the first 6 days were not included in previous studies (Cheloni et al., 2014;D'Agostino et al., 2012).
In contrast to the existing data-driven coseismic and postseismic slip inversion methods, we propose a novel logarithmic model-based coseismic and postseismic slip joint inversion algorithm (LogSIM). Compared with the existing inversion approaches, our proposed method can be used to jointly invert for coseismic slip and early afterslip by fully integrating different InSAR data sets. This paper is organized as follows. Section 2 illustrates the rationale of the proposed method that links individual unwrapped interferograms computed from separate synthetic aperture radar (SAR) data subsets in different orbits to the logarithmic model. The proposed method is validated by synthetic data sets, and the uncertainties of the parameters are presented in section 3. Section 4 is dedicated to present the obtained coseismic and postseismic slip models of the 2017 Mw 7.3 Sarpol Zahāb earthquake, Iran. In section 5, we discuss the presented algorithm and the results. Finally, conclusions and further developments are addressed in section 6.

Methodology
Let us consider that a number of SAR images m 1 , m 2 , ⋯ , m N obtained from N orbits that are available in a study area. M i (i = 1,2,…,N) interferograms are generated based on given spatial and temporal baseline thresholds satisfying the following inequality (Berardino et al., 2002): These interferograms might be generated using the traditional two-pass differential InSAR data processing method (Wegmuller, 1997). An unwrapped interferogram typically contains millions of pixels, which makes the direct inversion of data difficult. To conserve computational resources and time without losing details of the deformation, we use the quadtree method to downsample the data (Jónsson et al., 2002). The fault geometry is usually estimated from coseismic surface displacements based on a rectangular dislocation model (Okada, 1985) using nonlinear optimization search methods (Cervelli et al., 2001). Alternatively, the fault geometry can be inferred from previous geological or seismic studies in the same study area. For large thrust events, the global subduction geometry might be used as the fault geometry (Hayes et al., 2018).
The fault plane can be divided into P rectangular patches. The relationship between the observed surface displacement vector D from multiple data sets and the unknown spatiotemporal fault slip vector X is described as follows: where G is the matrix of Green's functions constructed using the elastic dislocation model (Okada, 1985) that relates the unit slip along the rake direction on each fault patch to surface displacements, B is a design matrix connecting G and X, T is the total postseismic number of SAR images, and x i co and x i;j post i ¼ 1; 2; ⋯; P; j ¼ 1; 2; ⋯; T ð Þ are the coseismic slip vectors and cumulative afterslip vectors on the fault patch i at each SAR image acquisition epoch j, respectively.d ti coþpost i ¼ 1; 2; …; T represents the interferograms that include both coseismic and postseismic deformation, d ti post;tj i ¼ 1; 2; ⋯; T; j ¼ 0; 1; 2; ⋯; T−1 ð Þ represents the postseismic displacement induced by afterslip between SAR image acquisition epochs i and j, and ε is the model misfit.
To overcome the rank deficiency problem of the above functional model (2), we introduce the following logarithmic decay function for each fault patch: where A i is the decay amplitude of fault patch i, τ is a uniform time decay constant, t j is the acquisition date of SAR image j, and t 0 is the onset time of the earthquake. We assumed that the afterslip decay time constant τ is uniform for all fault patches to reduce the number of unknown parameters and improve the computational efficiency (Gonzalez-Ortega et al., 2014;Johanson et al., 2006). Therefore, by combining equations (2)  6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 ½ =τ ½ 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 The number of unknown parameters in equation (4) is 2P+3 (including 2P, which denotes the number of coseismic slip patches x i co and afterslip decay amplitudes A i , one uniform decay time constant τ, one coseismic slip rake angle and one uniform afterslip rake angle). The number of unknown parameters in the proposed method depends only on the number of patches on the fault plane. Therefore, the number of unknown parameters to be estimated is much less than that in the traditional inversion method (P×(T+1)), as shown in equation (2). We use the nonlinear least squares curve fitting function (NLSCF) routine, a built-in function of MATLAB software, to solve the nonlinear equation (4) and search for the global optimal solutions (Levenberg, 1944;Marquardt, 1963). A flowchart diagram illustrating the main steps of the proposed method is demonstrated in Figure 1.

Data Simulation
We simulate observations from both ascending and descending orbits with fundamental parameters similar to C-band Sentinel-1 SAR data ( Figure 2). The ascending data include one preseismic epoch and five postseismic epochs, while the descending data include only five postseismic epochs. In total, 15 ascending and 10 descending unwrapped interferograms are generated. These interferograms are downsampled by using quadtree algorithm (Jónsson et al., 2002), after which random Gaussian white noise is added with a certain standard deviation reaching up to 5% of the maximum deformation. These interferograms are jointly inverted for the coseismic and spatiotemporal postseismic slips in the experiment.
We choose a strike angle of 351°and a dip angle of 15°to mimic the seismogenic fault plane of the 2017 Mw 7.3 Sarpol Zahāb earthquake . The coseismic and postseismic slip rake angles are set to 150°and 120°, respectively. The fault plane is divided into 120 patches with a patch size of 5 km × 5 km. The synthetic coseismic and postseismic slip distributions are assumed to take adjacent circular crack shapes with maximum values of 10 and 2 m at depths of 10 and 6.5 km, respectively. The synthetic coseismic slip and the cumulative afterslip on each fault patch are simulated using the logarithmic decay function, in which the decay constant τ is set to be 10 days ( Figure 4). We also carry out another synthetic experiment with a finer fault plane, and the estimated fault slips show good agreement with the inputs ( Figure S1 in the supporting information), but the experiment with a finer fault plane is more time consuming than the original experiment.

Estimating the Uncertainties of the Parameters
LogSIM uses the NLSCF routine to solve the optimal parameters X by minimizing the L2-norm of the misfit between the model D ′ = GBX and the data D:

10.1029/2019JB017953
Journal of Geophysical Research: Solid Earth where W is a weight vector. We set a searching boundary for each parameter according to prior knowledge (Tables S1 and S2).
The NLSCF routine uses the Levenberg-Marquardt method (Levenberg, 1944;Marquardt, 1963) to compute updates ΔX to a vector of initial parameter X 0 by solving the following linear equation (6): provided that there are N1 InSAR observations and N2 parameters, where J(X) is the N1×N2 Jacobian matrix of |WD ′ − WD|, I is an identity matrix, and α is a damping factor controlling the magnitude and direction of ΔX. The fitting iteratively updates the values of the parameters until one of several stopping criteria (e.g., the maximum number of iterations or function tolerance) is met.
Upon the successful completion of the NLSCF routine, the output results are the optimal values of the parameters, the residuals R N1 × 1 and the Jacobian matrix J N1 × N2 . Then, the covariance matrix of parameters C N1 × N2 can be calculated with the equation (7): The uncertainty of each parameter is calculated as the square root of its variance.

Accuracy Evaluation
As shown in Figure 4, the results estimated with LogSIM are in good agreement with the synthetic data. The estimated slip patterns match well with the inputs although no spatial or temporal smoothing was applied in the joint inversion. The root-mean-square (RMS) of the residuals between the synthetic coseismic slip and the estimated results is~0.9 m, and the max RMS of the residuals between the synthetic spatiotemporal postseismic slips and their corresponding estimated results is~0.07 m. In addition, the residual histograms shown in Figures 3p-3r follow a standard normal distribution with small variances.
The estimated decay time constant, coseismic, and postseismic rake angles are 10.6 ± 0.9 days, 150.1°± 0.4°, and 119.4°± 0.8°, respectively, which match well with the input values. The correlation coefficient between the input and estimated spatial and temporal fault slips on 100 fault patches is 0.92 ( Figure 5). Therefore, the synthetic experiment demonstrates that LogSIM is robust in jointly estimating coseismic and spatiotemporal postseismic slips from multiple InSAR data sets.
We designed further synthetic experiments to test the sensitivity of LogSIM to different levels of data noise, fault geometry, the first postseismic acquisition times and densities of postseismic data (Tables S3 and S4). The estimated rake angles, afterslip decay amplitude and decay time constant are independent of the fault geometry in the synthetic experiments (Table S3). Furthermore, we find that the density of postseismic data does not play an important role in the afterslip estimation, and the estimated coseismic slip and afterslip become less accurate if the early postseismic deformation data are absent (Table S4). This indicates the importance of early postseismic data in the inversion of fault slip.   earthquake to date in the Zagros Mountains (Barnhart et al., 2018;Chen et al., 2018). Geologically, the earthquake occurred within the fold-thrust belt of the northern Zagros Mountains, one of the most active orogens in the world ( Figure 6). The northern Zagros Mountains have been uplifting since the late Eocene due to the convergence between the Eurasian and Arabian continents with a current convergent velocity of~3 cm/year (DeMets et al., 2010;Mouthereau et al., 2012). Previous studies revealed that there is a set of thrust faults distributed across the northern Zagros Mountains due to the tectonic collision between Eurasia and Arabia (i.e., NE-SW convergence; Barnhart et al., 2018;Nissen et al., 2019).

InSAR Data Processing and Model Setup
The coseismic and postseismic ground deformation of the 2017 Sarpol Zahāb earthquake were imaged by the C-band Sentinel-1 satellite and the L-band ALOS-2 satellite in six different orbits (see Chen et al., 2018, for a detailed description of the data). The major steps of InSAR data processing are as follows: coregistration of slave images to the corresponding selected master images in different orbits, generation of raw interferograms, removal of the flat-earth phase effects and the topographic phase from the interferograms, interferogram filtering and phase unwrapping, conversion of the unwrapped phase to the radar line-of-sight displacement, and geocoding to World Geodetic System 1984 coordinates. Considering that unwrapped interferograms are often contaminated by atmospheric noise, we applied the Generic Atmospheric Correction Online Service for InSAR tools to correct the atmospheric effects on the interferograms (Yu et al., 2018). A total of 497 interferograms, including four coseismic and 493 postseismic displacement fields, was selected for the joint estimation of coseismic slip and spatiotemporal afterslip (Table 1). It should be noted that the Sentinel-1 SAR data acquired on ascending track 72 and descending track 79 and two ALOS-2 data cover both the coseismic deformation and the first several days of postseismic deformation (Figure 7). The seismogenic fault is set to be planar with a dip angle of 15°and a strike angle of 351°according to a previous study . The top of the fault is extended to the surface. The dimensions of the fault are assigned to be 120 km × 120 km (Figure 6), and the fault is divided into 15 × 15 fault patches. There is still no best way to automatically determine optimal weights of different InSAR data sets; we used the inverse of data variance (see Table S5) as the relative weights in the inversion. Following the approach described in section 2, we estimated the coseismic slip and spatiotemporal afterslip of the earthquake simultaneously using the downsampled Sentinel-1 and ALOS-2 InSAR data from six different tracks.

Results
Significant coseismic ground deformation signals caused by the mainshock were captured by both Sentinel-1 and ALOS-2 satellites (Figures 8a-8d). The maximum coseismic surface displacement of~0.8 m is observed Note. SAR = synthetic aperture radar.

10.1029/2019JB017953
from the ascending orbit. In contrast, only up to~0.5 m of coseismic ground deformation is observed from the descending orbit because of its less favorable viewing geometry. The coseismic displacements show apparent symmetric double-lobe patterns, including the southwestern region characterized by uplift and the northeastern region characterized by subsidence; these patterns, which were captured by both ALOS-2 and Sentinel-1 data, imply a blind NW-SE oriented seismogenic fault system. Four postseismic displacement maps spanning the first 120-day period from two ascending and two descending orbits of Sentinel-1 are shown in Figure 9. The postseismic deformation patterns are similar to the coseismic ground deformation maps but are shifted toward the SW. Up to 8 cm of postseismic ground movement toward the satellite is observed in the ascending data (Figure 9a). The forward model prediction results effectively reproduce the main features in the observations (Figures 8 and 9). The RMS is 3.0 and 2.4 cm for the Sentinel-1 ascending and descending orbits, respectively. The RMS is 2.7 and 3.2 cm for the ALOS-2 coseismic residual displacement fields, respectively. The RMS is 1.1, 0.7, 0.8, and 0.6 cm for four postseismic residual displacement fields. The histograms of the residuals (Figures 8m-8p and Figures 9m-9p) roughly follow normal distributions with small variations.
The coseismic slip and spatiotemporal afterslip on the seismogenic fault were estimated using LogSIM ( Figure 10). The maximum coseismic slip of~8 m is located close to the hypocenter at a depth of~16 km. Two major slip asperities were estimated by the proposed method near the epicenter, which agrees well with other studies (Barnhart et al., 2018;Chen et al., 2018;Feng et al., 2018). The kinematic spatiotemporal evolution of postseismic slip shows that the afterslip migrated updip of the seismic rupture ( Figure 10). The maximum afterslip of~0.5 m is localized between 11 and 14 km at the top edge of the coseismic slip zone. Afterslip is observed updip of the coseismic slip zone, which might reflect the difference between the frictional and material properties in and around the coseismic ruptures.
The temporal fault slip evolution of four selected fault patches ( Figure 10) is shown in Figure 11. These logarithmic decay afterslip patterns are consistent with the rate-and-state friction law (Dieterich, 1979;Ruina, 1983) applied widely in studies of the entire earthquake cycle. The coseismic rupture of the Sarpol Zahāb earthquake caused~1.8 m of movement on fault patch 100 (marked by the blue lines in Figure 11), which was followed by the rapid decay of afterslip reaching a cumulative value of~0.43 m within 120 days after the mainshock (marked by the black lines and red dots with one standard deviation error bar in Figure 11).

Performance of LogSIM
Although no regularization constraints were applied during the estimation of the spatial distribution of fault slip, LogSIM is effective at capturing the relatively smooth distributions of coseismic and spatiotemporal postseismic fault slip ( Figure 10). LogSIM is able to estimate slip in the very early period of postseismic deformation by combining data acquired in different orbits ( Figure 11). The temporal resolution of the estimated afterslip can be increased by up to one day using LogSIM (Figures 10 and 11). Postseismic deformation follows two major phases of deformation: an initial short-term initiation of afterslip that lasts for approximately a few days and a subsequent phase consisting of relatively long-term steady-state relaxation (Figures 10 and  11). The first initiation phase is critical for distinguishing whether a rate-and-state friction law or ratedependent friction law is preferred to explain the postseismic surface displacement. It is also important to assess the risk for future large aftershocks (Twardzik et al., 2019).
Overall, the spatial distribution of the estimated coseismic and postseismic fault slip in our study is in good agreement with that estimated in other studies ( Table 2). The coseismic moment release estimated to be 0.97 × 10 20 Nm (corresponding to Mw 7.3) is consistent with other studies (Barnhart et al., 2018;Chen et al., 2018;Ding et al., 2018;Feng et al., 2018;Kobayashi et al., 2017;Kuang et al., 2018;Nissen et al., 2019;Vajedian et al., 2018;Wang et al., 2018). The first 120-day postseismic moment release of 1.24 × 10 19 Nm is only~13% of the coseismic moment. If the first five days of data are not considered (Barnhart et al., 2018;Feng et al., 2018), the estimated postseismic moment released reported in the existing studies is 0.9 × 10 19 Nm, which is nearly 31% smaller than our estimation. Similarly, early postseismic signals can explain the discrepancy between our estimated maximum postseismic slip of~0.53 m and the slip of 0.4 m estimated from Barnhart et al. (2018). Our estimated maximum coseismic slip and slip patterns also show some differences from the existing studies (Table 2). Our maximum coseismic slip of~8 m is the largest value, suggesting that other studies may underestimate the maximum slip due to introduced smoothness constraints (Dettmer et al., 2014). In addition, this difference may also be explained by differences among the inversion strategies, input data and fault geometry in different studies. Therefore, compared with the existing fault slip inversion methods, LogSIM has three advantages: (1) LogSIM can be used to jointly invert for coseismic slip and early afterslip by considering the nature of postseismic slip following a logarithmic decay function. (2) LogSIM can improve the temporal resolution of the estimated afterslip by fully integrating different InSAR data sets. (3) Spatial regularization might not be required in the inversion if the temporal evolution of afterslip is well defined.

The Predominant Process of Postseismic Deformation
Aseismic afterslip, poroelastic rebound, and viscoelastic relaxation are three fundamental mechanisms responsible for postseismic transients (Zhao et al., 2017). Many existing studies show that afterslip often dominates early postseismic deformation spanning from several days to months (Ding et al., 2015; Segall Figure 10. The estimated coseismic slip and cumulative spatiotemporal afterslip for each time step using logarithmic model-based method. Black contour with a 1.5-m interval represents coseismic slip on the seismogenic fault. Star represents the hypocenter determined by the Iranian Seismological Center. The detailed fault slip evolution of the selected four patches marked by arrows in the last four subpanels is shown in Figure 11. Shrivastava et al., 2016). The spatiotemporal afterslip estimated in this study effectively explains the observed early postseismic transient deformation. Therefore, it is suggested that afterslip is the most likely driving process causing the observed short-term transient deformation.
A previous seismic imaging study reported that the average crustal thickness beneath the Main Recent Fault in the Northern Zagros Mountains is 56 ± 2 km (Paul et al., 2010), indicating that the seismogenic fault of the 2017 event is located in the brittle upper crust. The effective elastic plate thickness is about 60 km in the northwestern Zagros Mountains range (Audet & Bürgmann, 2011), and the geothermal gradient is relatively low (Saein, 2018), which suggests a rather cold/strong crust and lithosphere. Barnhart et al. (2018) found that postseismic relaxation models employing both Newtonian rheology and power law rheology with viscosities in the lower crust and upper mantle ranging from 10 21 to 10 18 Pa s predict surface deformation on the order of 10 −6 to 10 −3 m surface deformation, respectively, which is much lower than theobserved early postseismic deformation in this study. Therefore, we rule out viscoelastic relaxation as a significant contributor to the early postseismic transient deformation of the 2017 event.
Although poroelastic rebound may have occurred due to ubiquitous fluids in the upper crust (Mahsas et al., 2008), it is rare to find cases of poroelastic rebound dominating early postseismic deformation. Poroelastic deformation was previously modeled by using different Poisson's ratios for undrained and drained crustal rocks in a homogeneous half-space , in which the pattern and magnitude of the estimated poroelastic deformation were found to be inconsistent with the InSAR observations. Therefore, the contribution of poroelastic rebound seems to be negligible.

Relationship Between Afterslip and Aftershocks
There appeared to be a close relationship between seismic aftershocks and aseismic afterslip following the 2017 Mw 7.3 Sarpol Zahāb earthquake. More than 1000 Mw > 2.5 aftershocks have been documented by the Iranian Seismological Center. As shown in Figure 12, the cumulative number of aftershocks suggests a logarithmic increase indicating the rapid release of energy after the mainshock. Similarly, the cumulative moment calculated from the estimated spatiotemporal afterslip is also characterized by a logarithmic increase, demonstrating an aseismic adjustment after the earthquake. These similar patterns imply that rapidly decaying aseismic afterslip may have triggered the aftershocks (Helmstetter & Shaw, 2009). Figure 11. The temporal fault slip evolution of selected fault patches (see the last four subpanels in Figure 10). The blue y axis represents the coseismic slip, and the orange y axis represent the accumulated afterslip with uncertainty of one standard deviation. The uncertainties of temporal afterslip are propagated from the uncertainties of parameters according to the law of error propagation based on equation (3). a − b represents the frictional parameter that is discussed in the section 5.4. Moreover, the moment released by updip aftershocks within the study period is approximately 4.38 × 10 17 Nm, equivalent to Mw 5.73, which is only 3.5% of that released by afterslip in this study, indicating that most of postseismic deformation occurs aseismically (Burgmann et al., 2002). In addition, the inverted accumulative afterslip agrees well with the spatial distribution of the aftershocks on the updip fault plane (see Figure 1). We suggest that most of the aftershocks were driven by aseismic afterslip (Avouac, 2015;Barbot et al., 2009;Hsu et al., 2006;Perfettini & Avouac, 2007). The logarithmic patterns of afterslip and aftershocks with small decay time constants indicate that the postseismic activities decreased rapidly, possibly due to the fast healing of the fault system (Bedford et al., 2016). Some  aftershocks occurred on the downdip fault plane, but our inversion results and the existing studies do not show obvious postseismic fault slip at such depths (Barnhart et al., 2018;Feng et al., 2018). This probably suggests that InSAR data have a poorer resolution and inversion stability for downdip fault slip.

Frictional Properties of the Seismogenic Fault
The well-documented logarithmically decaying afterslip ( Figure 11) is similar to that derived from laboratory rock tests and other earthquake studies based on the rate-and-state friction law (Dieterich, 1979;Marone, 1998;Marone et al., 1991;Ruina, 1983;Scholz, 2002;Wennerberg & Sharp, 1997). The logarithmic model S(t) = Aln(1+t/τ) used in our study has the same form as S t Scholz, 2002), where S(t) denotes the afterslip evolution with time t after the mainshock, σ is the normal force, k is the crustal stiffness, v max is the maximum slip rate, and a − b is the difference between the rate and state parameters that determines whether part of a fault is a stable velocitystrengthening zone (positive a − b values) or an unstable velocityweakening zone (negative a − b values; Scholz, 2002). If the crustal stiffness and the normal force are known, the spatially varying parameter a − b for each fault patch can be inferred from the corresponding amplitude A and decay time constant τ. Assuming zero fluid pressure, the normal force can be estimated by σ = ρgh × cos(dip) (Zhou et al., 2018), where ρ is the crustal density (2.7 g/cm 3 in this study), g is 9.8 m/s 2 , and h and dip are the depth and dip angle of arbitrary fault patch, respectively. The crustal stiffness can be calculated by k = G/H (Marone et al., 1991), where H is the constant thickness of a velocity-strengthening zone (15 km is used in this study). Hence, the spatial parameter a − b can be calculated by a − b = Ak/σ. Our estimated values of a − b ranging from 10 −4 to~2 × 10 −3 agrees well with those of other studies in different tectonic settings and laboratory experiments ( Figure 13 and Table 3). Hsu et al. (2006) found frictional afterslip after the 2005 Mw 8.7 Nias-Simeulue megathrust earthquake in Sumatra with an a − b value of about 5 × 10 −4 , which is one to two orders of magnitude lower than that estimated in the laboratory. The high pore pressure at the subduction zone interface after an earthquake may explain this difference. Barbot et al. (2009) reported that the postseismic deformation after the 2004 Mw 6.0 Parkfield earthquake can be best explained by a rate-strengthening model with a − b of about 7 × 10 −3 near the upper limit value of laboratory experiments. Zhou et al. (2018) obtained a rate-and-state friction parameter a − b of about 3 × 10 −3 based on time-dependent postseismic slip following the 1978 Mw 7.3 Tabas-e-Golshan earthquake.
A low friction parameter a − b may explain why both coseismic slip and afterslip has occurred in the overlap region (see Figures 1 and 13), where a − b takes values on the order of 0.0005, allowing dynamic ruptures to propagate into velocity-strengthening regions (Kanu & Johnson, 2011;Zhou et al., 2018). Afterslip below the seismic region in this study is expected to occur in steady-state velocity-strengthening regions controlled by frictional properties. The variations in the spatial frictional properties suggest that the physical or  Sarpol Zahāb Thrust lithological properties of the fault rocks vary with depth (Barbot et al., 2009;Hsu et al., 2006;Miyazaki et al., 2004;Song & Simons, 2003). Laboratory data suggest that a − b decreases with depth in the upper 5 km or so (Blanpied et al., 1998), which reconciles with our estimates (see Figure 13a) indicating that the frictional parameter changes gradually from approximately 10 −4 to~2 × 10 −3 .
Interestingly, the afterslip of the 2017 event was primarily concentrated updip of the coseismic slip zone, while the postseismic deformation following other large earthquakes in similar orogenic belts thrust events such as the 1999 Mw 7.3 Chi-Chi earthquake (Hsu et al., 2002) and 2015 Mw 7.9 Gorkha earthquake (Zhao et al., 2017) showed that afterslip mainly occurred downdip of the coseismic rupture zone. According to the rate-and-state friction law (Dieterich, 1979;Scholz, 2002), the absence of afterslip downdip of the coseismic rupture may indicate velocity-weakening behavior in the deep crust, where the thick-skinned tectonics of the crystalline basement hosted the Sarpol Zahāb earthquake (Barnhart et al., 2018) and responded to coseismic slip by a rapid decrease in frictional stress while remaining strongly coupled in the lower crust (Mahsas et al., 2008). In contrast, a thin-skinned tectonic model was found to explain the behavior of the 1999 Chi-Chi earthquake (Hsu et al., 2002;Wang et al., 2000). The mainshock of the Chi-Chi earthquake occurred on a shallow-dipping (24°) ramp fault, and aftershocks were observed on the lower decollement with an shallower dip angle (6-10°; Hsu et al., 2002). Geologically, the Zagros Mountains are covered by 8 to 13 km thick Phanerozoic sedimentary strata folded into the characteristic anticlines (Barnhart et al., 2018). The updip migration of afterslip following the Sarpol Zahāb earthquake can be described by shallow thinskinned processes enhanced by increased stress after the mainshock. The depth (0-15 km) of afterslip is consistent with the depth of the base of the sedimentary units within the Zagros Mountains (Casciello et al., 2009). Therefore, the afterslip found in this study indicates low frictional properties in the shallow crust, where the folding of sediments may modulate part of the coseismic stress disturbance (Donnellan & Lyzenga, 1998;Mahsas et al., 2008). This stable thin-skinned process is in accordance with the velocitystrengthening behavior found in this study. In addition, Barnhart et al. (2018) found that a weak basal decollement with relatively low frictional resistance beneath the flat potion of the Mountain Front Fault, along which the residual strain after the mainshock was allowed to extend aseismically upward; in contrast, this kind of decollement was not observed for the Gorkha earthquake. Thus, we suggest that the location of afterslip with respect to the corresponding coseismic rupture zone depends mainly on the structure and frictional properties of the fault.

Conclusions
The proposed LogSIM joint inversion method makes use of all available interferograms to improve the temporal resolution of the estimated fault slip history. The results of synthetic experiments showed that LogSIM can efficiently and accurately recover the spatiotemporal fault slip, especially when the signal-to-noise level of the data is high. We applied the method to the real data from the Mw 7.3 Sarpol Zahāb earthquake that occurred along the border zone between Iran and Iraq in 2017. We found that the maximum coseismic slip reached 8 m and the first four months of afterslip reached~0.5 m. The inverted afterslip occurred primarily in the updip region of the coseismic rupture zone, which suggests that afterslip is an important mechanism of postseismic deformation and can trigger aftershocks. Based on the inverted afterslip, we infer that the friction parameter a − b ranges from 10 −4 at depth to 10 −3 close to the surface, which is consistent with the findings derived from laboratory experiments and other earthquakes. The use of LogSIM will benefit studies of earthquake rupturing, fault slip evolution, frictional heterogeneity, seismic cycle budget, and future earthquake risk assessments, especially for studies using only InSAR data due to their coarse temporal resolution.
Although we applied LogSIM to only InSAR data, the proposed method can be easily implemented for other geodetic observations. Challenges facing the application of our method include trade-offs between the computational efficiency and the number of discrete fault patches, the roles of other potential postseismic mechanisms, such as poroelastic rebound and viscoelastic relaxation, and the effects of larger aftershocks on the spatiotemporal evolution of afterslip.