Determination of the dipping direction of a blind reverse fault from InSAR: case study on the 2017 Sefid Sang earthquake, northeastern Iran

Determining the fault parameters of an earthquake is fundamental for studying the earthquake physics, understanding the seismotectonics of the region, and forecasting future earthquake activities in the surrounding area. Dense crustal deformation data such as Interferometric Synthetic Aperture Radar (InSAR) are useful for fault parameter determination, but determining the dipping direction of a blind fault is often challenging when the size of the earthquake is not large (M < 7) or when the coverage of the deformed area is limited to capture the details of rupture. The 5th April 2017, Mw 6.1 earthquake occurred near the city of Sefid Sang, northeast of Iran, provides an excellent case for exploring the potential of InSAR data for determining the dipping direction of a blind reverse fault. Using Advanced Land Observing Satellite-2 (ALOS-2) and Sentinel-1A interferograms of four different observation directions and a fault slip inversion method that allows thorough exploration of the fault geometry led to two candidates of reverse fault models, dipping either to the northeast or the south. The results show that the fault models of both dipping directions explain the data well, with a slight advantage in the northeast-dipping fault model in terms of the misfit when the atmospheric corrections were applied. The northeast-dipping fault model is, in addition, more consistent with the strike of the mapped active faults in the region and with the aftershock distribution, from which we infer that the 2017 Sefid Sang earthquake occurred on a northeast-dipping dextral-reverse fault. The preferred fault model has a strike angle of 314.8°, dip angle of 47.4° and rake angle of 130.3°, and a slip distribution of maximum 1.35 m at depth of 5 km equivalent to Mw 6.0. This study illuminates the difficulty of determining the dipping direction of blind faults even with InSAR measurements from multiple directions, but also that correcting for the atmospheric noise and comparing with other kinds of data can help infer the fault dipping direction.


Introduction
On 5 April 2017, an Mw 6.1 earthquake occurred in the northeastern region of Iran, about 90 km southeast of the city of Mashhad, the second most populated city in Iran with a population of over 3 million people (Fig. 1). The epicenter (35.776 N, 60.436 E, United States Geological Survey (USGS)) was at a remote mountainous area, about 30 km northeast of the Sefid Sang city where approximately 6000 people lived. This earthquake occurred in a region where no major fault with significant displacement has been recognized (Aflaki et al. 2019).
In northeastern Iran, the Kopeh Dagh and Binalud Mountains are the main deformation belts situated at the northeast boundary of the collision zone between Ghayournajarkar and Fukushima Earth, Planets and Space (2020) 72:64 Arabian and Eurasian plates (Fig. 1). The plate convergence rate is approximately 25 mm/year in the direction N10° E (Mcclusky et al. 2003;Vernant et al. 2004). The source mechanism solutions of the 2017 Sefid Sang earthquake show that the earthquake ruptured either a northwest-or an east-striking oblique reverse fault consistent with the regional transpressional stress regime (Shabanian et al. 2010).
There remains a question as to the dipping direction of the fault. Indeed, the difficulty in determination of the fault dipping direction for blind strike-slip, reverse and normal faults has been known (e.g., Biggs et al. 2006;Gardonio et al. 2018;Gong 2019;Materna et al. 2019). Xu et al. (2018) used Sentinel-1A SAR data and conducted two sets of inversions for different dipping directions. Based on the degree of misfit, they inferred that the 2017 Sefid Sang earthquake occurred on a southwest-dipping fault plane. Su et al. (2019) used the Sentinel-1A SAR data and conducted a set of inversion assuming northeast-dipping fault plane as the seismogenic fault. Based on the consistency with the Iranian Seismological Center's (IRSC) focal-mechanism solution, they concluded that their model was acceptable. Aflaki et al. (2019) also conducted a set of inversion assuming northeast-dipping fault plane using seismic and geodetic (Sentinel-1A interferograms) data. Based on the northwest-trending hypocentral distribution of aftershocks and northeastward increase in their focal depths, they preferred the Location of the 2017 Sefid Sang earthquake in northeastern Iran. Colored dots represent the aftershocks with magnitude larger than 2.5 located by the Iranian Seismological Center (IRSC) between 5 April 2017 and 31 December 2018, plotted together with the focal mechanisms of the mainshock determined by different organizations. The red lines represent faults from the USGS World Energy Project (https ://catal og.data.gov/datas et/major -fault s-in-iran-flt2c g). The cities of Mashhad and Sefid Sang are also indicated. Black boxes outline the used ascending and descending ALOS-2 SAR frames and white boxes outline the used ascending and descending Sentinel-1A SAR frames Ghayournajarkar and Fukushima Earth, Planets and Space (2020) 72:64 northeast-dipping fault plane with dextral-reverse motion for the Sefid Sang mainshock. All these three studies adopted a two-step approach for the inversion, first determining the fault geometry with a simplified slip assumption and then solving for the slip distribution on the obtained geometry, which is not an ideal procedure to search for the optimal geometry. In this study, we attempt to clarify the dipping direction of the seismogenic fault using an additional data set, atmospheric correction model, and a different modeling approach compared to previous studies. We used not only the Sentinel-1A, but also the ALOS-2 interferograms. We applied the atmospheric noise correction obtained using the Generic Atmospheric Correction Online Service (GACOS) for InSAR (Yu et al. 2018), after checking the consistency between the Sentinel-1A and ALOS-2 interferograms. We explored the possibility of different dipping directions by conducting inversions assuming both northeast-and south-(or southwest-) dipping fault planes by applying a Bayesian inversion method (Fukushima et al. 2013) that simultaneously solves for the fault geometrical parameters and slip distribution. Finally, we compared the candidate models with the active fault traces in the region and the aftershock distribution, and inferred the dipping direction of the ruptured fault.

