The 2019-2020 Khalili (Iran) Earthquake Sequence-Anthropogenic Seismicity in the Zagros Simply Folded Belt?

We investigate the origin of a long-lived earthquake cluster in the Fars arc of the Zagros Simply Folded Belt that is colocated with the major Shanul natural gas field. The cluster emerged in January 2019 and initially comprised small events of M n ∼ 3-4. It culminated on 9 June 2020 with a pair of M w 5.4 and 5.7 earthquakes, which was followed by >100 aftershocks. We assess the spatiotemporal evolution of the earthquake sequence using multiple event hypocenter relocations, waveform inversions, and Sentinel-1 Interferometric Synthetic Aperture Radar (InSAR) measurements and models. We find that the early part of the sequence is spatially distinct from the 9 June 2020 earthquakes and their aftershocks. Moment tensors, centroid depths, and source parameter uncertainties of 15 of the largest (M n ≥ 4.0) events show that the sequence is dominated by reverse faulting at shallow depths (mostly ≤4 km) within the sedimentary cover. InSAR modeling shows that the M w 5.7 mainshock occurred at depths of 2-8 km with a rupture length and maximum slip of ∼20 km and ∼0.5 m, respectively. Our results suggest that the 2019-2020 Khalili earthquake sequence was likely influenced by operation of the Shanul field, though elevated natural seismicity in the Zagros makes the association difficult to prove. Understanding how to distinguish man-made from natural seismicity is helpful for hazard and risk assessment, notably in the Zagros, which is both seismically active and rich in oil and gas reserves.

Introduction

Anthropogenic earthquakes, defined as those induced or triggered by human actions, have now been identified in many different regions across the globe (Foulger et al., 2018). Activities known or suspected to cause anthropogenic seismicity include subsurface fluid injection or extraction—through hydraulic fracturing, geothermal energy exploitation, and gas storage—as well as mining operations and water reservoir impoundment (Grigoli et al., 2017; Foulger et al., 2018; Keranen & Weingarten, 2018; Lei et al., 2020). These activities can introduce pore pressure transients and alter the local stress field, consequently promoting (or inhibiting) earthquake occurrence (Ellsworth, 2013; Dahm et al., 2013). Fluid injection‐induced seismicity (IIS) has become particularly widespread in recent years due to increased shale gas exploitation and wastewater disposal, geothermal stimulation, and gas storage (Ellsworth, 2013; Foulger et al., 2018; Lei et al., 2013). To date, IIS has reached moderate magnitudes—for example, the 2017 M w 5.5 Pohang earthquake (e.g., Grigoli et al., 2018), the 2016 M w 5.1 Fairview, M w 5.7 Prague, and M w 5.8 Pawnee, Oklahoma, earthquakes (Ellsworth, 2013; Keranen et al., 2014; Yeck et al., 2017), and the 2013 M w 4.3 earthquake at the Castor gas storage reservoir (Cesca et al., 2014)—sufficient that there are often strong socioeconomic impacts (Grigoli et al., 2017). Another example is the Sichuan Basin, where long‐term injection for disposal of unwanted water and short‐term injection for hydraulic fracking have induced significant earthquakes with magnitudes of up to M w 5.7 (Lei et al., 2020).

Recognizing anthropogenic earthquakes is particularly challenging in regions of naturally elevated seismicity with detailed source analyses essential in order to discriminate between the two (Dahm et al., 2015). The Zagros fold‐and‐thrust belt within the Arabia‐Eurasia collision zone (Figure 1a) offers an excellent example, comprising one of the most seismically active mountain belts as well as one of the greatest loci of oil and gas production in the world. The outer part of the range, known as the Simply Folded Belt, is characterized by a thick (averaging ∼10 km) sedimentary cover that contains hidden reverse faults that host frequent large, damaging earthquakes (e.g., Nissen et al., 2011; Talebian & Jackson, 2004). The folded and faulted sediments contain 90% of Iran's proven hydrocarbon reservoirs including the world's second largest gas reserves, estimated at ∼32.0 trillion cubic meters or ∼17% of Earth's total. Having started gas production in 1990, Iran now produces more than 1 billion cubic meters of gas per day from 36 gas fields, most of which are located in the Fars arc in the southeastern Zagros (Figure 1b; Esrafili‐Dizaji & Rahimpour‐Bonab, 2013; Vergés et al., 2011). Despite the intense hydrocarbon production, there have so far been no unequivocal cases of earthquakes linked to gas/oil extraction or wastewater disposal in the Zagros—though a few earthquakes have been attributed to reservoir impoundment (Kangi & Heidari, 2008), mining (Mansouri‐Daneshvar et al., 2018), and groundwater pumping (Kundu et al., 2019).

Figure 1

(a) Iranian seismicity, showing the location of the Zagros Mountains at the leading edge of the Arabia‐Eurasia collision zone. Red circles are M > 5.0 earthquakes from 1900 to 2019 from the USGS catalog. The most active, outer part of the Zagros (the Simply Folded Belt) can be subdivided into four tectonostratigraphic domains: from SE to NW, the Fars arc (F.A.), Dezful embayment (D.E.), Lurestan arc (L.A.), and the Kirkuk embayment (K.E.). (b) A zoom‐in of the Fars arc. A large number of anticlines are evident in the topography, several of which contain active gas fields (red ellipses). Black lines show major mapped active faults, including the right‐lateral Kazerun Fault (K.F.). Focal mechanisms from published waveform modeling studies are plotted at relocated epicenters and colored according to focal depth (Karasözen et al., 2019, and references therein). Notable earthquake sequences include those at Khaki‐Shonbe (KH‐SH; Elliott et al., 2015), Ghir (e.g., Berberian, 1995), Fin (Roustaei et al., 2010), and Qeshm (Lohman & Barnhart, 2010; Nissen et al., 2010, 2014). The black rectangle shows our study area (Figure 2).

Figure 2

