Constraints on the 2013 Saravan intraslab earthquake using permanent GNSS, InSAR and seismic data

SUMMARY

On 16 April 2013, an M_w = 7.7 earthquake struck the border of Iran and Pakistan in the central part of the Makran subduction zone with a reported depth of 80 km by USGS. This rare event in this poorly instrumented region helps to shed light on the kinematics of the subducting slab. We investigate source parameters of the Saravan intraslab normal earthquake using RADARSAT-2 SAR images in three ascending tracks, nine permanent GNSS sites and teleseismic data. The maximum coseismic displacement occurred at the SRVN GNSS station with 54.1 mm southeast horizontal and 42.7 mm upward vertical displacements. The coseismic ascending InSAR displacement maps illustrate a continuous and smooth NE-trending elliptical shape deformation pattern with a maximum of ∼29 cm of displacement away from the satellite. We use 25 broad-band teleseismic P-waveforms to estimate the focal mechanism of the main shock. A joint uniform inversion of InSAR, GNSS and teleseismic data reveals a NW-dipping SW-striking fault and a primarily normal-faulting earthquake with a minor right-lateral strike-slip component. The static slip distribution of the InSAR coseismic maps localizes variable slip at depths between 50 and 81 km with a maximum amplitude >3 m at 60–75.5 km depth, rupturing the oceanic crust of the subducted slab. The kinematic slip distribution exhibits a well-constrained slip pattern with a nucleation depth of 65 km. The source time function indicates that the earthquake reaches its maximum moment tensor release at ∼8 and ∼16 s. The NE-trend of the Saravan earthquake slip pattern, the orientation of the volcanic arc, and the distribution of the intraslab intermediate-depth normal earthquakes provide new insights into slab geometry in the central Makran subduction zone. We suggest that the slab bending at the hinge of subducting Arabian Plate is oblique along a NE–SW direction parallel to the volcanic arc rather than the shoreline or deformation front, and it is likely to be the reason for an oblique volcanic arc in the Makran subduction zone. These new constraints on the Makran slab geometry will help further studies in establishing realistic coupling maps for seismic hazard assessment.

1 INTRODUCTION