InSAR observations
We used two pairs of ALOS-2 images acquired in Strip Map (SM) mode on the ascending track 171 and descending track 64, respectively, and two pairs of Sentinel-1A images acquired in Terrain Observation with Progressive Scan (TOPS) mode on the ascending track 13 and descending track 93, respectively (Table 1, Fig. 1).
All SAR data were processed using the GMTSAR software (Sandwell et al. 2011). The effects of topography were removed using the 90-m Shuttle Radar Topography Mission digital elevation model (SRTM DEM) (Farr et al. 2000). The interferograms were then filtered to suppress the phase noise with a power spectrum filter (Goldstein and Werner 1998), phase unwrapped using SNAPHU package (Chen and Zebker 2000) and geocoded to the geographic coordinates. The interferograms were corrected for atmospheric delay using GACOS. The residual orbital errors were then modeled using a polynomial function and removed from the interferograms.
The resulting interferograms all show a similar NW-SE elongated bull's-eye patterns (Figs. 2 and 3). The ascending and descending interferograms with atmospheric corrections show the maximum line-of-sight (LOS) coseismic displacements of 10.7 and 14.2 cm, respectively, for the ALOS-2 SAR data, and 10.3 and 13.9 cm, respectively, for the Sentinel-1A SAR data.
By comparing the interferograms without ( Fig. 2) and with (Fig. 3) the atmospheric corrections, we notice that the fringe patterns in and to the north and south of the bull's-eye patterns of the descending interferograms were modified by the corrections (Figs. 2b,d and 3b,d). The atmospheric corrections using weather models are not always effective especially for wavelengths shorter than ~ 75 km (Murray et al. 2019), but the consistency in the corrected descending interferograms (Fig. 3b,d) indicates the effectiveness of the corrections for this particular case. As we shall see in the model interferograms (Figs. 8 and 9), the displacements in these areas are important for constraining the fault dipping direction. In the ascending interferograms (Figs. 2a, c and 3a, c), the atmospheric corrections do not seem to be effective. The inversion results presented later in this paper were basically obtained using the data with the atmospheric corrections.

Preparation of data
The data used in the inversions were prepared as follows. First, we transformed the atmosphere-corrected interferograms ( Fig. 3) into local kilometric Cartesian coordinates with the center of the deformed area set as the origin, and subsampled the data points in concentric grids (e.g., Fukushima et al. 2013). This subsampling method has been shown to be effective by a previous study (Fukushima et al. 2005) that compared with other subsampling methods such as the quadtree algorithm (Jonsson et al. 2002). Subsampling was conducted in such a way that the interval between the data points in the central area of deformation was set to be 1 km, smaller than spatial resolution of the slip distribution we aimed to solve for. The density of the data points became lower at areas farther from the deformation center. We set the areas used for inversions wide enough (Fig. 4). Areas with decorrelated signals were not used. The subsampling resulted in 639 and 825 data points for the ascending and descending pairs of ALOS-2 SAR data, respectively, and 926 and 926 data points for the ascending and descending pairs of Sentinel-1A SAR data, respectively (Fig. 4).