(a) Topography, faults (modified from the GSI), background seismicity, the Shanul, and neighboring Homa, Varavi, and West Bavush anticlines. The white ellipsoid shows the spatial extent of the Shanul gas field. The inset photo shows the Shanul gas field from ICOFC. Black squares show the location of active wells in the Shanul, Homa, and Varavi gas fields. Red circles are relocated earthquakes before 2019 (Karasözen et al., 2019). Blue triangles show IRSC broadband stations (KHL1 and LMD1). Purple lines show locations of the seismic profiles (this study and Motamedi et al., 2012). The A–A’ and B–B’ seismic sections are presented in Figures S1 and S2 in Supporting Information S1. (b) The structural cross section across the Shanul and West Bavush anticlines (C‐C’ profile), constructed by integrating seismic, field, and remote‐sensing data.

Here, we investigate a prominent cluster of felt earthquakes near Khalili, in the central Fars arc, starting in January 2019. The swarm‐like activity and its spatial association with the major Shanul gas field raised legitimate concerns of an anthropogenic cause. The sequence culminated in mid‐2020 with an M w 4.7 earthquake on May 31, M w 5.4 and 5.7 earthquakes on 9 June, and a sustained aftershock sequence. The largest event, at 17:18 UTC on 9 June, was responsible for several injuries.

In this study, we present a detailed analysis of the Khalili sequence to gain insights into its mechanisms and origins. By utilizing stations of the Iranian Seismological Center seismic network (IRSC; see Data availability), which are denser here than in many other parts of Iran, we relocated the 18‐month‐long sequence and inverted moment tensors and centroid depths for 15 largest (M w > 4.0) events. We also estimated the coseismic slip distribution of the 9 June 2020 M w 5.7 mainshock using Interferometric Synthetic Aperture Radar (InSAR) measurements and elastic dislocation models. We compared the results with subsurface geology constructed using 2D seismic profiles. Our results reveal a close spatial correlation between seismicity and operations in the Shanul gas field as well as a number of anomalous source characteristics for the largest events. We suggest a possible case of anthropogenic earthquakes related to gas extraction in the Zagros.

Background

Active Tectonics, Structure, and Seismicity of the Fars Arc

The Fars arc refers to the arcuate part of the southeastern Zagros between the Kazerun fault in the west and the Bandar Abbas syntaxis in the east (Figure 1b). GPS measurements indicate 10 mm/yr of NNE‐directed convergence across the central Fars arc (e.g., Tatar et al., 2004). This shortening is manifested at the surface in symmetric, range‐parallel folds with amplitudes of up to a few kilometers and wavelengths of ∼10–20 km (e.g., Edey et al., 2020), and at depth in frequent earthquakes on steeply dipping (30°–60°), blind reverse faults (Berberian, 1995; Nissen et al., 2011; Talebian & Jackson, 2004). There are no known examples of coseismic surface rupture in the Fars arc, and the mechanical relationship between buried faults and surface folds remains a matter of debate.

The sedimentary cover of the Fars arc is detached from the underlying basement by a layer of Ediacaran–early Cambrian Hormuz salt, which also surfaces in numerous diapirs (e.g., Barnhart & Lohman, 2012; Edey et al., 2020; Jahani et al., 2009, 2017). Estimates of the depth of this interface vary from as little as ∼6–8 km (e.g., Sherkati et al., 2005) to as great as ∼14–20 km (Jahani et al., 2009, 2017). In the central Fars arc, closest to our study area, orogen‐scale geological cross sections interpret the basement depth to be ∼8–12 km (e.g., Allen et al., 2013; Najafi et al., 2014). Analysis of local and teleseismic earthquakes recorded in 1997 by a temporary (∼2 months) dense seismological network in the Ghir region, ∼100‐km north of Khalili, resolved thicknesses of 11 and 46 km for the sedimentary cover and crust, respectively (Tatar et al., 2004).

InSAR and teleseismic waveform modeling studies suggest that many of the larger (M w > 5) earthquakes of the Fars arc are located within the so‐called “Competent Group” of mechanically strong platform carbonates that make up the middle‐to‐lower sedimentary cover at depths of ∼5–10 km (Barnhart et al., 2013; Elliott et al., 2015; Lohman & Barnhart, 2010; Nissen et al., 2010, 2011, 2014; Roustaei et al., 2010). At the same time, a number of microseismic studies have indicated concentrations of small earthquakes at probable basement depths of ∼10–20 km (e.g., Nissen et al., 2011; Tatar et al., 2004). Helping to reconcile these differences, a recent relocation of the 70‐year catalog of well‐recorded, moderate‐to‐large earthquakes indicated a focal depth range of 4–25 km (Karasözen et al., 2019). Till now, the largest instrumental earthquakes in the Fars arc have not exceeded M w 6.7, reflecting that the seismogenic layer is segmented vertically by the Hormuz salt and other weak evaporitic or shale horizons within the cover, across which seismic rupture cannot propagate (Nissen et al., 2010). This mechanical segmentation also manifests itself in coseismic slip planes with characteristically narrow (small width‐to‐length ratio) dimensions (Elliott et al., 2015; Roustaei et al., 2010).

Geologic Structure, Production History, and Background Seismicity of the Study Area

The Shanul field is part of a concentration of natural gas reservoirs in the central and western Fars arc (Figure 1b). This region is characterized by symmetric to weakly asymmetric “whaleback” folds with characteristic wavelengths of ∼10–20 km and amplitudes of ∼2–4 km, which are controlled primarily by décollement along the Hormuz salt at ∼8–12 km depth (Allen et al., 2013; Motamedi et al., 2012; Najafi et al., 2014). Published seismic reflection imagery shows that many of the anticlines exhibit “pop‐up” geometries accompanied by pairs of opposite‐verging, high‐angle reverse faults on both flanks, originating either from the Hormuz décollement at the base of the cover (Najafi et al., 2014) or a secondary detachment within Triassic Dashtak evaporites of the middle cover (Figure 2) (Motamedi et al., 2012).

The Shanul gas reserves are contained beneath the broad, symmetric Shanul anticline, NW of the small settlement of Khalili (Figure 2). This anticline is outlined by resistant carbonates of the Miocene Mishan formation, while its close neighbor to the east—the West Bavush anticline—is expressed in the Oligocene Asmari limestone (Figure 2). A cross section of the Shanul anticline published by the Geological Survey of Iran (GSI) depicts the Shanul anticline as flanked by steep reverse faults that originate in Paleozoic strata in the anticline core. This view is supported by our own interpretation of newly available National Iranian Oil Company (NIOC) seismic reflection imagery (Figures S1 and S2 in Supporting Information S1). In contrast, the West Bavush anticline has a tighter and more asymmetric (southward divergent) shape, reflecting that its flanking reverse faults originate at shallower (∼3 km) depths in Triassic Dashtak evaporites (Figure 2 and Motamedi et al., 2012). The faults underlying both anticlines emerge at the surface as longitudinal reverse faults trending ∼N100o–105° (Figure 2). A combination of remote‐sensing, field, and seismic data permits us to construct a structural cross section across these anticlines that characterize the faulting and folding of the sedimentary cover (Figure 2b).