The 2013 April 16 M_w 7.7 Saravan earthquake occurred 82 km NW of Saravan city in the vicinity of the border of Iran and Pakistan (Fig. 1) with a reported depth of 80 km (USGS catalogue, https://earthquake.usgs.gov). The earthquake was followed by only eight aftershocks with magnitudes between 4 < M_w < 5.7 until 2023 based on the ISC-EHB Bulletin (Engdahl et al. 1998, 2020; Weston et al. 2018), and the most recent one occurred on 25 December 2013 (depicted by circles in Fig. 1b). The Saravan earthquake is the largest event in Iran in the last 50 yr. It happened in a sparsely populated area and caused one death and 12 injuries in Iran, and 40 deaths and 180 injuries in Pakistan (Zare et al. 2013). The reported USGS depth and focal mechanism suggest a normal-faulting earthquake ruptured the subducted Arabian oceanic lithosphere in the Makran subduction zone. Barnhart et al. (2014) used InSAR data and teleseismic waveforms independently to derive causative fault parameters and slip distribution. They determine a N-dipping fault plane, with variable slip between 30 and 90 km depth and maximum slip patch at 50–80 km depth (Barnhart et al. 2014). They proposed that half of the oceanic lithosphere, consisting of the mantle and the crust, was ruptured during this earthquake with the occurrence of the slip distribution in the oceanic mantle. Using the intersection of the 2013 Saravan event slip pattern with the thermal model of Smith et al. (2013), Barnhart et al. (2014) suggested that the dehydration of hydrous minerals and slab bending forces are the likely deriving mechanisms of the Saravan event.

It's noteworthy that the seismic hazard associated with events such as the Saravan earthquake differs from the ongoing discussions about the megathrust's ability to host large-magnitude events, potentially tsunami-triggering earthquakes. Although the Saravan earthquake does have implications for seismic hazard, particularly given its depth and the potential for its effects to be felt over a wide region, the main focus in examining this event lies in its seismotectonic significance. Furthermore, the east–west-oriented instrumental thrust earthquakes, parallel to the subduction front, exhibit a distinct pattern from the northeast-trending normal-faulting earthquakes and volcanic arc (Fig. 1a). This complex deep structure underscores the need for comprehensive seismic studies of the area. The Saravan earthquake occurred in the border region of Iran and Pakistan, where the crustal structure associated with the subduction zone remains largely unconstrained due to the lack of local instrumentation and detailed studies. It is practically impossible to install seismic stations to study the deep structure and slab geometry in detail due to the logistical inaccessibility in the vicinity of the Iran–Pakistan border. Recently, the seismic sites located in the western Makran subduction zone have been used to investigate the crust and uppermost mantle structure (e.g. Motaghi et al. 2020; Priestley et al. 2022). Therefore, the intraslab Saravan event that occurred on the border with a normal-faulting mechanism offers more information about the slab interface geometry in this region, and it converts the 2-D geometrical information of the subduction zone from seismic data in the western part (Motaghi et al. 2020; Priestley et al. 2022) into a 3-D view.

In this study, we initially use the permanent GNSS network run by the National Cartographic Center (NCC) of Iran to estimate the coseismic surface displacements related to the 2013 Saravan earthquake. Secondly, we take advantage of the newly provided orbital information by the Satellite operator (MDA Corporation) and updated GAMMA software to generate coseismic interferograms using RADARSAT-2 images. Thirdly, we use Teleseismic P-waveforms to determine the focal mechanism of the earthquake. Fourthly, we combine all data sets to derive source parameters of the Saravan earthquake causative fault with a Bayesian bootstrap optimization non-linear inversion approach. Fifth, we estimate the static slip distribution along the retrieved causative fault by inverting all three coseismic InSAR displacement maps and kinematic slip distribution by inverting InSAR and teleseismic P-waveforms. Particularly, these steps are different from Barnhart et al. (2014) in five ways. First, we present the GNSS time-series of nine sites located in the Makran subduction zone and their coseismic offset due to the Saravan event. Secondly, the orbital error contribution is reduced in coseismic interferograms using newly provided orbital information and updated features of GAMMA software. Thirdly, we combine GNSS and InSAR data along with teleseismic waveforms to invert them jointly, considering NW-dipping and SE-dipping fault planes based on the probabilistic earthquake source inversion approach. Fourthly, we use the joint inversion technique to determine both static and kinematic slip distribution. This involves analysing three InSAR coseismic interferograms for static slip and utilizing InSAR and teleseismic waveforms for kinematic slip. Finally, this detailed analysis of the coseismic slip characteristics of the Saravan earthquake combined with updated information on the deep structure and geometry of the western Makran subduction zone based on new seismic images (Motaghi et al. 2020; Priestley et al. 2022) helps to improve our knowledge of the slab geometry in the central part of the subduction zone, where there is logistical inaccessibility for installing the local seismic station.

2 TECTONIC SETTING

The Makran subduction zone results from the convergence of the Arabian oceanic plate and the Eurasia in southeast Iran and south of Pakistan. Global navigation satellite system (GNSS) measurements indicate that the subduction zone accommodates ∼20 mm yr⁻¹ of shortening between Muscat in Oman and Chabahar, located at the Iranian coastline (Fig. 1a), while their sparse distributions make it hard to study the interseismic coupling of the Makran subduction zone interface (Vernant et al. 2004; Frohling & Szeliga 2016; Penney et al. 2017; Khorrami et al. 2019). Using recently installed GNSS sites in the western Makran subduction zone, Abbasi et al. (2023) determined that the interseismic coupling between the subducting Arabian Plate and the overriding continental plate in the western segment of the Makran subduction zone is significantly less, by a factor of four, in comparison to its eastern segment. Lin et al. (2015) studied the interseismic coupling of the eastern Makran megathrust by applying time-series to Envisat images. They discovered high coupling in the central section of the Eastern Makran, indicating the potential for significant earthquakes (magnitude 7+), and lower coupling aligns spatially with the subduction of the Sonne Fault Zone (Kukowski et al. 2000). The northern extent of the EW-trending subduction zone is limited by two Quaternary covers of the Jazmurian depression in Iran and the equivalent Mashkel depression in Pakistan (Fig. 1a). The NE-trending volcanic arc consists of the Bazman and Taftan volcanoes in Iran as well as Soltan in Pakistan (green triangles in Fig. 1a). The region north of the Makran subduction zone is divided into two parts by the NS-trending Sistan Suture Zone that implies different kinematic boundary conditions in western and eastern Makran (Byrne et al. 1992, SSZ in Fig. 1a), although these structures terminate at the northern margin of the accretionary prism (Jackson & McKenzie 1984; Walker et al., 2013; Penney et al. 2017). Pajang et al. (2021) investigated the seismic hazard of the western Makran by analysing the deformation of the accretionary wedge using three N–S seismic profiles and the frictional properties of the megathrust, suggesting the existence of seismic asperities in the western Makran, which could lead to a large earthquake. The instrumental thrust earthquakes (blue focal mechanisms in Fig. 1a) extend from west to east along the coast parallel to the subduction front and are primarily concentrated in eastern Makran. It's worth noting that these thrust earthquakes could potentially be aftershocks of the only recent megathrust earthquake in the region, which occurred on 27 November 1945. The normal earthquakes (red focal mechanisms in Fig. 1a) extend inland with a rough northeast trend, aligning parallel to the volcanic arc (green triangles in Fig. 1a). This pattern suggests the presence of a complex deep structure and slab geometry in the region, further emphasizing the need for detailed studies and investigation into the seismotectics of this area.

Geodetic and seismic data have been used to assess the mean dip of the subduction interface and the deep structure in the Makran, but there are still open questions. The early studies based on seismic data proposed ∼3° as the dip of the oceanic plate at ∼62.5°E (White & Louden 1982; Kopp et al. 2000) and a thick accretionary prism in the Gulf of Oman and coastal Makran. Penney et al. (2017) suggest the average dips of ∼11° in the western Makran (beneath Bazman) and ∼9° in the central region (beneath Taftan and Soltan volcanos) using subduction-related earthquakes and a shear-wave velocity profile beneath Chabahar. It's worth noting that these dip angles, although steeper than those typical for the oceanic Moho or the base of deformable sediment near the trench, are not surprising. This is because they cover a substantial inland distance where the subducting slab is anticipated to steepen. In other words, as we move farther inland from the trench, it is expected that the dip angle of the subducting slab becomes progressively steeper. A regional Pn-tomography study for the deep structure of the Middle East and Iranian plateau shows a contrasting Pn velocity anomaly under the eastern and western parts, leading to the claim of a different characteristic of the eastern and western Makran subduction zone (Al-Lazki et al. 2014), separated by the SSZ. Using teleseismic data, Motaghi et al. (2020) investigate the deep structure and the geometry of the velocity interfaces beneath the Makran subduction zone through six seismic stations located in the coastal Makran (at the latitude of 25.5°–26.5°N). They observed an interface at a depth of 9–15 km dipping northward at a dip of ∼2° with a sharp change in dip south of the Qasr-e Qand thrust fault (QQF in Fig. 1a). With analysis of teleseismic body waves reflected and converted in the oceanic crust, they estimated a 4–7 km thickness for the oceanic crust. Priestley et al. (2022) investigate the crust and uppermost mantle structure of the western Makran based on the analysis of teleseismic body-wave and local and regional surface wave recordings from a network of 26 temporary and permanent seismometers operated from 2016 to 2020 (black triangles in Fig. 1a). They found a 22–28 km sedimentary cover near the coast. Their estimated dip of oceanic Moho in the coastal Makran is 4° ± 2°. At the southern edge of the Jazmurian depression, the dip of subduction steepens, and the slab dives down with a dip of ∼18° beneath the continental crust of the southern Iranian Plateau, which has an almost flat Moho at a depth of 4045 km (Priestley et al. 2022).

3 DATA AND PROCESSING

Geodetic observations (GNSS and InSAR) help to estimate the static parameters of seismic sources, while teleseismic waveforms with good temporal resolution offer valuable information about the depth of the rupture and its time evolution. Therefore, geodetic and seismic data exhibit complementary information for studying seismogenic fault parameters and source rupture processes. Here, first, we present the GNSS data that captured the Saravan earthquake; secondly, we establish the coseismic displacement map by InSAR; thirdly, we estimate the focal mechanism solution using the broad-band P-wave data and finally, we combine all the information to investigate the source parameters and the rupture process.

3.1. GNSS coseismic displacement

The Iranian Permanent GNSS Network (IPGN) was established by the National Cartography Center (NCC) after the devastating 2003 Bam earthquake in southeastern Iran. All GNSS stations were set up on geodetically designed concrete pillars deeply rooted in stable ground (Khorrami et al. 2019). The observed GNSS data of nine permanent stations of IPGN (NEHB, ZABL, FHRJ, GLMT, SRVN, ANGN, NKSR, JASK and CHBR) situated in the study area have been analysed by the GAMIT/GLOBK software, version 10.6 (Herring et al. 2015). To stabilize the reference frame, the IPGN network is combined with 37 GNSS stations of the International GNSS Service (IGS) network in the daily solutions. We calculate daily positioning solutions along with the tropospheric zenith delays, adjusting IGS final orbits and using coherent Earth Orientation Parameters (EOP). We use the tropospheric mapping function VMF1 (Boehm et al. 2006) to calculate the tropospheric zenith delay every 2 hr. Then, we use the FES2004 ocean tide model for ocean loading (Lyard et al. 2006) and the absolute antenna phase centre model (IGS_05) for variation of antenna phase centres in the daily solutions. To obtain a robust reference frame for the coordinate time-series of the Iranian stations, we combine the resultant daily solutions with global network solutions using the Kalman filter implemented in the GLOBK package (Herring et al. 2015). An individual noise model for each permanent station is calculated based on first-order Gauss Markov (FOGM) process noise (Gelb 1974; Reilinger et al. 2006; Mousavi et al. 2013; Walpersdorf et al. 2018; Khorrami et al. 2019). We use this FOGM model to predict the uncertainty of each site by considering the length and the noise of the time-series. Finally, we estimate the coordinate time-series of each station with respect to Eurasia using the Eurasia-fixed coordinates and velocities of ITRF2014 (Altamimi et al. 2016) for the IGS reference stations. Fig. 2 illustrates the NS component of the time-series spanning from 2006 to 2021 for nine stations. Each data has been detrended using its interseismic displacement rates predating the Saravan earthquake on 16 April 2013, under the assumption that the interseismic rate follows a linear trend over time.

To calculate the coseismic offset, we subtract the mean position of each site five days before and after the earthquake. The standard deviations of coordinates of five days before and after the earthquake are used to estimate the uncertainty for the co-seismic offset following the error propagation rule. The time-series presents coseismic offsets during the M_w 7.7 Saravan event at all stations with maximum values at SRVN, ZABL, NEHB, GLMT and NKSR. Fig. 3 and Table 1 present the coseismic offsets and their uncertainties for the nine permanent GNSS stations at a distance of up to 450 km from the epicentre. Stations at larger distances (such as FHRJ, CHBR and JASK) were slightly displaced. The ANGN station was not affected by the Saravan event but by the 11 May 2013 M_w 6.1 Minab earthquake (Figs 2 and 3). The maximum coseismic displacement occurred at the SRVN station with 54.1 ± 1.9 mm horizontal displacement in a southeast direction (17 ± 1.3 and 51.4 ± 1.3 mm towards east and south, respectively) and 42.7 ± 2.6 mm upward vertical displacement. The second-largest horizontal displacement of about 12 ± 1.6 mm northwestward was observed at the ZABL site, located 300 km north of the epicentre (Table 1).

3.2. InSAR coseismic displacement map

Space-borne Interferometric Synthetic Aperture Radar (InSAR) has been extensively used in recent decades to monitor and measure ground displacement due to geophysical phenomena, including earthquakes, subsidence, salt diapir, and volcano eruptions (e.g. Colón et al. 2016; Amiri et al. 2020; Davis et al. 2021; Mohammadnia et al. 2021; Park & Hong 2021; Vajedian et al. 2023). InSAR-based surface displacement measurements require careful interpretation due to potential error sources, including topography, orbital inaccuracies, phase unwrapping and atmospheric disturbances (Massonnet & Feigl 1998). Accurate correction and removal of these unwanted signals during data processing are essential.

To measure the 2013 Saravan coseismic deformation, we use SAR images acquired by the RADARSAT-2 satellite in Wide Multi-Look Fine Resolution Beam mode in two ascending tracks, MF21W and MF22W (Table S1). Three coseismic interferograms are generated by using a set of images from both the MF21W track and the MF22W track. For the first interferogram (Int. 1), the pair covers 384 d from 29 March 2012 to 17 April 2013, for the second interferogram (Int. 2), it covers 408 d from 29 March 2012 to 11 May 2013, and for the third interferogram (Int. 3) it covers 408 d from 22 March 2012 to 04 May 2013 (Table S1). The preprocessing of SAR images and the interferogram generation procedure are performed with the GAMMA software (Wegmüller & Werner 1997). We remove topographic phase contributions from the interferograms using the 30 m ASTER Digital Elevation Model (DEM, Welch et al. 1998). To use InSAR displacement maps in an earthquake source parameter inversion, it is necessary to reduce the residual orbital ramp and tropospheric phase contributions. We remove the orbital ramp contribution using newly provided orbital information by the Satellite operator (MDA Corporation) and updated features of GAMMA software to generate coseismic interferograms, which is one of the improvements in the processing step compared to previously generated interferograms (Barnhart et al. 2014). We remove tropospheric phase delay for each track using Generic Atmospheric Correction Online Service for InSAR (GACOS, Yu et al. 2017, 2018a, b) implemented in the Kite package (Isken et al. 2017). We apply adaptive filtering (Goldstein & Werner 1998) to interferograms and unwrap the interferograms using a minimum-cost flow method (Costantini 1998). The final unwrapped and geocoded interferograms are represented in Figs 4(a)–(c). Finally, we establish a reference point located significantly away from the deformation zone to ensure minimal displacements in the far field. We present the reference point with pink hexagons for Int.1 and Int. 2 in Figs 4(a) and (b). For Int. 3, the reference point is located at latitude 29.6°N.

LOS displacement maps in three ascending interferograms (panels a, b and c in Fig. 4) show the coseismic displacement with a maximum value of ∼29 cm moving away from the satellite. The continuous and smooth NE-trending elliptical shape of the deformation pattern indicates a buried source with a normal mechanism roughly parallel to the volcanic arc rather than to the EW-oriented deformation front of the Makran subduction. The triangle in Fig. 4(a) shows the location of the SRVN GNSS station. Furthermore, the coseismic displacement maps (Figs 4a, b and c) reveal no surface rupture due to this earthquake.

Figs 4(a) and (b) also exhibit the epicentre of the 2013 Saravan earthquake from ISC-EHB, CMT, USGS, and GFZ catalogues, as well as the distribution of its aftershocks based on the ISC-EHB catalogue in the four days after the earthquake up to 20 April 2013. The greatest aftershock occurred about 4 hr and 45 min after the main shock with M_w 5.6, at a depth of 60 km based on the CMT catalogue. The pattern of the aftershock distribution shows a NE-trend. The coseismic interferograms cover seven aftershocks with moderate magnitude occurring in the time interval of the interferograms. Given that Int. 1 and 2 were generated from the same primary SAR images with 24 d difference in time but covering similar regions, their difference in line-of-sight (LOS) displacement is anticipated to reveal the early post-seismic displacement linked to the aftershocks. However, subtracting Int. 2 from Int. 1 displays near-zero LOS displacement, as shown in Fig. 4(d). This finding suggests that early post-seismic deformation can be neglected in comparison with the coseismic displacement of the main shock in earthquake source parameter inversion.

3.3. Teleseismic seismograms

We also gather broad-band teleseismic recordings obtained from Incorporated Research Institutions for Seismology (IRIS, Trabant et al. 2012) and GEOFON seismic network selected within epicentral distances of 30°–90° equivalent to 3000–10 000 km (red triangles in Fig. 5a). By excluding noisy records, overall, 25 stations are collected to have an evenly azimuthal coverage with high signal-to-noise ratio. We apply a Butterworth bandpass filter between 0.01 and 0.05 Hz to the vertical components of the seismic waveforms. The traces were windowed for 25 s, starting 5 s before the P-wave arrival time, to include in the data fitting during the inversion. In the following, we use these seismic data sets to estimate the full, deviatoric, and double-couple moment tensor of the Saravan main shock. Then, they are jointly inverted with geodetic data to derive the earthquake rupture process.

4 ESTIMATION OF EARTHQUAKE SOURCE PARAMETERS

In recent years, joint inversion of seismic and geodetic data has been extensively used to explore fault source parameters (e.g. Delouis et al. 2002; Wang et al. 2018; Zheng et al. 2020). The seismic waveforms include the fault's rupture history, in particular, origin time and velocity (Ide 2007), and on the other hand, geodetic data provide information about fault geometry and its static slip distribution (Zhang et al. 2012; Steinberg et al. 2020). In the following, we use teleseismic waveforms to derive a focal mechanism solution (point source modelling) for the 2013 Saravan earthquake. Then, we jointly utilize teleseismic waveforms, GNSS and InSAR data to derive the fault parameters of the 2013 Saravan earthquake in the first step, and then, in the second step, we discretize the fault into multiple patches to constrain the static and kinematic slip distribution along the fault plane through an inversion of InSAR data and joint inversion of teleseismic waveforms and InSAR data, respectively.

We use Grond (Heimann et al. 2018), an open-source Python software package for probabilistic earthquake source inversion based on the Pyrocko package (Heimann et al. 2017), for moment tensor inversion and the joint inversion of seismic and geodetic data. The Grond implements a Bayesian bootstrapping optimization (BABO, Rubin 1981; Heimann 2011) to explore the full model space and map model parameters with relevant uncertainties and trade-offs by offering a flexible design of misfit functions. Grond uses the pre-computed Green's functions to accelerate forward modelling, which is handled by a command line tool, Fomosto, in the Pyrocko package (Heimann et al. 2017). To derive the elastic parameters of the layered medium, we calculate dynamic Green's functions for far-field data (i.e. teleseismic waveforms) using the QSSP code (Wang et al. 2006) and static Green's functions for near-filed data (i.e. geodetic data) using PSGRN/PSCMP code (Wang et al. 2017) implemented in the Pyrocko toolbox (Heimann et al. 2019).

The misfit between observed and modelled data is represented by the objective function, which is the global minimum searched during the optimization. Furthermore, the objective function has rules for handling different observed data sets, which Lp-norm is applied, and how data errors are considered in optimizations. BABO starts exploring model space, taking samples from a uniform distribution within given parameter bounds. All bootstrap realizations are evaluated in parallel. Based on the performance of the tested models, the proposal distributions for new models are gradually and independently updated for every bootstrap realization until they become smaller or equal to the continuously updated estimates of the parameter distributions (Kühn et al. 2020). This Bayesian bootstrap-based method has been successfully applied in several recent studies in different regions (e.g. Dahm et al. 2018; Jamalreyhani et al. 2019, 2020, 2021, 2023; Büyükakpınar et al. 2021; Cesca et al. 2021; Petersen et al. 2021) and described in previous publications (e.g. Dost et al. 2020; Kühn et al. 2020).

4.1. Focal mechanism solution

We test for full, deviatoric and double couple (DC) moment tensor inversion to estimate focal mechanism solution. The parameters of the moment tensor inversion include centroid location, depth, time, magnitude, six relative moment tensor component ranges and a duration parameter controlling a half-sinusoid shape source time function. We use the ISC-EHB solution as a prior for the centroid north and east variation values and also allow depth, time, duration and magnitude (M_w) to vary between 40 and 100 km, −10 to 10 s relative to hypocentre origin time, 10–50 s and 5.7–8, respectively, according to the CMT solution. Overall, we test 55 000 models for the moment tensor inversion, with 10 000 models drawn randomly based on a uniform distribution, while 45 000 models were drawn using a direct search for selected high-score models. The velocity structure in the Saravan region lacks robust constraints. Additionally, when using centroid moment tensor inversion with teleseismic waveforms, choosing a global green function has no impact on the results. Therefore, we define Earth layering based on the AK135 1-D velocity model (Kennett & Engdahl 1991) with source depths between 1 and 150 km and spacing of 5 km. The synthetic seismograms are calculated with a grid spacing of 5 km in the epicentral distance of 1–10 000 km (Fig. 5b). For determining the weight of each station in teleseismic data analysis, we consider the influence of both data noise and the distance from the epicentre on the signal amplitude. The best solution of this framework is the minimum misfit function between synthetic and observed seismic waveforms.

For full, deviatoric and DC moment tensor inversion, we plot seismic station distribution and decomposition of the best-fitting solution to Isotropic (ISO), compensated linear vector dipole (CLVD) and double-couple (DC) components. Fig. 5(a) exhibits the observed and modelled traces of DC moment tensor inversion for the 2013 Saravan main shock. Fig. 5(a) shows the target information from top to bottom in each trace, station name, distance to the source on the left, as well as the azimuth of the station with respect to the source and target misfit on the right-hand side. The waveforms presented in Fig. 5(a) include the restituted to displacement and filtered observed trace with tapering and processing (black) along with the filtered, tapered, shifted and processed synthetic target trace (red). In each waveform, the coloured boxes on the upper right-hand side illustrate the relative weight of the target within the entire data set of the optimization (top box, orange) and the relative misfit contribution to the global misfit of the optimization (bottom box, red). Fig. S1 presents the station distribution, solution of DC moment tensor inversion and decomposition of it in ISO, CLVD and DC parts, the Hudson plot with the ensemble of bootstrap solutions, the location plot and the convergence sequence of each parameter for the DC solution, indicating that all parameters converged very well during the moment tensor inversion. Fig. S2 illustrates the distribution of parameters and the trade-off between them. Moreover, we illustrate observed and synthetic traces of each station for full (Fig. S3) and deviatoric (Fig. S4) moment tensor inversion. Finally, we depict the station distribution, solution of DC moment tensor inversion and decomposition of it in ISO, CLVD and DC parts, the Hudson plot with the ensemble of bootstrap solutions, location plot, and convergence sequence of each parameter for the full (Fig. S5) and deviatoric (Fig. S6) solution. The Hudson plots of full, deviatoric, and DC moment tensors indicate that the double-couple component is significant, and we consider the result of the DC moment tensor in the following.

The results of DC moment tensor inversion reveal that the earthquake is located at 61.94°E and 28.01°N with a moment magnitude of M_w 7.7 (red beach ball in Fig. 4a). The depth of the main shock is 79 ± 5 km, deeper than CMT, ISC and GFZ solution values and close to the USGS solution value (Table 2). We derive a normal focal mechanism solution for the 2013 Saravan earthquake with the strike, dip, and rake equal to 237° ± 5°, 63° ± 2° and −101° ± 7° for the first nodal plane and 83° ± 16°, 28° ± 4°, −67° ± 10°, for the second nodal plane, respectively, based on only teleseismic data, which are consistent with mechanisms derived from other common institutes (Table 2). In the next section, we estimate the earthquake source parameters by combining the near-source high spatial resolution InSAR data with far-field high temporal resolution teleseismic observations.

4.2. Joint uniform slip modelling

To obtain the parameters of the causative fault for the 2013 Saravan earthquake, we use jointly three interferograms presented in Fig. 4 and teleseismic waveforms prepared at the moment tensor calculation using the Pyrocko toolbox (Heimann et al. 2017). Regarding the selection of GNSS stations for inversion, we precisely evaluated the reliability of each station's data. Out of the nine available stations near the epicentre, only four (SRVN, ZABL, NEHB and NKSR) were included in the inversion due to their substantial displacements exceeding uncertainties. We exclude GLMT, CHBR, ANGN, FHRJ and JASK stations from the inversion due to their low signal-to-noise ratio. Their coseismic displacements are smaller than the associated uncertainties, making the data unreliable for this analysis. However, tests were conducted with all stations included in the inversion to assess their reliability, confirming minimal impact on the modelled fault parameters and justifying the exclusion of stations with poor signal-to-noise ratio. This behaviour is attributable to the Grond software package, which incorporates the data set's variance-covariance matrix during inversion.

The source parameters include fault geometry (length, width and depth), fault orientation (strike and dip), location of the fault plane (easting and northing of the fault centre) and the rake and amount of uniform slip. Furthermore, we add dynamic parameters, including rupture nucleation point, time and velocity. To decrease the number of observed InSAR data points to reasonable quantities and increase the inversion speed, we use the quadtree subsampling method (Jónsson et al. 2002). Consequently, we downsample three InSAR coseismic displacement maps (Int. 1, Int. 2 and Int. 3) to 493, 438 and 415 observation points, respectively. We use spatial sampling methods and construct semi-variograms and covariograms to assess noise in regions of the displacement map unaffected by coseismic deformation. These metrics are weights in the Grond inversion process, following the methodology outlined in Jonsson & Sudhaus (2009). We consider a rectangular dislocation fault plane in a homogenous and layered medium (Okada 1985) to model the surface displacement.

Due to the high uncertainty of global Earth models in the Makran subduction zone, we use a local velocity structure retrieved by receiver function analysis of the GOSH seismic station (see Fig. 1b for the location, Mokhtarzadeh 2022) to derive the near-field static Green's functions (Fig. S8b). The source depth range is 1–120 km, with a grid spacing of 0.5 km. Furthermore, we consider the source–receiver horizontal distance equal to 600 km to cover all geodetic data with a spacing of 0.5 km. The pre-calculated far-field dynamic Green's functions for teleseismic waveforms are the same as the moment tensor inversion step.

The symmetrical shape of the coseismic displacement pattern (Figs 4a–c) suggests either SE-dipping or NW-dipping for the causative fault plane. Additional data are required to solve this ambiguity in the interpretation and modelling of the causative fault plane. Therefore, we test both NW- and SE-dipping fault planes using all available data in our modelling. The root mean square (RMS) of the InSAR interferograms for the NW-dipping fault models proved to be lower than those for the SE-dipping fault (Table S2). Based on the lower RMS, indicating a better fitting, and the aftershocks deepening to the northwest (Fig. 1c), the NW-dipping fault plane is suggested as the causative fault plane for the Saravan main shock. The detailed modelling procedure of the NW-dipping fault is thoroughly reported in the subsequent section dedicated to the inversion of source parameters. To introduce the reasonable search range for the earthquake source parameter, we take advantage of the LOS displacement maps and focal mechanism solutions. We constrain the search area of all parameters, including depth of 30–70 km, dip 40°–80°, strike 200°–360°, length 50–150 km, width 0–40 km and rake −130° to −80°. We also allow the nucleation points X and Y, the relative horizontal and along-dip position of the rupture nucleation point on the fault to the centre location, to vary between −1 to + 1, with 0 being in the centre, −1 being at the top fault edge and 1 at the bottom fault edge to the hypocentre location found from moment tensor inversion.

We apply a Bayesian bootstrap-based probabilistic approach to the available data to estimate the source parameter following a non-uniform inversion approach. After some tests, we observe a trade-off between slip and fault width in different tested models; thus, we use the information about the moment magnitude derived from moment tensor inversion to fix the slip to 4 m in this modelling step. The slip distribution will be estimated in the next linear modelling step. We consider 20 000 models for the random search space and use 80 000 iterative models to find the best solution for each parameter. We include uncertainties through a Bayesian bootstrap-based probabilistic joint inversion scheme. For determining the weight of each station in teleseismic data analysis, we consider the influence of both data noise and the distance from the epicentre on the signal amplitude. In the case of InSAR and GNSS data, the weight vectors (⁠$\mathrm{w}$⁠) and weight matrices (⁠$\mathrm{W}$⁠) are determined by considering the data error correlations expressed in the data error variance–covariance matrix $\mathrm{\Sigma }$⁠:

where $\sigma$ is the N × 1 vector of the standard deviation of the GNSS coseismic displacement vectors, $\mathrm{\Sigma }$ is the M × M matrix and M represents the number of InSAR data points (Isken et al. 2017; Heimann et al. 2018). In this inversion process, we initially assume equal weights for all data sets. Throughout the iterations, these weights are dynamically adjusted based on target balancing, noise characteristics within the data, and potentially user-defined preferences (Heimann 2011). Finally, all weights are normalized to ensure they contribute appropriately to the inversion. Fig. S7 shows the top twenty data contributions to the inversion process. It reveals that InSAR interferograms, particularly MF22 and MF21, along with SRVN GNSS station, have a substantial contribution compared to other data sets.

The comparison between observed and modelled InSAR displacement maps and teleseismic P-waveforms exhibits excellent consistency and low residual (Figs 6 and S8). Moreover, the comparison between observed and modelled GNSS displacement vectors demonstrates a reasonable fit between them (Fig. 6j) with a horizontal residual of 5.5 mm at the SRVN station, 2.4 mm for ZABL, 2.0 mm for NEHB and 4.3 mm for NKSR. The convergence plots and distribution of parameters show that most of them converged well except for the y component of the nucleation point and rupture velocity (Figs S9 and S10). Fig. S11 exhibits trade-offs between source parameters.

Table 2 presents the retrieved seismogenic fault parameters from InSAR, the joint inversion of InSAR and GNSS and their joint inversion with teleseismic data, besides the result of Barnhart et al. (2014). We find that the results retrieved following the three approaches are in good agreement by considering the uncertainties of each parameter (Table S2), and here, we refer to the results of the joint inversion of all GNSS, InSAR and teleseismic data. The resultant seismogenic fault parameters and rupture information (Table 2) from the joint inversion in this study are consistent with other solutions such as CMT, USGS, GFZ and Barnhart et al. (2014). We obtain a rectangular fault plane from the joint inversion with a top depth of ∼56 km, a length of ∼97 km and width of ∼17 km (the black rectangle in Fig. 6). The results reveal that the seismogenic fault with a strike of 234° dips to the northwest at an angle of 64°, and its surface projection (the black dashed line in Fig. 6) is located ∼50–60 km southeast of the epicentre. The rake is −111°, indicating that the slip on the fault plane has a significant normal component with a right-lateral strike-slip component. The rupture velocity is 4.5 km s⁻¹, and the origin time is equal to −1.8 s relative to the earthquake origin time of the CMT solution. The resulting RMS value from this step for Int. 1, Int. 2 and Int. 3 are 2 cm, 2.7 cm and 2.2 cm, respectively. The normalized global misfit value of all data sets is 0.3. Considering a rigidity modulus of 30 GPa, a length of 100 km, a width of 20 km and 4 m of slip, the moment release is 2.5E + 20 Nm, equivalent to an earthquake of M_w 7.6.

4.3. Slip distribution modelling

4.3.1. Static slip distribution

To retrieve a non-uniform slip model allowing the slip to vary across the fault plane, we use the fault plane derived in the previous step and we assume that all the parameters are fixed except the slip along the fault plane. As the SRVN site is the only GNSS site situated within the areal extent of the coseismic displacement map, we neglect all GNSS stations in the slip distribution inversion. We use three subsampled interferograms jointly to derive static slip distribution (Fig. S12a, d and g). In this step, we use the Bayesian Earthquake Analysis Tool (BEAT) (Vasyura-Bathke et al. 2020), which utilizes the Bayesian finite fault framework with Sequential Monte Carlo (SMC) sampling (Del Moral et al. 2006). BEAT uses Bayesian inference to estimate the parameters of a deformation source model from geodetic and seismic data. We extend fault length and width to 180 and 120 km to prevent unwanted slip at the fault plane borders and subdivide it into 600 patches with a 6 km × 6 km patch size. For static slip distribution, the source parameters include two slip components parallel and perpendicular to the fault rake for each fault patch. In static slip distribution, the relationship between the slip of the fault segments and the resultant displacement on the Earth's surface is linear, expressed by $d(m)=Gm$⁠. Here, $G$ indicates the sensitivity kernel which is a matrix comprising displacement Green's functions. Each function represents the response at the earth's surface to a unit slip in either the slip-parallel or slip-perpendicular direction on a subfault, as detailed by Aki & Richards (2002). Here, $G$ is constructed given the source geometry and the desired subfault discretization. We used Green's functions utilized in uniform slip modelling to compute Green's functions slip in parallel and perpendicular to the fault rake direction. Using Laplacian smoothing under the assumption of smooth slip on adjacent patches, we utilize the SMC algorithm to solve for the two slip parameters on each patch. Additionally, we estimate the smoothing parameter and hierarchical scaling for each interferogram following the approach explained by Vasyura-Bathke et al. (2020).

In the static slip distribution step, there are a total of 1203 parameters to be determined based on the number of fault patches. This includes a pair of slip parameters for each segment (two slip components for 600 segments), one smoothing parameter, and a pair of hierarchical scaling parameters. The initial conditions for each Markov Chain within the high-dimensional parameter space are set by selecting a random smoothing weight from its prior (uniform) distribution. With these smoothing parameters, we compute the rake-parallel slip distribution using the regularized non-negative least-squares approach (Fukuda & Johnson 2008). We set the retrieved slip vectors in the slip-parallel direction and zero in the slip-perpendicular direction as the starting points for the first stage of the SMC process (Vasyura-Bathke et al. 2020).

The observed and modelled surface displacement maps, as well as the residual maps related to each interferogram, are shown in Fig. S12. The retrieved RMS values for distributed slip modelling are 1.6 cm for Int. 1, 1.8 cm for Int. 2 and 1.6 cm for Int. 3, showing an improvement compared to the uniform slip modelling. Fig. 7 presents the static slip distribution along the causative fault and the Marginal posterior distribution for the Laplacian smoothing weight, with the black vertical line marking the maximum likelihood value. The slip distribution is characterized by a significant slip area >0.5 m (red dashed polygon in Fig. 7) and a maximum slip area >3 m (continuous red polygon in Fig. 7) at the depths of ∼50–81 km and ∼60–75.5 km, respectively (see depth axis in the right-hand side of Fig. 7). Fig. 7 reveals that all reported main shock hypocentres, except CMT, lie to the SE of the maximum slip distribution area, and the slip pattern has a NE–SW orientation. The distribution of aftershocks from the ISC-EHB catalogue (circles in Fig. 7) aligns almost perfectly with the significant slip area. Considering a rigidity modulus of 30 GPa, the seismic moment is 4.84E + 20 Nm, equivalent to an earthquake with a M_w 7.8, in agreement with the uniform slip moment magnitude and the catalogue values (Table 2).

4.3.2. Kinematic slip distribution

In order to get the kinematic slip distribution of the 2013 Saravan earthquake, we utilize both InSAR interferograms and teleseismic waveforms as observations. We use the posterior slip distribution from the static slip distribution result as a prior slip to decrease the parameter space. Using the same velocity model and approaches in the static slip distribution, Green's functions for the joint inversion were computed for both the teleseismic waveforms and the InSAR data. In this step, the kinematic rupture parameters include the hypocentre location, time and the rupture velocity and duration on each fault patch. We also estimate the smoothing parameter and hierarchical residual scaling parameters. We set the initial rupture nucleation point from the uniform slip model, and we let it vary along the fault strike and dip based on the fault length and width in the inversion. We estimate the noise of each seismic trace five seconds before the synthetic P-wave arrival as a covariance matrix in the inversion. We consider four time windows four each subfault with a rise time of 5 s. We consider the posterior distribution from the static slip distribution as prior information for initializing the Monte Carlo samples. Fig. 8 represents the kinematic slip distribution results indicating that the hypocentre location is well constrained at the centre of the fault at a depth of 65 km (see star at Fig. 8a and depth axis in right) and propagated to the surface. The comparison of the slip patterns between the kinematic and static slip distributions reveals a concentration of maximal slip (>3 m) at a depth of 60–70 km close to the main shock's hypocentre. The source time function, which shows the rate of moment tensor release with time after earthquake origin, indicates that the earthquake lasted for 40 s (Fig. 8b). It is worth mentioning that the dashed red line (∼28 s) represents the end of the moment release based on the slip distribution pattern and rupture fronts. The synthetic waveforms of this kinematic slip distribution model fit the observed waveforms (Fig. S13).

5 DISCUSSION

In this study, we initially estimated the coseismic displacement offset of the 2013 Saravan Mw 7.7 earthquake using the GNSS time-series. The time-series reveal a maximum coseismic displacement at the SRVN station, the closest site to the epicentre, about 54.1 and 42.7 mm southeastward horizontal and upward vertical displacement, respectively (Fig. 2). Moreover, the GNSS data reveals far-field coseismic displacements extending up to 300 km from the epicentre, with stations like ZABL and NEHB in the north exhibiting particularly significant displacements. Although the GNSS measurements at stations farther from the epicentre were of smaller magnitude, the detection of far-field displacement underscores the critical need to incorporate deep earthquakes like this one into future seismic hazard assessments. Furthermore, these observations highlight the importance of deploying a denser network of GNSS stations within the region.

We generated coseismic interferograms using RADARSAT-2 images and utilizing newly provided orbital information by the Satellite operator (MDA Corporation) and updated GAMMA software. InSAR coseismic displacement maps show a maximum of ∼29 cm moving away from the satellite in the LOS direction (Figs 4a–c). We complete our data set by including broad-band teleseismic recordings to investigate the coseismic characteristics and source parameters of the 2013 Saravan M_w 7.7 earthquake. We use a Bayesian bootstrap optimization non-linear inversion approach to test NW and SE dipping fault scenarios. The retrieved geometry of the causative fault plane from joint inversion of GNSS, InSAR, and teleseismic data for our best-fitting model corresponds to a plane located ∼50–60 km southeastof the epicentre with a strike of 234° and a dip of 64° towards the northwest. It is characterized by a normal movement with a minor right-lateral strike-slip component compatible with USGS, CMT and GFZ focal mechanisms, and the results of Barnhart et al. (2014). To determine the spatial slip distribution of the Saravan earthquake and its temporal evolution, including the nucleation point and rupture duration, we apply a Bayesian finite fault inversion. We first use InSAR displacement maps to derive static slip distribution and then combine them with teleseismic P-waveforms to find the kinematic slip distribution. Both the static and kinematic slip distribution maps (Figs 7 and 8) reveal an oblique alignment pattern for the spatial slip distribution on the causative fault plane. The kinematic slip distribution constrains the maximum slip values of (>3 m) at the depth of 60–70 km (Fig. 8a). The kinematic slip distribution indicates that the rupture nucleates at a depth of 65 km, and according to the source time function (Fig. 8b), the earthquake reaches the maximum moment tensor release at ∼8 and ∼16 s. Both the kinematic slip distribution (Fig. 8) and the aftershocks from ISC-EHB (marked by black circles in Fig. 7) show a northeastward propagation. The slip distribution of Barnhart et al. (2014) using InSAR displacement map and teleseismic waveforms separately, considering NW dipping fault, depicts a top and bottom of significant (>0.5 m) and maximum (>3 m) slip values at 27–88 and 59–86 km (dashed and continues blue polygons in Fig. 7) which the bottom depth of it is located deeper than our slip pattern results. We relate this discrepancy to the newly processed interferograms used in our study and the higher residual of Barnhart et al. (2014) in their modelling.

To derive the driving mechanism for the Saravan earthquake, information about slab geometry is necessary for this part of the Makran subduction zone. However, this remains undocumented due to logistical constraints for instrumental deployments, particularly along the border between Iran and Pakistan. Barnhart et al. (2014) assumed a gentle dipping slab (∼0.5–8°) extending under a wide accretionary prism (400–500 km) and then a sharp dip variation (27–30°) in the district of the Saravan earthquake based on the earthquake locations of the USGS Catalog (http://www.earthquake.usgs.gov) and offshore seismic reflection profiles (White & Klitgord 1976). Thus, they concluded that the maximum slip patch location (see continuous blue polygon in Fig. 7) retrieved from InSAR and teleseismic analysis is deeper than the slab thickness of the oceanic crust by introducing a slab with a dip <30°. Based on this geometry and slip characterization on the seismogenic fault plane, Barnhart et al. (2014) proposed that the earthquake occurred in the mantle and slip propagated to or near the top of the oceanic crust.

The new seismic images of the deep structure of the western Makran subduction zone, using the local seismic stations (Motaghi et al. 2020; Priestley et al. 2022), suggest a slab geometry with a low-dip angle of 4° ± 2° near the coast which steepens to an 18° dip at the southern edge of the Jazmurian depression. It is worth mentioning that normal-faulting earthquakes on continents typically occur along planes dipping between 30° and 60° (Jackson 1987). By assuming a flat subducting slab and a normal faulting earthquake, we can combine the typical normal fault dip angle (around 45°) with the proposed 18° dip of the Makran subducting slab. This theoretical model predicts a total dip of approximately 63° for the overall fault plane, which closely matches the observed dip angle (64°) of the Saravan earthquake causative fault. Moreover, Penney et al. (2017) suggested that the occurrence of a normal faulting earthquake at a depth of 74 km beneath Bazman volcano, suggests that the subduction interface needs to increase its deep north of the line of normal-faulting earthquakes.

Taking into account the shallow thrust earthquakes (Fig. 9a) that align east–westwards dominantly along the shoreline and deeper normal earthquakes (Fig. 9a) to the north that occurred with a northeastward trend parallel to the volcanic arc (green triangles in Fig. 9a), we propose to extend this updated slab geometry towards the east up to the location of the Saravan earthquake with a northeastward trending of bending axis at the hinge in the subducting Arabian Plate, instead of a classical bending axis parallel to the trench. Indeed, this northeastward trend is also visible in the slip propagation of the Saravan earthquake (Fig. 9a), and extending the coseismic slip patch to the west generates an intersection with the seismic profile at the eastern edge of the Jazmurian basin where the slab bends beneath the seismic profile (Motaghi et al. 2020; Priestley et al. 2022).

Moreover, we propose the retrieved depth of the maximum slip pattern, ∼60–70 km (Fig. 8), as the rough value for the depth of the oceanic Moho (Fig. 9b) at the location of the Saravan earthquake in the central Makran subduction zone. This is compatible with the Saravan earthquake being generated by slab bending and with the bending axis of the Arabian oceanic slab being oblique along a NE–SW direction. This oblique bending probably also explains the oblique volcanic arc in the Makran subduction zone. This inference is also suggested by Penney et al. (2017) using the ∼ENE–WSW distribution of normal-faulting earthquakes at depths greater than 50 km. Finally, the NE-trend of the normal earthquakes (Fig. 9a) also points to a continuous NE–SW orientation of the bending axis at a hinge in the subducting Arabian Plate. Thus, we suggest a schematic slab geometry extrapolating the new findings in the western Makran subduction zone based on teleseismic data (Motaghi et al. 2020; Priestley et al. 2022) up to the central part of the subduction zone with a continuous slab characterized by a NE-trending bending axis (Fig. 9b).

Moreover, Barnhart et al. (2014) suggested that the Saravan earthquake nucleated within temperatures of 700–900 °C using the thermal model of Smith et al. (2013). Therefore, the earthquake should have occurred on a pre-existing fault due to slab-bending forces and dehydration of hydrous minerals (Barnhart et al. 2014). The new thermal model of Khaledzadeh & Ghods (2022) indicated that the previous model of Smith et al. (2013) for eastern Makran is incompatible with the relocated seismicity pattern of the Makran subduction zone. This discrepancy is caused by the steep dip of the subducting slab considered in the geometry of Smith et al. (2013). The new thermal model (Khaledzadeh & Ghods 2022) highlighted that the Saravan main shock occurs in the isotherm of 450–550 °C. Therefore, the dehydration embrittlement of the oceanic mantle, as suggested by Barnhart et al. (2014), is not required. Our slip distribution of the Saravan earthquake, new geometry information on the deep structure of the Makran subduction zone (Priestley et al. 2022), and the updated thermal model (Khaledzadeh & Ghods 2022) altogether underline that the earthquake occurred as a result of slab-bending forces with a dominant slip distribution within the oceanic crust.

Finally, the apparent stress (⁠${\tau }_a$⁠) of the Saravan earthquake can be estimated by dividing the total seismic energy released by the product of the causative fault area and slip (Choy & Kirby 2004). Considering the seismic energy of the Saravan earthquake of about 2.4E+15 J, according to the National Earthquake Information Center (NEIC), and the causative fault plane with a length of 80 km, width of 50 km and slip of 3 m (Fig. 8), the calculated apparent stress for the 2013 Saravan earthquake is 0.25 MPa. This apparent stress is consistent with a group of deeper earthquakes with ${\tau }_a<1\mathrm{M}\mathrm{P}\mathrm{a}$ (Choy & Kirby 2004). Choy & Kirby (2004) proposed this type of intraslab earthquakes occurring in warm-slab environments at considerable depths tend to release less energy per moment unit compared to those in cold-slab settings. Nevertheless, their occurrence inland from the trench and increased probability of being situated beneath densely populated areas should be recognized as potential predictors of seismic impact.

6 CONCLUSION

We studied the 2013 Saravan M_w 7.7 earthquake, the largest event in the last 50 yr in Iran, using both geodetic (GNSS and InSAR) and seismological observation. We observe a maximum coseismic displacement at the SRVN GNSS station with 54.1 and 42.7 mm southeast horizontal and upward vertical displacements, respectively. RADARSAT-2 coseismic displacement maps in ascending geometry reveal a NE-trending elliptical deformation pattern, reaching a maximum of ∼29 cm in the LOS direction of the satellite. The location of the aftershocks presents the same NE trend of the displacement towards the Iran–Pakistan border. The focal mechanism of the main shock is estimated using 25 broad-band teleseismic P-waveforms. We applied a joint uniform inversion on geodetic (InSAR and GNSS) and seismic data to evaluate the earthquake source parameters. The causative fault plane is a SW-striking, NW-dipping plane with dominant normal kinematics and a minor right-lateral strike-slip component. We retrieved the static and kinematic distributed slip pattern with a NE-trend at 50–81 km depth on the causative fault plane. The kinematic slip distribution reveals that the earthquake nucleates at the depth of 65 km, which aligns with the maximum slip concentration (>3 m) at the depth of 60–70 km and continues to release moment 28 s after the origin time.

The location of the Saravan earthquake close to the Iran–Pakistan border allows us to place constraints on the central Makran subduction geometry, where too few seismic stations are available due to logistical difficulties. The retrieved coseismic slip distribution in the central part of the Makran subduction zone in this study has the same northeastward trend as the distribution of preceding intraslab intermediate-depth normal earthquakes and the volcanic arc of this subduction zone. It proposes that the slab bending at the hinge of the Arabian subduction plate has a NE-trend that can explain the oblique volcanic arc in the Makran subduction zone. Thus, the suggested geometry can reproduce different seismotectonic characteristics and tectonic configurations in the Makran subduction zone. The precise analysis of the Saravan earthquake provides relevant new information on the Makran subduction zone geometry and its mechanism and fills an important data gap in central Makran. The new slab geometry will be helpful for the interpretation of surface deformation measurements and constraining a coupling map of the Makran subduction interface that is urgently needed for a realistic seismic hazard assessment in the region. A precise and well-constrained coupling map could reveal the real seismogenic potential of western part of Makran. Indeed, while the eastern part of the Makran subduction zone is able to produce large shallow thrust events (e.g. the 1945 M = 8 Pasni earthquake), none have been recorded in instrumental and historical times in the western Makran. Furthermore, the apparent stress of the strong Saravan earthquake aligns with a type of earthquake more likely inland and beneath populated areas with a high potential for seismic impact. This highlights the necessity of a comprehensive study based on seismic and particularly geodetic data to continue improving our knowledge of the geodynamics and the seismic hazard present in this area.

Acknowledgement

We thank the National Cartographic Center (NCC) of Iran for its precious help in conducting the GNSS measurements presented in our work. We would like to express appreciation to the Canadian Space Agency (CSA) for providing RADARSAT-2 images. We would like to thank the reviewers, Dr Camilla Penney and one anonymous reviewer, for taking the time and effort necessary to review the paper. We sincerely appreciate their insightful comments and suggestions, which helped us to improve the quality of the manuscript. We would like to express our appreciation to Hannes Vasyura-Bathke for his valuable insights and assistance in analysing the slip distribution through the utilization of the BEAT software. We are thankful to Mahtab Aflaki for her assistance in preparing the 3-D schematic illustration of the Makran subduction zone. This work was supported by the CNES (French Space Agency) through an APR project providing a research fellowship for Zahra Mousavi's research visit and the computational means for data processing. We would like to express our appreciation to the International Seismological Centre (2024), ISC-EHB data set (https://doi.org/10.31905/PY08W6S3).

DATA AVAILABILITY

We utilized teleseismic waveforms, obtained via FDSN (International Federation of Digital Seismograph Networks) services from IRIS (Incorporated Research Institutions for Seismology) and the GEOFON data centre of the GFZ German Research Centre for Geosciences. Earthquake data for our study was sourced from various catalogues, including ISC (International Seismological Centre) accessible at http://isc.ac.uk, CMT (Global Centroid Moment Tensor) available at https://www.globalcmt.org, USGS (United States Geological Survey) accessible at https://earthquake.usgs.gov and GFZ catalogue, as referenced in Quinteros et al. 2021. Most of the figures presented in this study were generated using the publicly available Generic Mapping Tools (GMT) software, developed by Wessel & Smith (1995). In our InSAR processing, we relied on the ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) Digital Elevation Model (DEM). This work uses the open-source libraries Bayesian Earthquake Analysis Tool (https://pyrocko.org/beat) and Pyrocko (https://pyrocko.org). KITE was used for post-processing the InSAR data (https://pyrocko.org/kite). The interferograms, GNSS (Global Navigation Satellite System) time-series, and relevant codes used for computation and plotting are available upon reasonable request. Interested parties may contact Meysam Amiri at m.amiri@iasbs.ac.ir to obtain access to these materials.

References