Model assumptions and parameters
We assumed the Earth to be a homogeneous and elastic half-space, with Poisson's ratio of 0.25. The top of the fault was assumed to reach the ground surface. The length (30 km) and the width (20 km) of the fault were set large enough that the slip can be modeled with zero displacement close to the assumed fault edges. Slip was allowed at the top of the fault. The fault was discretized with triangular elements in such a way the size of each triangular element was 2 km. We solved for both nonlinear and linear parameters. The nonlinear parameters consisted of the dip angle, rake angle, strike angle, horizontal location of the fault top side and weight of the smoothing constraint. As for the rewrapped with an interval of 3 cm. a, b ALOS-2 interferograms from ascending track 171 and descending track 64, respectively. c, d Sentinel-1A interferograms from ascending track 13 and descending track 93, respectively. The white and black rectangles represent the northeast-and south-dipping fault models, respectively, determined by the inversions in this study. The black arrows represent satellite heading (azimuth) and LOS directions Ghayournajarkar and Fukushima Earth, Planets and Space (2020) 72:64 horizontal location of the fault top side, we did not need two parameters because moving the fault parallel to the fault strike would lead to practically an identical result if the fault is set long enough. Therefore, we only introduced the location of the mid-point of the fault top side along the northing axis as the location parameter. We fixed the location of the mid-point of the fault top side along the easting axis after preliminary inversions. The linear parameters were the amount of slip on the 280 triangular fault elements and the displacement offsets for each InSAR data set. The offsets were required because the InSAR data were not referenced to any point.
We conducted two sets of inversions assuming both northeast-and south-dipping fault planes, but wide ranges of the nonlinear parameters were set to explore all the possible fault geometries (Table 2). We assumed these two fault planes based on the moment tensor solutions (Table 3), indicating the source fault was either an oblique reverse fault striking to the northwest or to the east.

Likelihood function
We took the maximum likelihood model as the optimum model, where the likelihood function p was defined as follows (Fukuda and Johnson 2008): (1) as adopted by previous studies (Fukushima et al. 2013).
Here, C is an arbitrary constant, d is the data vector whose length is N, s is a model vector composed of linear parameters whose length is M (284 in our problem), m is a model vector composed of the parameters defining the geometry of fault (nonlinear model parameters except α 2 ), σ 2 is the variance of the data noise, α 2 is the weight of the smoothing matrix L(m), G(m) is the design matrix that relates the slip and the surface displacements, and Σ −1 d is the inverse of the data covariance matrix. Subsampled data points from the atmospheric corrected InSAR data that were used in the inversions, shown in kilometric Cartesian coordinates. a, b ALOS-2 ascending and descending subsampled data points, respectively. c, d Sentinel-1A ascending and descending subsampled data points, respectively Page 7 of 15 Ghayournajarkar and Fukushima Earth, Planets and Space (2020) 72:64 Note that the matrices G(m) and L(m) are functions of the geometrical model parameters. The elements of the design matrix G(m) were computed using the solution of angular dislocations in an elastic half-space (Thomas 1993). We followed the definition of the smoothing matrix adopted by Fukushima et al. (2013).
The data covariance matrix Σ d accounts for the uncertainties in the observed displacements and their correlation. We assumed that the data noise follows an exponential correlation function (Fukushima et al. 2013). The (i,j) element of the covariance matrix for each of InSAR data was assumed to be: where σ 2 is the variance, a is the correlation length and r ij is the distance between the ith and jth data points. The variance for each InSAR dataset was set as calculated value in the non-deformed areas of the interferograms (Table 1, Additional file 1: Figures S1 and S2). As for the correlation length, using the values calculated from the noise (4-7 km, Additional file 1: Table S1, Figure S2) led to appearance of spatially correlated residual south of the bull's-eye deformation patterns (Additional file 1: Figure  S3), probably due to less penalty given to relatively shortwavelength but spatially correlated misfit. Since the displacements in this area were important for constraining the fault dipping direction, we set the correlation length to 1 km (Table 1) after trials and errors so that the spatially correlated residual was suppressed.

Inversion algorithm
The nonlinear and linear parameters were solved using the Particle Swarm Optimization (PSO) and the damped least-squares method with non-negativity constraints, respectively, and simultaneously, following the method of Fukushima et al. (2013). With each set of nonlinear parameters selected by the PSO, non-negative least-squares (NNLS) inversion was performed to determine the slip distribution and to obtain the likelihood (Eq. (1)) corresponding to the set of model parameters.
The likelihood values obtained in each iteration step is then used by the PSO to choose the next sets of nonlinear parameters. This procedure was repeated until convergence.
The uncertainties of the nonlinear model parameters were evaluated by the nearest-neighborhood interpolation of the likelihood function (Sambridge 1999). The uncertainties of the linear model parameters were evaluated by the model covariance matrix in the least-squares inversion scheme (Menke 1986). The inversion composed of the processing chain explained above was conducted separately assuming the northeast-and south-dipping faults.