The Shanul reservoir was discovered in 1995 and the first well was drilled in 2004, with gas extraction starting in 2006 from Permo‐Triassic Dehram Group carbonates, capped by Dashtak evaporites at a depth of ∼3–4 km (Motamedi et al., 2012; Esrafili‐Dizaji & Rahimpour‐Bonab, 2013). The gas field belongs to the Iranian Central Oil Fields Company (ICOFC), one of the five major production companies of the NIOC, while the Southern Zagros Oil and Gas Production Company is responsible for its operation. So far, 18 wells have been drilled in the Shanul gas field (Figure 2). According to the ICOFC, 35 million cubic meters per day of gas are extracted from the Shanul field and the neighboring Homa reservoir, which together have a capacity of 220 billion cubic meters. Gas extraction from the 16th well of the Shanul field commenced in 2016, producing 600,000 cubic meters per day. We have no access to detailed production data, such as extraction or injection time series, volumes, or fluid pressures. We investigated whether InSAR time series data could reveal the timing of production activities by analyzing Sentinel‐1 time series from 2018/01/01 to 2020/06/01 using the LiCSBAS software (Morishita et al., 2020). Unfortunately, prior to the 2020 mainshock, we could not distinguish any clear surface deformation signals from above the level of atmospheric noise; a full analysis would be required to extract any signal from the time series but lies beyond the scope of this study. However, considering the overall gas capacity and extraction rate, gas reserves from both reservoirs are likely to become depleted within about 3 years. Fluid injection, which is typically applied in the aging gas fields when production is waning, is therefore likely to have started in both reservoirs.

There are no historical records of any earthquake unambiguously linked to faults within our study area (Ambraseys & Melville, 1982; Berberian, 1995). Modern seismicity in Iran is monitored and reported by permanent networks of the IRSC and the International Institute of Earthquake Engineering and Seismology (IIEES), which have been densified with time, especially since around 2012. Before this time, seismic coverage in the Zagros was restricted to sparse broadband stations (Hosseini et al., 2019), limiting route detection thresholds and location accuracies and preventing focal mechanisms from being determined for most small‐to‐moderate earthquakes. Despite this, global networks have provided a magnitude of completeness of less than 4.5 for the period after 1990 in Zagros (Shahvar et al., 2013). Nevertheless, the relocated catalog of Karasözen et al. (2019) indicates only two events of m b 4 (on 10 August 2009 and 1 October 2010) that are colocated with the Shanul anticline (Figure 2).

Source Characteristics of the 2019–2020 Khalili Seismic Sequence

Multiple‐Event Relocation

We start by assessing the overall spatiotemporal evolution of the 2019–2020 sequence using a multiple‐event calibrated hypocentral relocation. We used the Mloc implementation (Bergman & Solomon, 1990) of the hypocentral decomposition algorithm (Jordan & Sverdrup, 1981), consistent with several earlier regional studies (Elliott et al., 2015; Karasözen et al., 2019; Nissen et al., 2010, 2019; Roustaei et al., 2010). Among the ∼300 events (M n ≥ 2.5) reported by IRSC, we relocated 115 events (M n ≥ 3.0). IRSC station coverage is sufficient (Figures S3 and S4 in Supporting Information S1) that we could employ a “direct calibration” (Karasözen et al., 2019), yielding epicentral uncertainties of ≤5 km for most events and ≤3 km for 84 of the best‐recorded earthquakes (Figure S5 in Supporting Information S1). All the events in the cluster have a maximum Pg azimuthal gap of less than 120° (Figure S4 in Supporting Information S1), ensuring a strong control on epicenter locations. The rather thick crust, 43 km, ensures that the Pn phases beyond 200 km, so all phases within 1.8° are still Pg arrivals. The hypocentroid has been calculated only using good‐quality Pg and Sg phases recorded at less than 1.8°, which were sufficient to avoid azimuthal gaps greater than 100°.

We used a slightly modified version of the 1D‐layered velocity model (Figure S6 in Supporting Information S1) of Karasözen et al. (2019) to predict theoretical travel times (Figure S7 in Supporting Information S1) at local and regional distances. Owing to insufficient close by station coverage, we were unable to solve the focal depths of most events. The focal depths of 19 events were constrained with phase reading from a very nearby station (KHL1 station, Figures 2 and 3 panels a), but focal depths for the remaining 96 events were fixed to a representative cluster average of 7 km (Table S3 in Supporting Information S1). Therefore, the relocated seismicity cannot be used to infer the dip of the causative faulting at depth. From experience, fixing the focal depths of shallow (<15 km) earthquakes in this way has a negligible effect on epicenter accuracy (Ghods et al., 2012; Karasözen et al., 2019). In other words, actual focal depths could be several kilometers deeper or shallower than the cluster default depth of 7 km, and this would not influence our obtained epicenter locations.

Figure 3

(a) Relocated epicenters of M n ≥ 3.0 events and focal mechanisms of M n ≥ 4.0 events from January 2019 to October 2020 colored by centroid depth. Black circles are events in Phase 1 (prior to 31 May 2020) and gray circles are those in Phase 2. Red lines show faults in the region (modified after the GSI). The white ellipsoid shows the spatial extent of the Shanul gas field and black squares show the location of active wells in the Shanul gas field. (b) Structural cross section across the Shanul and West Bavush anticlines (C–C’ profile and color same as Figure 2) with our relocated focal mechanisms at their centroid depths. (c) Temporal evolution of seismicity from IRSC catalog (M n ≥ 2.5) with events plotted by the magnitude and colored as in (a).

The spatiotemporal evolution of seismicity clearly depicts two phases of the sequence (Figure 3). Phase 1 started in January 2019 and continued through early 2020 and is swarm like, lacking a dominant mainshock or clear taper of aftershocks. Phase 1 events follow a WNW–ESE‐oriented trend centered on the southern limb of the Shanul anticline. Phase 2 commenced with the M w 4.7 foreshock on 31 May 2020 and included the 9 June 2020 M w 5.4 and 5.7 earthquakes and their aftershocks. Phase 2 seismicity lies along a separate WNW–ESE‐oriented trend located between the Shanul and West Bavush anticlines (Figure 3).

Furthermore, we calculated the local magnitude of completeness as M c = 2.9 using the IRSC catalog for 2008–2021. Therefore, by restricting our analysis to events with M n > 3.0, we avoided any possible temporal change of completeness magnitude (Figures S8 and S9 in Supporting Information S1). The temporal variation of the M c shows that it starts from 4.5 in 1996 but reduces to ∼3.0 in 2008 (Figure S9 in Supporting Information S1).

Focal Depth of the 9 June 2020 M w 5.7 and 5.4 Earthquakes

Well‐constrained focal depths are an important potential discriminator of induced earthquakes. We used independent teleseismic data to estimate the focal depth of the two largest events of the Khalili sequence, that is, the M w 5.7 and M w 5.4 earthquakes on 9 June 2020.

We modeled the delay between the direct P arrival and the surface‐reflected pP phases at teleseismic distances recorded at different arrays, which depends on the source depth and the average P‐wave velocity above the hypocenter. For this analysis, we used the Array Beam Depth Tool (see Code Availability). The qualitative comparison of the observed and synthetic beams for different source depths (Figure 4) allows us to estimate focal depths, complementing the other observations made in this study. The synthetic beams are constructed assuming source‐ and receiver‐specific crustal models (after Karasözen et al., 2019, and CRUST2.0 model, Bassin et al., 2000) and a global AK135 model (KennettEngdahl, Buland, & Buland, 1995) for the propagation of the seismic waves in between. This method is extensively used to estimate the focal depth of earthquakes in several other regions (e.g., Braun et al., 2018; Büyükakpınar et al., 2021; Cesca et al., 2021; Gaebler et al., 2019; Negi et al., 2017). Results from all of the analyzed arrays are consistent with a focal depth of ∼4–5 km for both the M w 5.7 earthquake (Figure 4) and the M w 5.4 earthquake (Figure S10 in Supporting Information S1).

Figure 4

Estimation of the focal depth of the 9 June 2020 M w 5.7 Khalili mainshock using teleseismic records in six seismic arrays. (a) Yellow triangles show the location of the analyzed arrays (BCA, IMAR, YKA, ILAR, GERES, and Alice). The black star shows the 9 June 2020 M w 5.7 Khalili earthquake. (b) Black lines show synthetic waveforms based on bespoke velocity models and source mechanisms fixed at varying depths. Red lines show the depth phases (P and pP). Blue waveforms represent observed stacked array beams corresponding to each array. A focal depth of 4–5 km offers the best visual coherency between observed and synthetic traces.

Regional Moment Tensor Solutions

Full moment tensor (MT) solutions obtained through regional waveform inversions are a key tool for induced seismicity studies, providing critical information on the source geometry and the rupture process (e.g., Dahm et al., 2015). Observations of relevant non‐double‐couple (non‐DC) components through MT decomposition have been used as an indicator for seismicity induced by mining (Cesca et al., 2013) and fluid injection (Zhang et al., 2016). For very specific earthquakes, for example, those involving an underground collapse, a full MT inversion can be directly used to detect specific induced events (e.g., Cesca et al., 2013). The presence of intrinsic non‐DC terms for IIS and their interpretation in terms of rupture processes remains debated (Lei et al., 2017, 2020). While non‐DC terms have been repeatedly reported for IIS (Dost et al., 2020; Wang et al., 2018; Zhao et al., 2014; Zhang et al., 2016), they were rarely confirmed for the largest IIS earthquakes (Grigoli et al., 2017). In the frame of geothermal stimulation, events with significant isotropic components have been found for earthquakes occurring during the injection at Basel, Switzerland, but not for those occurring after the stimulation (Zhao et al., 2014). At the Groningen gas field, the Netherlands, where gas was extracted over many years, the seismicity accompanying the reservoir depletion exhibits both DC terms, controlled by local fault structures, and predominantly negative isotropic non‐DC terms, interpreted as a consequence of compaction processes (Dost et al., 2020). A limited, but consistent amount of non‐DC terms have been resolved for IIS in Alberta, Canada, interpreted as a consequence of fracture growth and/or complex noncoplanar faulting accompanying hydraulic‐fracturing stimulations (Wang et al., 2018). Similarly, at Pohang, South Korea, an M w 5.5 earthquake induced during geothermal stimulation showed a clear non‐DC component, which could be explained by the quasi simultaneous activation of two faults with different geometries and rupture types (Grigoli et al., 2018). Perhaps, the clearest evidence is provided by Zhang et al. (2016), who showed that induced seismicity in Western Canada exhibits significantly higher non‐DC moment tensor components than natural earthquakes when estimated using the same procedures and velocity models and who specifically ruled out artifacts due to seismic noise or simplified velocity models. However, most triggered earthquakes are characterized by shear fracturing on preexisting faults, and full MT inversions and decomposition results are more useful for the inference of the rupture geometry than for discrimination. Furthermore, full MT inversions are challenging for small‐to‐moderate magnitude events, requiring a robust assessment of any resulting non‐double‐couple (DC) source terms. Probabilistic waveform inversion techniques, which provide estimations of the parameter uncertainties and trade‐offs, provide the best approach to assess reliable non‐DC components (Kühn et al., 2020; Zahradnik et al., 2008). Among the parameters that are estimated by a centroid full MT inversion (scalar moment, centroid depth, fault plane angles, and percentages of decomposed MT terms), the centroid depth is a particularly important discriminator between anthropogenic and natural seismicity in the region where both are probable (Dahm et al., 2013; Grigoli et al., 2017).