Results
Figure 5 (top) shows how the PSO search for the case of the northeast-dipping fault model, converged toward a maximum in the likelihood function. A total of 462 iterations were made with Np (number of particles) = 12, meaning that the search involved the evaluation of 5544 sets of nonlinear parameter values. After carrying out the nearest-neighbor approximation and redistributing 10,000 points in such a way that the density of points is proportional to the posterior probability density function (Sambridge 1999), 95% confidence intervals of the nonlinear parameters were computed by discarding the top and bottom 2.5% in the one-dimensional marginal probability density distributions (Fig. 5 bottom, Table 3). The distributions are close to Gaussian, indicating that the nonlinearity of the parameters is weak. The distributions and resulting confidence intervals indicate that the parameters are well constrained. For the case of the northeast-dipping fault, the slip was distributed on an elongated area of 12 × 4 km 2 (Fig. 6a), with a maximum of about 1.35 m at a depth of about 5 km. The main slip was localized between depths of 4 to 7 km (Fig. 6c). For the case of the south-dipping fault, the dominant slipped area was 14 × 6 km 2 (Fig. 6b), with a maximum of about 0.80 m at a depth of about 7 km. The main slip was localized between depths of 3 to 10 km for this case (Fig. 6d). The estimated uncertainties of the slip distributions calculated after the atmospheric corrections indicate that the slip might extend deeper, but only the obtained slip shallower than 7 km is significant (Fig. 7). As Fig. 7 shows, the levels of standard deviation of slip are larger for the northeast-dipping fault model than the south-dipping one. This can be attributed to the lack of data east to the fault location on the ALOS-2 ascending interferogram (Fig. 2a). The missing area is important for constraining the slip on the deeper parts of the northeastdipping fault plane, leading to a larger standard deviation of slip for this fault plane.
The values of the released geodetic moment of the northeast-and south-dipping fault models were 1.36 × 10 18 Nm (Mw 6.02) and 1.31 × 10 18 Nm (Mw 6.01), respectively, assuming a shear modulus of 30 GPa. These values are consistent with the Global Centroid Moment Tensor (GCMT) solution (Table 3).
For the slip models of the two nodal planes, we found that they both show a good fit with the observed coseismic displacements (Figs. 8 and 9). Both the northeastand south-dipping fault models gave the root mean square (RMS) of the residuals of approximately 0.75 cm for the ascending and descending ALOS-2 datasets and approximately 0.35 cm for the ascending and descending Sentinel-1A datasets, resulting in the average RMS residuals of 0.53 cm and 0.55 cm for the northeast-and south-dipping fault models, respectively (Table 3). The fitting was hence slightly better for the northeast-dipping fault model. Since the strike of the south-dipping fault deviates from the strike of the nearby mapped active fault traces (Fig. 1), we also conducted another inversion assuming a southwest-dipping fault with the strike fixed to those of the mapped active faults in the region to examine whether this assumed fault strike would give an acceptable fit with the observed coseismic displacements. The results of this inversion indicated systematic residual patterns (Fig. 10) and a larger RMS residual (Table 3), which led us to conclude that this southwest-dipping fault is unlikely.