We performed full MT inversions for four earthquakes in the Khalili sequence; the M w 5.4 and 5.7 Phase 2 events on 9 June 2020 and two events with normal mechanisms in Phase 1. We also undertook deviatoric MT inversion—using the standard decomposition between compensated linear vector dipole (CLVD) and double‐couple terms—for 11 moderate events of M w 4.0–4.8, including two Phase 1 events and seven additional Phase 2 events. We used a probabilistic MT inversion method (Heimann et al., 2018), which provides ensembles of best‐fitting MTs, which are used to estimate uncertainties and trade‐offs for all inverted source parameters. This technique has been successfully applied to other earthquakes in the Zagros (Jamalreyhani et al., 2019) and in other regions (e.g., Kühn et al., 2020; Petersen et al., 2021).

We set up the MT inversion to simultaneously fit three‐component waveforms in the time (full displacement waveforms) and in the frequency domains (full amplitude spectra). Synthetic seismograms were computed using precalculated Green's functions (Heimann et al., 2019) based on a regional velocity model by Karasözen et al., (2019). For events smaller than M w 5, we adopted the frequency band of 0.02–0.07 Hz; for the larger pair of events, we used the frequency band of 0.015–0.05 Hz. To avoid systematic error in the MT solutions due to sensor misorientation, we applied the sensor orientation corrections (Braunmiller et al., 2020) for the IRSC stations. The resolved focal mechanisms are in good agreement with the IRSC, GCMT, and GEOFON solutions, for the few cases where they are available, but we estimate in most cases shallower centroid depths (mostly ≤4 km with estimated uncertainties of 0.5 km) noting also that Global CMT has a very limited resolution at shallow depths of 15 km or less (Maggi et al., 2002; Wimpenny & Watson, 2021). The centroid depths obtained by IRSC using completely different methods also highlight the shallow depths. All obtained source parameters together with their uncertainties (68% confidence intervals) are listed in Table 1, and depth probability distributions for each event studied are shown in Figure S11 in Supporting Information S1.

Table 1

Moment Tensor Solutions of the 15 Events in the 2019–2010 Khalili Seismic Sequence Obtained in This Study

No	Date and time (UTC)	Latitude°	Longitude°	M w	Depth (km)	Strike1°	Dip1°	Rake1°	Strike2°	Dip2°	Rake2°
1	2019‐06‐24 15:14:08	27.677	53.231	4.2	4.0 ± 1.0	244 ± 13	62 ± 3	−85 ± 16	53 ± 14	28 ± 5	−99 ± 20
2	2019‐06‐28 09:08:54	27.594	53.175	4.2	6.0 ± 1.0	335 ± 24	83 ± 2	170 ± 46	66 ± 2	80 ± 2	7 ± 2
3	2019‐07‐16 12:02:24	27.686	53.284	4.0	4.0 ± 2.0	265 ± 56	64 ± 7	−56 ± 25	28 ± 14	41 ± 12	−138 ± 42
4	2019‐11‐13 17:57:45	27.737	53.073	4.2	1.0 ± 0.5	6 ± 20	37 ± 5	163 ± 4	110 ± 5	80 ± 3	53 ± 4
5	2020‐05‐31 23:59:00	27.756	53.309	4.7	3.0 ± 0.5	292 ± 5	58 ± 4	99 ± 5	96 ± 5	33 ± 3	76 ± 8
6	2020‐06‐09 16:08:48	27.704	53.363	5.4	3.0 ± 0.5	278 ± 5	46 ± 3	90 ± 8	98 ± 6	44 ± 3	90 ± 7
7	2020‐06‐09 17:18:12	27.669	53.411	5.7	3.0 ± 1.0	97 ± 7	58 ± 6	79 ± 8	297 ± 8	33 ± 6	107 ± 12
8	2020‐06‐13 22:04:14	27.691	53.361	4.8	3.0 ± 0.5	294 ± 4	33 ± 2	100 ± 6	102 ± 4	57 ± 1	83 ± 4
9	2020‐06‐13 23:15:03	27.701	53.342	4.6	3.0 ± 0.5	288 ± 9	43 ± 6	85 ± 13	115 ± 9	47 ± 6	95 ± 12
10	2020‐06‐14 18:06:00	27.698	53.373	5.2	2.0 ± 0.5	296 ± 10	35 ± 11	110 ± 14	92 ± 11	57 ± 10	77 ± 15
11	2020‐07‐02 00:29:00	27.739	52.885	4.1	1.0 ± 0.5	304 ± 6	43 ± 10	100 ± 9	112 ± 8	48 = 11	82 ± 8
12	2020‐07‐10 20:14:04	27.662	53.525	4.5	4.0 ± 0.5	276 ± 5	26 ± 2	83 ± 8	104 ± 4	64 ± 2	93 ± 4
13	2020‐08‐25 12:16:00	27.788	53.213	4.2	3.5 ± 0.5	280 ± 6	48 ± 4	90 ± 7	100 ± 6	41 ± 4	90 ± 8
14	2020‐08‐31 03:36:50	27.810	53.222	4.8	3.0 ± 0.5	284 ± 3	45 ± 2	91 ± 2	103 ± 3	45 ± 2	89 ± 3
15	2020‐09‐08 01:34:17	27.791	53.246	4.3	1.5 ± 0.5	305 ± 9	77 ± 11	144 ± 21	44 ± 20	55 ± 11	15 ± 21

Note. Columns refer to the event number, date and time in UTC (yyyy‐mm‐dd hh:mm:ss), relocated latitude and longitude, magnitude (Mw), centroid depth, and strikes, dips, and rakes of the two nodal planes with estimated uncertainties.

We observe distinct patterns of focal mechanism in the two phases of the 2019–2020 seismic sequence (Figure 3a). Phase 1 events exhibit diverse mechanisms and depths, comprising one very shallow (1 km centroid depth) reverse faulting event, two normal faulting events at 3–4 km, and a slightly deeper (7 km) strike‐slip earthquake. The normal faulting events appear linked to a series of short, shallow, ∼NE‐trending faults mapped along the crest of the Shanul anticline (Figures 2 and 3), likely the consequence of bending stresses within the upper layer of the fold. Figure S12 in Supporting Information S1 shows waveform and amplitude spectra fits for the 24 June 2019 M w 4.2 event with its normal mechanism. We observe that for the pair of normal mechanisms, the non‐DC part is larger than the DC part (Figures S13–S14 in Supporting Information S1).

In contrast, Phase 2 seismicity follows a typical foreshock‐mainshock‐aftershock pattern. Its similar ENE–WSW‐oriented thrust faulting mechanisms and consistently shallow (≤4 km) centroid depths suggest rupture along a single fault or fault zone, parallel to both the local fold axes and the overall seismicity trend. Figure 5 shows waveform and amplitude spectra fits for the 9 June 2020 M w 5.7 mainshock. Our full moment tensor decomposition into isotropic (ISO) and CLVD and DC components reveals a relatively large CLVD component, similar to that resolved independently by GCMT and GEOFON. The M w 5.7 mainshock centroid depth is slightly shallower (3  1 km) than the focal depth (∼5 km, resolved by teleseismic pP‐P delays in Section 3.2), suggesting upward rupture directivity. In any case, both results point to a shallow source within the middle‐to‐upper sedimentary cover.

Figure 5

Full moment tensor solution of the 9 June 2020 M w 5.7 earthquake. (a) Waveform fits in time domain and amplitude spectra for the M w 5.7 earthquake. Red and gray waveforms/spectra show synthetic and observed records, respectively. Information on the top of the waveform fits gives station names with transverse (T), radial (R), or vertical (Z) components. Numbers within the panels describe the time window and the frequency band. (b) The centroid depth probability distribution function of the M w 5.7 event from waveform inversion. The red solid vertical line gives the median of the distribution and the dashed red vertical line gives the mean value. Dark gray vertical lines show reference parameter values (8 km). The overlapping red‐shaded areas show the 68% confidence intervals (innermost area), the 90% confidence intervals (middle area), and the minimum and maximum values (widest area). The plot ranges are defined by the given parameter bounds and show the model space (1–15 km). (c) The decomposition of the full moment tensor in ISO, CLVD, and DC parts. The symbol size indicates the relative strength of the components. The Global Centroid Moment Tensor solution is shown for comparison. (d) The fuzzy full MT solution illustrates the uncertainty of the solution.

By implementing a Bayesian bootstrap approach in this study, the used optimization method explores the full model space and maps model parameter uncertainties and trade‐offs. Figure S15 in Supporting Information S1 illustrates the range of solutions for each parameter pair of the M w 5.7 mainshock, which helps to identify potential trade‐offs. No significant trade‐offs are found among the depth and other source parameters, which confirm the reliability of the depth estimate and overall of the moment tensor solution. We also fixed the best moment tensor solution and magnitude and compared synthetic seismograms for different centroid depths with the observation for the M w 5.7 mainshock. Figure S16 in Supporting Information S1 shows how the waveforms fit changes as a function of depth. The best fit for the whole waveform is found for the shallowest source (3 km) for the M w 5.7 mainshock (Figure S16 in Supporting Information S1).

Fault Geometry and Slip Distribution of the 9 June 2020 M w 5.7 Earthquake From InSAR Modeling

We used Sentinel‐1 InSAR imagery and elastic dislocation modeling to characterize the 9 June 2020 M w 5.7 mainshock fault geometry and coseismic slip distribution. Using the earliest available postseismic acquisitions, we constructed one 12‐day interferogram on descending‐track D64 and two 24‐day interferograms on ascending tracks A130 and A28 (Figure 6a). The interferograms also each capture the M w 5.4 foreshock that occurred 70 min before the mainshock. All three interferograms exhibit a WNW–ESE‐oriented pattern of deformation, containing 3–4 fringes, equivalent to ∼10 cm of displacement toward the satellite. The close similarity of fringe patterns in the descending‐ and ascending‐track interferograms implies that this deformation is predominantly uplift, centered upon a tight syncline between the Shanul and West Bavush anticlines. A single fringe (∼3 cm) of deformation away from the satellite is also evident to the North of the main fringe ellipse, colocated with the northern limb of the West Bavush anticline.

Figure 6

(a) From top to bottom: coseismic interferograms on tracks D64, A130, and A28. From left to right: observed, model, and residual interferograms. Results are shown rewrapped in order to accentuate deformation gradients (2π radians = 2.77 cm displacement). The black line is the surface projection of the model faults. Red and blue stars are relocated epicenters of the M w 5.7 and M w 5.4 earthquakes, respectively, and white dots are relocated aftershocks (Phase 2). (b) Coseismic slip distribution from modeling the InSAR data. The model fault is divided into 2 km2 patches. Red and gray stars are hypocenters of the M w 5.7 and M w 5.4 earthquakes, respectively. Circles show the relocated 20 aftershocks, the focal depths of which were constrained with phase reading from a very nearby station and colored according to time.

To characterize the fault geometry and slip distribution responsible for the observed fringe patterns, we followed routine elastic dislocation modeling procedures (Funning et al., 2005; Okada, 1985; Wright et al., 1999); that have been applied to other earthquake sequences in the Fars arc (e.g., Elliott et al., 2015; Nissen et al., 2010; Roustaei et al., 2010); full details are provided in Supplementary text S1. We found that both NNE‐ and SSW‐dipping model faults could reproduce the overall InSAR deformation pattern, but that in either case, a kinked, two‐fault geometry offered noticeable visual and numerical improvements over a single‐slip plane. Based on the adjacency of the mainshock hypocentral locations as well as the narrow aftershock cloud, our preferred configuration is the SSW‐dipping, two‐fault model (Figure 6a and 6b), though alternative models are also provided in Figures S17–S21 in Supporting Information S1 . We also produced a variety of model slip distributions across a range of Laplacian smoothing factors (Funning et al., 2005), selecting for our preferred model one with a relatively high degree of smoothing but a relatively low model misfit (Figure S22 in Supporting Information S1).

The preferred model faulting parallels the West Bavush anticline but projects to the surface close to its axial trace and is thus offset southward from the SSW‐dipping fault interpreted in Figure 3b. The ∼15 km‐long western InSAR model fault segment strikes 108°, dips 64° N, and has a rake of 84°; the shorter ∼5‐km‐long eastern model fault segment strikes 95°, dips 66° N, and has a rake of 110°. Near the mainshock hypocenter, slip is concentrated at depths of ∼2–8 km, but elsewhere it falls within an even narrower range of ∼2–6 km (Figure 6b), consistent with the 3 1 km centroid depth resolved using seismological data and the independent depth estimate based on array data at teleseismic distances. The hypocenter is located close to the bottom of the slip patch and at about midway along its length, indicating rupture propagation updip and bilaterally along strike. The upward directivity may be responsible for the extensive damage in the hanging wall of the mainshock fault and could help explain the widespread occurrence of land sliding evident in satellite photographs (Valkaniotis, 2020).