Fault dipping direction
The inversion results were consistent with the reported focal mechanisms of USGS, GCMT and IRSC (Table 3), either dextral-reverse rupture on a northeast-dipping fault or a sinistral-reverse rupture on a south-dipping fault. In our results, the northeast-dipping fault model had a smaller residual (0.53 cm) and was more consistent with the strikes of the mapped active faults in the region (Fig. 6a). The southwest-dipping fault model fixing the strike to the mapped active faults in the region had a larger RMS residual (0.68 cm) and gave systematic residual patterns (Fig. 10).
We also recognized that the aftershocks were more consistent with the strike and dip angles of the northeast-dipping fault model (Fig. 6e), although the aftershocks were systematically deeper than the fault plane determined by the inversion. Further, the strike of the northeast-dipping fault model was parallel to the mapped active faults in the region.
From the above reasons, the 2017 Sefid Sang earthquake likely ruptured a northeast-dipping fault. It can also be said that determining the dipping direction of a blind reverse (or normal) fault using InSAR data alone, even with multiple observation directions, is challenging at least with today's InSAR measurement accuracy. As we notice in the modeled interferograms in Figs. 8 and 9, the fringe patterns to the north and to the south of the bull'seye patterns are quite different for the northeast-and south-dipping fault models. In a future, when more effective mitigation of the atmospheric noise is realized, we will perhaps be able to determine the dipping direction of the blind faults using (atmospheric-noise free) InSAR data alone.
The relative locations of aftershocks have been shown to be capable of estimating the fault dipping plane of the mainshock, if the seismic network was dense enough (e.g., Waldhauser and Ellsworth 2000;Yang et al. 2009;Yukutake and Iio 2017). The aftershock distribution, even after relocation, shows only slight preference on the northeast-dipping fault (Fig. 6). Denser spatial coverage of seismic stations and usage of appropriate velocity structure models would be needed to sharpen the aftershock distribution, from which the fault dipping direction of the mainshock could be estimated with higher confidence.

Atmospheric delay correction
Atmospheric effects represent one of the major error sources of repeat-pass InSAR and could mask actual displacements due to tectonic or volcanic deformation. As we discussed in the previous section, atmospheric delay corrections can be a key factor for determining the dipping direction of the fault since the spatial wavelength and pattern of the important signals for the distinction are similar to the atmospheric delay. In this study, we applied the tropospheric delay estimated by the GACOS  Fukushima Earth, Planets and Space (2020)  service, which utilizes the Iterative Tropospheric Decomposition (ITD) model. Our atmospheric delay correction appeared to be effective partially, on descending interferograms, removing some atmospheric signals near the coseismic deformation consistently for both ALOS-2 and Sentinel-1A datasets (Fig. 3). All of the northeast-, southand southwest-dipping fault models without atmospheric corrections gave RMS residuals of 0.66, 0.70 and 0.82 cm, respectively, which were reduced to 0.53, 0.55 and 0.68 cm, respectively, with the atmospheric corrections. The correction effect may have been limited by the presence of short-wavelength turbulence in GACOS products ), which could have been the case for the ascending interferograms that appeared to show no improvements.

Seismotectonics implications
In this study, the aftershocks were systematically deeper than the fault plane determined by the inversion conducted on the InSAR data, although the aftershock distribution suggested northeastward fault dipping direction. This implies that the assumed crustal velocity in the studied area may have been systematically larger than the true value and thus may have overestimated the depth of the hypocenters.
Such deviation of the velocity structure from the reality may be due to the presence of thick deposits in this region. The thickness of the deposits in the Sefid Sang section of the Kashafrud formation is 1177 m (Poursoltani 2017). The presence of this thick deposit layer may be the reason why the rupture did not reach the ground surface. If so, we can hypothesize that blind faulting is ubiquitous in this region and the fault traces cannot be clearly identified from the surface morphology. This further implies that we cannot heavily rely on the mapped active faults to evaluate the seismic hazard in this region, although the strikes of the mapped active faults are useful for comparing with the candidate models.

Conclusions
In this study, we investigated the potential of InSAR data for determining the dipping direction of a blind reverse fault caused by the 2017 Sefid Sang earthquake. We used four available interferograms measuring displacements from multiple directions, two from Sentinel-1A data and two from ALOS-2 data. We applied the GACOS atmospheric noise correction, which partially removed the atmospheric noise in the descending interferograms. We explored the possibility of different dipping directions by conducting Bayesian inversions assuming both northeast-and south-dipping fault planes.
Our results showed the northeast-dipping fault model had a smaller RMS residual (0.53 cm) and was more consistent with the strikes of the mapped active faults in the region. We also showed that the aftershocks were more consistent with the strike and dip angles of the northeast-dipping fault model. From the above reasons, we conclude that the 2017 Sefid Sang earthquake likely ruptured a northeast-dipping fault. The preferred fault model has a strike angle of 314.8°, dip angle of 47.4° and rake angle of 130.3°, and a slip distribution of maximum 1.35 m at depth of 5 km equivalent to Mw 6.0.
This study illuminated that determining the dipping direction of a blind reverse (or normal) fault using InSAR data alone, even when displacements are measured from multiple directions, is challenging at least with today's InSAR measurement accuracy. However, correcting for noise and comparing with other kinds of data, such as aftershock distribution and active fault traces, can help infer the fault dipping direction.