The InSAR model moment of ∼6.6 × Nm (M w 5.8) is roughly 50% larger than that of the mainshock GCMT and USGS W‐phase solutions (4.2–4.3 × Nm, M w 5.7). The M w 5.4 foreshock can account for much of the difference, though it is possible that our models also include a small amount of postseismic afterslip.

To test the robustness of our inversion results, we performed two resolution tests following Elliott et al. (2015); full details are provided in supplementary Text S2. In the first, we tested the sensitivity of residual (observed minus modeled) displacements to variations in the top and bottom depths of a rectangular, uniform slip model fault plane. Results show that the top and bottom depths are well resolved at 2–3 km and 4–6 km, respectively (Figure S23 in Supporting Information S1), in good agreement with our preferred model slip distribution (Figure 6b). In the second test, we performed check‐board inversions to assess the slip distribution resolution (Elliott et al., 2015). Results show that our model slip distribution can be characterized at a 2‐km resolution along the updip edge of the slip plane, but that the deeper part of the model slip distribution might be subject to smoothing artifacts (Figure S24 in Supporting Information S1).

Discussion: Was the 2019–2020 Khalili Sequence Induced?

Observations that unequivocally link a seismic sequence to fluid injection time series in a gas field are rare. Here, we have no access to the kinds of production data—such as extraction/injection time series, volumes, and pressures of the injected fluid—that could confirm a causative link to trends in the temporal evolution of seismicity. However, the timing of the 2019–2020 Khalili sequence is at least consistent with an anthropogenic origin. The first well in the Shanul gas field was drilled in 2004, and gas extraction started in 2006. However, fluid injection started later, only after production had peaked. Seismicity in the Shanul gas field emerged several years after the start of extraction. It is reported that when injecting wastewater into a depleted gas reservoir, significant seismic activity would be initiated and it rapidly increases if the pressure exceeds the original pressure before production (Lei et al., 2013). However, time delays of several months or even several years between the start of production and the emergence of induced seismicity have been observed in other gas fields. For example, in the Groningen region of the Netherlands, the first earthquake occurred after 28 years of production (Richter et al., 2020), while earthquake sequences near the Oklahoma Wilzetta and Texas Cogdell oil fields commenced ∼20 years after the initiation of fluid injection (Keranen et al., 2014). In the Sichuan Basin, China, induced seismicity emerged several decades after the start of fluid injection (Lei et al., 2020).

The Fars arc of the Zagros is one of the world's most seismically active mountain belts in which M 5–6 mainshock‐aftershock sequences are common. Figure 1 demonstrates how seismicity is distributed broadly over a width of ∼100 km, that many earthquakes are colocated with prominent anticlines, and that in the central/western Fars arc, many of the anticlines contain gas fields. Therefore, on the basis of their location alone, the Khalili earthquakes do not appear to be unusual. We can compare source characteristics of the 2019–2020 Khalili sequence with those of natural, background seismicity in the Fars arc. Our calibrated relocation of the seismic sequence reveals two distinct clusters in space and time. The first phase is spatially localized in the southern part of the Shanul gas reservoir and resembles a swarm without a clear mainshock (Figure 3c). Three out of four centroid depths (1–4 km) are close to the probable production levels of the reservoir (∼3–4 km), and they differ markedly from earthquake depths in neighboring parts of the Zagros, which average ∼8–10 km for larger events modeled using teleseismic body waveforms and ∼10–15 km for smaller events recorded locally (Figure 7; e.g., Nissen et al., 2011; Tatar et al., 2004).

Figure 7

Earthquake depth histograms for the Fars arc (between 51.5° and 56.7° E and up to 29.5° N). (a) Centroid depths from waveform modeling. Red bars show the 2019–2020 Khalili sequence from regional waveform modeling (this study). Gray bars show teleseismic body‐waveform models of earlier earthquakes (Nissen et al., 2011, and references therein); note that there are no teleseismic centroid depths shallower than 4 km. (b) Focal depths from local and regional arrival time data. Red bars show calibrated relocations with free depth solutions for the 2019–2020 Khalili sequence (this study) with number of events multiplied by 10 for clarity. Gray bars show microseismicity recorded by dense local networks at Qeshm Island (Nissen et al., 2010), Fin (Roustaei et al., 2010), Ghir, Khurgu (Nissen et al., 2011), and Khaki‐Shonbe (Elliott et al., 2015).

Seismicity during the first phase also exhibits a wide variety of focal mechanisms, which are atypical in the region, including two normal faulting events with significant non‐DC components within the area of the Shanul gas reservoir. We interpret that the first phase was induced by the depletion of the Shanul gas reservoir, which likely produced local stress perturbations and/or pore pressure changes that could have triggered slip on preexisting faults (such as those identified in regional geological maps and seismic lines). Normal faulting events are favored to occur above a depleting reservoir (Segall, 1989). The second phase of the sequence commenced with the M w 4.7 foreshock on 31 May 2020, and culminated in the M w 5.4 and 5.7 earthquakes on June 9. This phase is centered upon the West Bavush anticline, ∼15‐km NE of the Shanul reservoir. Similar to the first phase, it is marked by shallow centroid depths of ≤4 km (Figure 7). However, unlike the first phase, the second phase exhibits typical foreshock‐mainshock‐aftershock patterns and reverse faulting mechanisms that are compatible with regional tectonic stresses. The 9 June 2020 M w 5.7 mainshock ruptured a steep, SSW‐dipping reverse fault within the core of the west Bavush anticline, releasing most of its moment at shallower depths of ∼2–6 km, within the upper half of the sedimentary cover. The fault width (∼4 km) is thus narrow with respect to the fault length (∼20 km), similar to patterns observed elsewhere in the Fars arc and presumably reflecting lithologic barriers to updip and downdip rupture propagation (Elliott et al., 2015; Roustaei et al., 2010).

Most well‐recorded Fars arc earthquakes of similar magnitude to the Khalili mainshock—such as those at Qeshm Island in 2005–2008, Fin in 2006, and Khaki‐Shonbe in 2013—were centered in the middle or lower cover (at depths of ∼5–9 km), and their aftershock sequences included concentrations of events at unequivocal basement depths (Barnhart et al., 2013; Elliott et al., 2015; Lohman & Barnhart, 2010; Nissen et al., 2010, 2014; Roustaei et al., 2010). The shallow depths of Phase 2 of the Khalili sequence are therefore unusual (Figure 7), but they are not unprecedented. The 2013 M w 6.2 Khaki‐Shonbe earthquake ruptured a subsidiary fault plane (southeast of the deeper, main rupture) at shallow depths of ∼2–4 km (Elliott et al., 2015), and the nearby 2014 M w 5.1 Bushkan earthquake slipped at depths of ∼2–6 km (Kintner et al., 2019). We note that both these shallow slip planes are associated with anticlinal structures containing active gas fields (the Kangan/Zireh fields and Dalan fields, respectively), though further study would be needed to assess whether they too might have been induced.

The depths of Phase 2 events are also similar to those reported for induced seismicity associated with hydrocarbon reservoirs in other regions (e.g., Cesca et al., 2014; Cesca et al., 2021; Dahm et al., 2015). We therefore interpret that while the second phase involved the release of background tectonic stresses, it might have been triggered by stress changes from the first stage or by fluid migration, by means of pore pressure diffusion or poroelastic stresses. Fluid migration may reach large distances of tens of kilometers (e.g., Goebel et al., 2017), which makes this hypothesis fully compatible with the location of Phase 2 events and can occur over relatively long time periods. The timing of such a triggering process is difficult to interpret here, due to the lack of knowledge on local diffusivity conditions and the potential presence of pathways controlling fluid migration. However, fluid migration of ∼3–4 km per month has been reported in other gas reservoirs (e.g., Cesca et al., 2021) and a delay of 4–5 months between phases 1 and 2 at ∼10 km distance appears to be consistent with previous observations of fluid‐driven seismicity (e.g., Hainzl et al., 2012).

The lack of subsurface fluid flow data and or constraints on geomechanical properties precludes us from quantifying the likelihood of induced or triggered seismicity through physics‐based probabilistic modeling (Dahm et al., 2015). However, as a check on our interpretation of the 2019–2020 Khalili sequence as being of anthropogenic origin, we applied qualitative discrimination approaches based on a series of questions. This qualitative discrimination approach has been used for the discrimination of induced seismicity in various other regions (e.g., López‐Comino et al., 2021). Application of the Frohlich et al. (2016) criteria supports our inference that the Khalili earthquake sequence is induced. Moreover, we applied a new framework proposed by Verdon et al. (2019), comprising a series of variably weighted questions with positive numerical scores assigned to characteristics of induced seismicity (+100%) and negative scores to those of natural events (−100%). Our seismological analysis yields an induced assessment ratio (IAR) of +40% with a high evidence strength ratio (ESR) of 94.8% (ranging between 0% and 100%) supporting the quality and quantity of information used in the assessment (Figures S25 and S26 in Supporting Information S1). Results indicate that the 2019–2020 Khalili seismic sequence might have been induced. Our interpretation is based on (a) the occurrence of swarm‐like activity very close to the active gas field (Phase 1 being inside the reservoir) and (b) the abnormally shallow depth of seismicity (resolved from relocations, waveform modeling, and InSAR) compared to previous well‐resolved earthquakes in the Fars arc. The latter is illustrated by Figure 7, comparing the Khalili focal and centroid depths with previous well‐studied earthquakes in the Fars arc.

Conclusions

We present a detailed analysis of the 2019–2020 Khalili seismic sequence in the Fars arc of Zagros Simply Folded Belt. The Zagros is a region of elevated natural seismicity but the proximity of this sequence to the Shanul gas field, from which a very large volume of gas has been extracted over ∼14 years, raised the possibility that these were induced events. We analyzed the sequence using local, regional, and teleseismic data and further constrained the largest earthquake with InSAR modeling. A comparison with previous background seismicity highlights the anomalously shallow depths of the 2019–2020 sequence, suggesting that human‐induced stress changes related to the Shanul gas field may have caused the Khalili seismic sequence. This inference is further supported by the application of a variety of qualitative indicators, but a more sophisticated, probabilistic assessment would require injection/extraction data, which are lacking. This is, to our knowledge, the first possible case of anthropogenic seismicity directly linked to gas field operations in the Zagros. Triggering of the M w 5.4 and 5.7 events highlights the need to identify large faults in the vicinity of active gas fields and to avoid pore pressure perturbations that could destabilize these faults. Iran already hosts a significant number of proven hydrocarbon reservoirs and has huge potential for new discoveries. For instance, in 2019, the NIOC discovered a gas reservoir ∼20 km north of the Shanul field, named Eram, with a very considerable gas capacity. Our results suggest that the exploitation of these reservoirs should be preceded by risk assessment studies and accompanied by the implementation of dedicated, sophisticated monitoring, which would allow seismic activity to be detected early and its evolution tracked.

Code Availability

Some of the maps were prepared using the Pyrocko toolbox (Heimann et al., 2017; https://pyrocko.org/) and GMT 5 software (Wessel et al., 2013; https://www.generic-mapping-tools.org/). For relocation, we used the mloc program (Bondár et al., 2008; https://seismo.com/mloc/). The probabilistic source inversion was performed with the Grond framework (Heimann et al., 2018; https://pyrocko.org/grond/docs/current/). InSAR data were processed using GAMMA software (Wegmüller & Werner, 1997; https://www.gamma-rs.ch/) and downsampled and inverted using codes developed by the Centre for the Observation and Modelling of Earthquakes, Volcanoes, and Tectonics (COMET) group (https://comet.nerc.ac.uk/), available from E. Nissen upon request. The mainshock focal depth was calculated using “Array Beam Depth Tool” tools (Kriegerowski, 2021; https://github.com/HerrMuellerluedenscheid/ArrayBeamDepthTool). The magnitude of completeness was calculated using the ZMAP package (http://www.seismo.ethz.ch/en/research-and-teaching/products-software/software/ZMAP) (Wiemer, 2001).

Supporting Information S1

Click here for additional data file.