Fault mechanics and post-seismic deformation at Bam, SE Iran

Summary

The extent to which aseismic deformation relaxes co-seismic stress changes on a fault zone is fundamental to assessing the future seismic hazard following any earthquake, and in understanding the mechanical behaviour of faults. Here we use models of stress-driven afterslip and viscoelastic relaxation, in conjunction with post-seismic InSAR measurements, to show that there has been minimal release of co-seismic stress changes through post-seismic deformation following the 2003 M_w 6.6 Bam earthquake. Our analysis indicates the faults at Bam remain predominantly locked, suggesting that the co- plus interseismically accumulated elastic strain stored downdip of the 2003 rupture patch may be released in a future M_w 6 earthquake. Our observations and models also provide an opportunity to probe the growth of topography at Bam. We find that, for our modelled afterslip distribution to be consistent with forming the sharp step in the local topography over repeated earthquake cycles, and also to be consistent with the geodetic observations, requires either (1) far-field tectonic loading equivalent to a 2–10 MPa deviatoric stress acting across the fault system, which suggests it supports stresses 60–100 times less than classical views of static fault strength, or (2) that the fault surface has some form of mechanical anisotropy, potentially related to corrugations on the fault plane, that controls the sense of slip.

Seismic cycle, Fault zone rheology, Dynamics and mechanics of faulting

1 INTRODUCTION

The 2003 December 26 Bam earthquake (M_w 6.6) in Kerman Province, SE Iran, occurred along a previously unknown strike-slip fault and killed at least 31 000 people (Jackson et al. 2006; Berberian 2014). Seismological and geodetic source models proposed that slip in the earthquake was concentrated between 2–8 km depth, with a maximum fault slip of 2.7 m (Fialko et al. 2005; Funning et al. 2005; Stramondo et al. 2005; Motagh et al. 2006; Peyret et al. 2007). However, the pattern of aftershocks extended well beneath the co-seismic slip patch, to nearly 20 km depth (Nakamura et al.2005; Tatar et al.2005), suggesting the bottom ∼10 km of the fault zone remained unruptured, yet is able to generate earthquakes. This observation raised serious concerns regarding whether the bottom half of the fault at Bam could produce another damaging earthquake by releasing accumulated interseismic strain in response to co-seismic loading (Jackson et al. 2006). Post-seismic deformation may relax some of the co-seismic stress changes through fault creep (afterslip) or viscoelastic relaxation. In this paper we study the patterns of post-seismic deformation at Bam to investigate the degree of stress relaxation following the earthquake.

Understanding how much afterslip occurred following the 2003 Bam earthquake is pivotal to assessing seismic hazard and the frictional properties of the faults. Recent studies suggest faults that creep extensively during the interseismic period also experience large amounts of post-seismic afterslip relative to co-seismic slip (Heki et al. 1997; Furuya & Satyabala 2008; Barbot et al. 2009; Thomas et al. 2014). This observation implies that the amount of afterslip is related to the degree of locking of the fault zone around the co-seismic slip patch, and the magnitude of the co-seismic stress changes. A previous study of Bam revealed that the post-seismic surface deformation between 2004 and 2007 is consistent with limited afterslip following the 2003 earthquake (Fielding et al. 2009), indicating the fault zone surrounding the co-seismic rupture may remain predominantly locked. We investigate the suggestion of fault locking at Bam further by making additional measurements of the post-seismic ground deformation and performing new modelling of afterslip and viscoelastic relaxation.

Following a review of previous work in Section 2, we develop a new interferometric synthetic aperture radar (InSAR) time-series of ground deformation between 2004 and 2009 at Bam in Section 3, with the intention of investigating the rheology of the fault zone. In Section 4, we present forward calculations of post-seismic afterslip and viscoelastic deformation that are consistent with the InSAR observations. We quantify the maximum possible relaxation of co-seismic loading on the down-dip extension of the fault, below the co-seismic rupture and above the base of the seismogenic layer (∼10–20 km depth range). Subsequently, in Section 5, we investigate post-seismic afterslip on the shallow fault zone (0–10 km) using stress-driven afterslip calculations (e.g. Barbot et al.2009) in conjunction with the post-seismic InSAR measurements.

Models of stress-driven afterslip can also be used to investigate the stress state within the shallow crust. Four kilometres east of Bam is a short-wavelength topographic ridge formed by motion along an underlying reverse fault (Jackson et al. 2006). In Section 6 of this study, we use models of stress-driven afterslip to calculate the pattern of motion along the reverse fault under the ridge, and so test under what conditions co-seismic stress changes can explain both the short-term InSAR observations and long-term growth of the topography along the ridge.

2 REGIONAL AND LOCAL TECTONIC FRAMEWORK

2.1 Regional tectonics

Bam lies on the south-western side of the Dasht-e-Lut (Fig. 1), a desert region devoid of seismicity in its interior. The regional tectonic framework is dominated by the transition from ocean-continent subduction east of 58 °E to continental collision between Arabia and Eurasia in central and western Iran. Convergence between Arabia and Eurasia occurs at 20–30 mm yr⁻¹ (Vernant et al. 2004; Reilinger et al. 2006), with 10–15 mm yr⁻¹ of this motion accommodated by shortening within the Zagros Mountains (Jackson & McKenzie 1984; Vernant et al. 2004). In eastern Iran the majority of the convergence is accommodated by slip along the subduction interface beneath the Makran Ranges (Vernant et al. 2004) . As a result, there is differential motion between eastern and central Iran that is manifest as ∼15 mm yr⁻¹ right-lateral shear on N–S oriented planes (Fig. 1b).

Figure 1.

Overview of the active tectonics of SE Iran. In (a) the location of the maps in (b) and (c) are shown on the country-wide map. (b) depicts the GPS velocities of Walpersdorf et al. (2014) relative to stable Eurasia, omitting the error ellipses (∼1–2 mm yr⁻¹) for clarity, and highlights the right-lateral shear on N–S planes in eastern Iran. (c) shows the distribution of hypocentres from the EHB catalogue (Engdahl et al. 2006) highlighting the concentration of events around the western edge of the Dasht-e-Lut. Also plotted are the locations of major active faults from Walker & Jackson (2004). Focal mechanisms of body-waveform modelled earthquakes along the Gowk Fault to the north of Bam are shown in blue (Berberian 2014). (d) Source models for the 2003 Bam earthquake from a number of geodetic (red), seismological (blue) and joint (green) inversions are shown as lower hemispherical projections for the Arg-e-Bam Fault (AEBF) and Bam-Baravat Fault (BBF), which are the two identified active faults at Bam (see Fig. 2).

Seismicity and fault mapping delineate a number of N–S trending strike-slip faults that accommodate the right-lateral shear around the eastern and western edges of the Dasht-e-Lut (Fig. 1c). The major faults on the west of the Dasht-e-Lut have an estimated slip rate of 2–5 mm yr⁻¹ (Walker et al. 2011; Walpersdorf et al. 2014) and have experienced extensive seismicity over the instrumental period (Berberian 2014), including four M_w > 5 earthquakes along the Gowk Fault Zone alone in the past 40 yr (Berberian 2005) (see Fig. 1c). The deformation around Bam is consistent with this regional strain field, with GPS studies indicating 1–3 mm yr⁻¹ of differential N–S motion (Walpersdorf et al. 2014; Fig. 1b, Supporting Information Fig. S1).

Interseismic measurements of ground deformation using InSAR between 1992 and 1999 revealed no resolvable displacements projected into the satellite line-of-sight (LOS; Fialko et al. 2005). This observation requires surface creep rates on the faults at Bam to be less than the InSAR detection threshold (≲2 mm yr⁻¹ of LOS motion), which is consistent with the slow interseismic loading inferred from GPS.

The lack of observed interseismic LOS displacement across the Bam-Baravat Ridge does not rule out creep on the underlying fault. North of Bam, in the Shahdad fold-thrust belt (Walker & Jackson 2002; Fielding et al. 2004), ∼20 km long, west-dipping reverse faults associated with more developed anticlinal ridges than at Bam are creeping at ∼2 mm yr⁻¹ in LOS (Copley & Jolivet 2016), which is close to the InSAR detection threshold. Therefore it is possible the Bam-Baravat Fault is creeping, but generates surface displacements smaller than the noise in the InSAR measurements over our observation period.

2.2 Local geomorphology and faulting

Prior to 2003 the only active fault identified at Bam was inferred from the presence of a 15 km long, N–S trending escarpment that lies 15–30 m above the surrounding desert plains, known as the Bam-Baravat Ridge (Berberian 1976; Figs 2a and b). The topographic asymmetry of the escarpment, presence of an uplifted water table, incised surface drainage in its western limb, sharp topographic step (Fig. 2b) and folded bedding (Hessami et al. 2005) suggest the ridge was formed by slip on a west-dipping reverse fault (Jackson et al. 2006). Deflection of contemporary drainage channels around its southern tip are indicative of the active southward propagating growth of the ridge (Jackson et al. 2006), but the lack of offset rivers crossing the front of the ridge, or offset qanat tunnels (underground irrigation channels that can be up to 2000 yr old (Jackson 2006)), suggest there is no resolvable long-term strike-slip motion.

Figure 2.

Geomorphology and topography in the Bam region, with the location of co-seismic surface ruptures and aftershocks. (a) is a Landsat 8 image (bands 4, 3, 2) of Bam and the surrounding region. Black dots show the location of co-seismic surface ruptures from the 2003 earthquake (from Jackson et al.2006) and the white dots represent aftershocks recorded by Tatar et al. (2005), which have been located with the HYPO71 algorithm (Lee & Lahr 1972). The surface ruptures form two distinct sets. Those in the west show an en-echelon left-stepping pattern, consistent with right-lateral motion, and lie above the Arg-e-Bam Fault (AEBF), but are not associated with any overlying topography. The fractures in the east occur along the front of the Bam-Baravat Ridge, and are associated with the surface projection of the Bam-Baravat Fault (BBF). (b) is an ALOS 30 m DEM highlighting the location and morphology of the Bam-Baravat Ridge. (c) shows the variation in along-strike elevation of the ridge’s crest relative to the adjacent plains to the east. (d) shows a topographic profile perpendicular to the ridge, with the regional trend removed. The patterns of relative uplift are consistent with formation due to slip along a west dipping reverse fault.

The length scale of the topography perpendicular to the Bam-Baravat Ridge is 3–4 km (Fig. 2d). Growth of the ridge is therefore sensitive to slip along the top ∼5 km of the underlying fault, as any fault motion at depths ≳5 km would lead to a broader pattern of uplift. The sharp step in elevation across the escarpment requires that dip-slip fault motion extends to near the surface. Similar structures have been recognized elsewhere in Iran such as at Tabas-e-Golshan (Walker et al. 2003), Ferdows (Walker et al. 2003), Sefidabeh (Berberian et al. 2000), Rigan (Walker et al. 2013) and Shahdad (Berberian et al. 2001).

2.3 Co-seismic observations from the 2003 Bam earthquake

Co-seismic deformation provided numerous opportunities to probe the fault geometry and slip at Bam, with studies using a combination of teleseismic body waveform modelling (Talebian 2004; Jackson et al. 2006; Poiata et al. 2012), InSAR (Talebian 2004; Fialko et al. 2005; Funning et al. 2005; Stramondo et al. 2005), optical-image pixel tracking (Fu et al. 2004; Binet & Bollinger 2005; Peyret et al. 2007; de Michele et al. 2008) and levelling (Motagh et al. 2006). These studies identified an N–S trending right-lateral fault that dips ∼85° to the east beneath a western set of surface fractures (Fig. 2a), with up to 2.7 m of slip concentrated between ∼2 and 8 km depth. There is also evidence for oblique co-seismic slip on a steeply dipping (∼60°) reverse fault beneath the Bam-Baravat Ridge, which is consistent with the forming a second set of surface fractures (Fig. 2a) and leads to a significant improvement in the fit of the co-seismic source models to the InSAR observations (Jackson et al. 2006). As a result, in this paper we use a two fault co-seismic slip model, with the geometry and slip distribution of Funning et al. (2005) used throughout. We refer to the strike-slip fault associated with the western set of surface fractures as the Arg-e-Bam Fault (AEBF), and the reverse fault that projects to the surface at the front of the Bam-Baravat Ridge as the Bam-Baravat Fault (BBF; see Fig. 2a).

3 OBSERVATIONS: POST-SEISMIC InSAR

An earlier study by Fielding et al. (2009) produced an InSAR time-series of ground deformation between 2004 and 2007 at Bam, to study the post-seismic compaction of the shallow fault zone directly over the surface projection of co-seismic rupture on the Arg-e-Bam Fault. However, as a result of their processing technique, the final time-series contained extensive decorrelation over the urban and vegetated regions of Baravat and Bam. In this study our focus is on the ground displacements associated with afterslip around the ends of Arg-e-Bam Fault and across the Bam-Baravat Fault, therefore we require improved coherence in the city of Bam in order to resolve ground deformation that is not present in the study of Fielding et al. (2009). To address this issue we produce a new, extended time-series of post-seismic deformation using SAR data collected by Envisat between 2004 and 2009, from ascending track 156 and descending track 120.

The SAR data were processed using a combined persistent-scatter and small-baselines approach implemented using StaMPS (Hooper et al. 2004, 2007). Combined processing optimizes the coherent pixel density in the radar scene, which is important in this study because the vegetated regions of Bam and Baravat will contain some pixels with stable phase characteristics through time (e.g. buildings, walls), but these might be surrounded by pixels that decorrelate rapidly, therefore would not be selected by small-baselines processing alone. Single-look complex images were focused using ROI_PAC (Rosen et al. 2004) and interferograms generated using DORIS (Kampes & Usai 1999), with Delft orbital data and topographic corrections performed using an SRTM 3 arc-second elevation model (Farr et al. 2007). Interferograms were then unwrapped using a statistical-cost network-flow algorithm (Chen & Zebker 2001). We attempted to remove any relic orbital effects and a component of long-wavelength atmospheric noise by estimating and subtracting a linear ramp from each interferogram. The resulting data set consists of an ascending-track time series that includes 25 SAR acquisitions and a descending-track time series that includes 28 acquisitions (Figs 3 and 4).

Figure 3.

Baseline-time plot of the interferogram pairs used in creating the small-baselines and persistent-scatter time series, and the chain-stack interferograms. The grey circles correspond to individual SAR acquisitions, with the black lines between each SAR acquisition representing a small-baselines interferogram pair. Grey circles without any interferograms represent SAR acquisitions that have a high component of atmospheric noise. The red circle is the single-master used in the persistent-scatter processing. The black dashed line is the time of the 2003 Bam earthquake. The red lines represent interefrogram pairs used in creating the chain stacks shown in Fig. 8. We use spatial baseline interferograms longer than is common because the stable and flat desert region around Bam leads to minimal temporal and geometrical decorrelation.

Figure 4.

Time-series of post-seismic LOS displacements derived from the combined persistent-scatter and small-baselines processing (PS+SB) at Bam from January 2004 to 2009. (a) shows the cumulative surface displacements at different times, with the time of each snapshot after the first acquisition in January 2004 shown in the top left corner. Black dots are the co-seismic surface ruptures, and the satellite look vector projected into the horizontal plane is shown in the bottom right of each image (the incidence angle is 23°). We use the convention that positive LOS change is motion towards the satellite. (b) Profiles of LOS displacement extracted across profile A–A^΄ from both the combined PS+SB time-series after 5.7 yr, and a set of chain-stacked interferograms covering the same time period (red lines in Fig. 3). The location of the profile’s intersection with the two faults at the surface are outlined as black vertical dashed lines where AEBF is the Arg-e-Bam Fault and BBF is the Bam-Baravat Fault. These profiles highlight the step change in LOS across the Bam-Baravat Fault, which we interpret to be a tectonic signal reflecting shallow post-seismic creep. The sharp negative LOS change over the Arg-e-Bam Fault was studied by Fielding et al. (2009) and has been interpreted to reflect compaction of the shallow fault zone following co-seismic damage.

In addition, we produced 94 descending and 54 ascending track interferograms with perpendicular baselines <500 m, processed with ROI_PAC (Rosen et al. 2004), using Delft orbits and the SRTM 3 arc-second elevation model to correct for the effect of topography on radar phase. The interferograms were unwrapped using the branch-cut algorithm (Goldstein & Werner 1998) and a linear ramp was removed. In order to map the full 5.7 yr of post-seismic deformation without using heavily decorrelated and long temporal baseline interferograms, we formed a series of chain-stacks by summing the LOS motions in multiple interferograms that share a common master or slave image (see red lines, Fig. 3), and have short temporal and perpendicular baselines (e.g. Biggs et al.2007). The chain stacks act as a quality control on the StaMPS time-series output, and provide us with measurements of the long-wavelength ground deformation (>10–20 km) associated with post-seismic processes occurring at depth. Further details of the interferograms used are provided in the Supporting Information. All InSAR measurements use the convention that positive LOS changes correspond to motion of the ground towards the satellite.

3.1 Post-seismic InSAR: results

The time-series of post-seismic ground deformation is shown in Fig. 4. Fielding et al. (2009) concentrated on the short-wavelength signals associated with compaction of the top 1–2 km of the fault zone in response to co-seismic dilation. We find similar patterns of LOS displacement to Fielding et al. (2009), and use our longer time-series to confirm that both the shallow fault zone compaction and afterslip signals share a similar temporal evolution.

In this study our focus is on the lobes of positive and negative LOS motion radiating from the tips of the Arg-e-Bam Fault and changes in LOS across the Bam-Baravat Fault (Fig. 4). The spatial pattern of the signals around the tips of the Arg-e-Bam Fault, and their temporal evolution (Fig. 5), are consistent with deformation caused by transient post-seismic strike-slip motion on the periphery of the co-seismic rupture (Fielding et al. 2009). The temporal evolution of ground deformation shows initial rapid motion in the first 2 yr of the post-seismic period, after which deformation slows significantly until the visible post-seismic transient is complete within our observation time-scale (Fig. 5).

Figure 5.

Relative LOS displacement evolution for the four boxes highlighted in Fig. 4. Each point represents the difference in the mean LOS displacement within the box highlighted in Fig. 4 and a collection of reference pixels at a point far from the active faults (point 5 in Fig. 10a). The error bars represent the combined standard deviation of pixel LOS in each box and the reference pixels. The solid black lines are the best-fitting functions from a 1-D spring-slider analogue model of the form: u(t) = at + bt_rlog [1 + d(exp (t/t_r) − 1)] (Perfettini 2004). The red line is the best-fitting stress-driven afterslip model when the whole fault surface around the compacted co-seismic rupture is able to creep (surface displacements shown in Fig. 10). The blue line is a stress-driven afterslip model in which only small patches around the edge of the compacted co-seismic slip patch are able to creep, with the remainder of the fault remaining locked. Both the fully creeping and locked models are driven by the same compacted co-seismic slip distribution. The surface displacements associated with the locking model are shown in Supporting Information Fig. S3(b).

Across the front of the Bam-Baravat Ridge a sharp, negative ∼1 cm step change in LOS is visible from west to east in both the descending and ascending interferograms (see Fig. 4, profile A–A^΄ at BBF), which increases in amplitude consistently through time. This signal is too short wavelength and temporally too consistent to be associated with atmospheric noise. The positive LOS displacement of the ridge relative to the adjacent plain is opposite to what would result from aquifer discharge through the local qanat tunnels at Baravat, as this would cause subsidence of the ridge and therefore a negative LOS change relative to the adjacent plains. The change in elevation across the ridge front (15–30 m) and DEM errors are too small to account for the large step in LOS change. We therefore conclude that this signal is tectonic and represents shallow, thrust-oriented creep on the Bam-Baravat Fault underlying the ridge. As there is no resolvable interseismic LOS displacement across the ridge front between 1992 and 1999 (Fialko et al. 2005), shallow creep on the Bam-Baravat Fault probably accelerated following the 2003 Bam earthquake in response to co-seismic stress changes.

We find no obvious signals associated with tectonic motions that occur on length-scales of 10–20 km, suggesting any viscoelastic deformation or fault creep at depth must lead to surface displacements smaller than the amplitude of the noise in the InSAR measurements (∼1.5 cm at 10–20 km wavelength). In the following section we use these observations, in conjunction with evidence from aftershocks, slip inversions, and forward modelling, to investigate the rheology of the fault zone at Bam, and quantify the contribution of post-seismic deformation mechanisms to the relaxation of co-seismic loading.

4 DEEP POST-SEISMIC DEFORMATION

4.1 Limited post-seismic moment release at Bam

Fielding et al. (2009) used their post-seismic InSAR results to perform a slip inversion that revealed there has been only 10–15 cm of shallow afterslip above 5 km depth on the Arg-e-Bam Fault, and 10–15 cm of afterslip localized to a 4 × 3 km patch along the base of the Bam-Baravat Fault. This slip model accounts for a post-seismic moment release of 1–3 × 10¹⁷ Nm, corresponding to 1–3 per cent of the co-seismic geodetic moment (Funning et al. 2005). Such a low post- to co-seismic moment release ratio (M_ps/M_cs) can be seen to lie outside of an empirically defined pattern produced from a compilation of fault slip inversions using similar data and techniques to that performed at Bam (Fig. 6). This compilation includes large magnitude subduction zone earthquakes that cause significant viscoelastic relaxation, with the associated surface deformation potentially being mapped into afterslip-only inversions, leading to an over-estimate of M_ps/M_cs. However, this effect is likely to be limited for smaller magnitude, shallow earthquakes, which induce far less viscoelastic relaxation. As the difference in M_ps/M_cs between Bam and other earthquakes is consistent across both subduction zone events and smaller intra-continental events, the inference that M_ps/M_cs is anomalously low for Bam appears robust.

Figure 6.

Compilation of co-seismic and post-seismic geodetically derived moment release from published studies (listed in the Supporting Information Table S1). Uncertainties in the measurements of moment release are estimated from the spread in values in a range of inversions, therefore points without error bars are earthquakes with only one published inversion that does not quote uncertainties. The moment release ratio at Bam lies outside of the empirical scaling shown by the majority of other earthquakes, where the scaling between co-seismic and post-seismic moment release is generally 1–0.1.

The estimated post-seismic moment release in the shallow fault zone at Bam is likely to be robust because: (1) the time-series in Fig. 5 shows the low M_ps/M_cs is not due to incomplete sampling of a long-lived post-seismic transient, and (2) extrapolating the time-series of LOS displacement in Fig. 5 back to the time of the earthquake using a logarithmic function suggests that only minor (0.5–2 cm) LOS displacement is likely to have been missed in the early post-seismic period, where there is no SAR data. This unsampled deformation does not account for the order of magnitude difference between M_ps/M_cs observed at Bam and the majority of other studied earthquakes (Fig. 6).

A low M_ps/M_cs, in conjunction with the observation of limited interseismic creep (Fialko et al. 2005), indicates the majority of slip on the faults at Bam occurs in earthquakes or through interseismic creep too slow to observe with present methods. Fault zones that show the opposite behaviour, in which M_ps/M_cs is ≳80 per cent, are consistently located near parts of the fault surface that have interseismic creep rates similar to the far-field loading rate (e.g. the Parkfield section of the San Andreas Fault (Barbot et al. 2009) and the Longitudinal Valley Fault in Taiwan (Thomas et al. 2014)). For the high M_ps/M_cs events there is limited locking around the co-seismic rupture, therefore post-seismic creep is not limited by the elastic resistance caused by adjacent locked zones.

One potential source of unaccounted post-seismic moment release at Bam is afterslip between 10 and 20 km depth that remains invisible in the InSAR observations due to the effects of long wavelength atmospheric noise and the low sensitivity of surface deformation measurements to slip at depth. Deep afterslip has primarily been found to be co-located with, and share the same temporal evolution as, deep aftershock sequences (Perfettini 2004; Savage et al. 2005; Hsu et al. 2006; Peng & Zhao 2009; D’Agostino et al. 2012; Mencin et al. 2016). In addition, aftershock sequences are often localized around the edges of the co-seismic rupture, and share at least one nodal plane with the mainshock (Bodin & Horton 2004; Tatar et al. 2005), indicating that some aftershocks represent failure on the same fault surface as the mainshock. Therefore local aftershock studies may provide some insight into whether deep afterslip is occurring at Bam.

4.2 Aftershocks at Bam

The aftershock sequence recorded at Bam consists of catalogues produced by a dense local array of 23 stations over January 2004 (Tatar et al. 2005) and a smaller local array of nine stations recording events between February and March 2004 (Nakamura et al. 2005). We are interested in the depth extent of aftershocks, therefore only consider the catalogue of Tatar et al. (2005), as the 5 km station spacing used in their study is small enough to reliably locate events within the required depth range with around ±2 km uncertainty. The aftershock distribution indicates that seismicity was focused on the down-dip extension of the co-seismic rupture, between 10 and 18 km depth, and almost absent from the region directly over, and adjacent to, the co-seismic slip patch. A similar pattern was seen following the 2011 Christchurch earthquake (Elliott et al. 2012) (http://quakesearch.geonet.org.nz/), and the 2004–2009 Qeshm Island earthquakes in Iran (Nissen et al. 2010). The aftershocks follow a decay in frequency as a function of time, but show no particular temporal evolution in their locations (Figs 7b and c). Finally, the majority of the well-constrained aftershock mechanisms have at least one nodal plane compatible with right-lateral motion along N-S planes, the same mechanism as the main shock (Tatar et al. 2005). All three observations are consistent with transient afterslip loading and breaking small asperties along the downdip edge of the co-seismic slip patch. An absence of aftershocks adjacent to the co-seismic rupture indicates this region is either too strong to break in response to co-seismic stress changes, or that the region is creeping and contains no locked patches.

Figure 7.

Characteristics of the aftershocks recorded over January 2004 by Tatar et al. (2005) and regional seismicity recorded by the IIEES network in Iran. (a) shows the hypocentral location of each aftershock projected onto the fault plane, colour coded by the time after the main shock, with contours of the co-seismic slip distribution of Funning et al. (2005). (b) highlights that there is no temporal evolution in the location of aftershocks. (c) shows the decay in aftershock frequency as a function of time following the main shock from Tatar et al. (2005) and (d) shows the decay in regional seismicity with time following the Bam earthquake in the region between 58.2°–58.5° longitude and 28.9°–29.2° latitude from the catalogue of Karimiparidari et al. (2013) (magnitude of completeness ∼4.5). There was no recorded seismicity in this area before December 2003.

4.3 Models of deep afterslip

Post-seismic InSAR measurements reveal no evidence for a long-wavelength deformation signal that might be associated with afterslip beneath the co-seismic rupture on the Arg-e-Bam Fault (Figs 8a and c). This suggests that any afterslip occurring below the co-seismic rupture cannot lead to surface displacements greater than the noise levels in the interferograms, which provides an upper bound on the total amount of afterslip at depth.

Figure 8.

A comparison of long-wavelength LOS displacements in chain-stack interferograms with calculations for viscoelastic relaxation and deep afterslip. (a) and (d) show a descending and ascending chain-stack interferograms over the period displayed in the title with the format YYYYMMDD. The dashed box in (a) is the area shown in Fig. 4, and the thin lines represent the surface projections of the faults from the slip model of Funning et al. (2005). (b) and (e) are results of the LOS displacements from a forward calculation in which the ground surface deforms in response 5.7 yr of viscous relaxation beneath a 15 km thick elastic lid with an underlying half-space of Newtonian viscosity 1 × 10¹⁹ Pa s. This is considered the weakest possible half-space viscosity that produces no viscoelastic deformation visible in the interferograms. (c) and (f) are the LOS displacements following 50 cm of afterslip along the base of the Arg-e-Bam Fault between 10 and 18 km depth, which we suggest is the maximum amount of afterslip that could occur without the signal becoming visible in the chain stacks (a,d).

To quantify this upper bound we perform a grid search of dislocation calculations in an elastic half-space (Okada 1985) to find the ground deformation resulting from afterslip in the region of aftershock activity, on the downward extension of the Arg-e-Bam Fault. The elastic properties of the crust are taken from the seismic velocity model of Tatar et al. (2005), which indicates the shear modulus of the upper crust at Bam is 30 ± 10 GPa, and we assume a Poisson ratio of 0.25. The along-strike length, dip and rake of the fault are fixed to those derived from the overlying co-seismic slip patch in the model of Funning et al. (2005) and we vary the amount of afterslip and down-dip width of the afterslip patch. Surface displacements for each calculation are converted into satellite LOS and are used to compute the signal-to-noise ratio (SNR) by finding the ratio between the maximum modelled LOS changes and the noise in the interferograms. An SNR > 1 means that afterslip-related ground deformation is visible in the interferograms. The contribution of noise to apparent LOS variations in the interferograms is estimated to be ∼1.5 cm for signals with a wavelength between 10 and 20 km, based on the mean amplitude of LOS variations in regions of the interferograms that are considered to be undeforming (see Figs 8a and c).

Assuming the region experiencing afterslip overlaps with the area of aftershock activity on the lower half of the Arg-e-Bam Fault (between 10 and 18 km), we find there must be <50 cm of deep afterslip for the resulting surface displacements to remain invisible in our InSAR measurements (Fig. 9a). This corresponds to a maximum ‘invisible’ geodetic moment release of 1.5 × 10¹⁸ Nm, which is 10 per cent of the co-seismic moment release. Hence, it is possible that afterslip at depth could account for the difference in M_ps/M_cs between Bam and the lower end of other studied earthquakes’ moment scaling (Fig. 6).

Figure 9.

Calculations of ground deformation for different amounts of afterslip and viscoelastic relaxation along the base of the Arg-e-Bam Fault, converted into a signal-to-noise ratio (SNR). An SNR = 1 contour marks the upper bound on the amount of deep afterslip and visco-elastic relaxation that could occur without being visible in the interferograms. (a) shows the SNR for varying amounts of afterslip on the area below the co-seismic slip patch over the post-seismic observation period. The purple circle corresponds to the model shown in Figs 8(c) and (f). (b) shows the SNR for a set of viscoelastic calculations with variable elastic layer thickness and half-space viscosity. The purple circle corresponds to the model shown in Figs 8(b) and (e).

The small amount of both co- and post-seismic fault slip (i.e. < 50 cm) between 10 and 18 km depth relative to the shallow co-seismic slip (2.7 m) suggests <20 per cent of the co-seismic shear stress changes on the bottom 10 km of the fault zone have been released through post-seismic fault creep. An alternative mechanism for stress release beneath the fault zone is viscous deformation of the lower crust (Freed & Bürgmann 2004; Ryder et al. 2007). In the following section we explore the extent to which viscoelastic relaxation can release the co-seismic stress changes imposed on the ductile part of the lithosphere, and determine the possible rheology of the region beneath the brittle fault zone at Bam.

4.4 Forward calculations of viscoelastic relaxation

We use a simple Earth model of an elastic layer overlying a uniform Newtonian viscous half-space that is perturbed due to co-seismic stress changes from the Bam earthquake. Calculations for the resulting deformation, and its effect on surface displacements, are performed using RELAX (Barbot et al. 2009; Barbot & Fialko 2010a,b) for various thicknesses of the elastic layer and viscosities of the underlying half-space. The surface displacements after 5.7 yr of deformation are projected into the satellite LOS and compared to the noise levels in the interferograms.

Based on a realistic estimate of the elastic layer thickness of 15 km, which is smaller than the seismogenic thickness of 18–20 km (Tatar et al. 2005), we find that the minimum possible viscosity of the underlying half-space is ∼1 × 10¹⁹ Pa s (Fig. 9b). If the viscosity were any lower, the surface deformation resulting from viscoelastic relaxation would become visible in our interferograms, which is inconsistent with the observations. Any viscoelastic relaxation that occurs will begin to reload the overlying fault zone. We find that after 5.7 yr viscoelastic relaxation may increase the average shear stress on the bottom 10–15 km of the Arg-e-Bam Fault by a maximum of 15 per cent of the co-seismic shear stress change on the same fault patch.

The calculations presented above are sensitive to the assumed thickness of the elastic layer (Fig. 9b), but nonetheless suggest the lack of observed long-wavelength post-seismic deformation reflects limited deep afterslip and limited viscoelastic relaxation beneath the fault that ruptured in the 2003 earthquake. We discuss this finding below in the context of seismic hazard.

The InSAR observations reveal short-wavelength ground deformation is focused around the tips of the co-seismic surface fractures, indicating there has been deformation around the edges of the Arg-e-Bam Fault in the shallow crust (<5 km). In the following section we assess whether this deformation is consistent with models of afterslip on a uniformly creeping fault in response to co-seismic stress changes, and probe the frictional properties of the fault zone.

5 SHALLOW POST-SEISMIC DEFORMATION

To relate the co-seismic stress changes on the fault zone following the Bam earthquake to the resulting afterslip inferred from the InSAR measurements, we employ the rate-and-state friction parametrization (Dieterich 1978; Ruina 1983; Scholz 1998). This empirical description of friction has been applied to both laboratory (Ruina 1983; Marone 1998) and natural faults (Perfettini & Avouac 2007; Barbot et al. 2009; Copley & Jolivet 2016), and appears to be a good approximation of the constitutive relationship that governs fault mechanics in the upper crust. A preliminary test of the application of this formulation to the faults at Bam is to compare the observed deformation with simple analogue models of rate-strengthening fault creep.

5.1 1-D fault analogue models

The temporal evolution of the LOS displacements in Fig. 4 can be compared to a simple 1-D fault analogue (Marone et al. 1991; Perfettini 2004; Barbot et al. 2009). One such model treats the fault mechanics as a spring-slider block system with friction along its base, which is governed by a simplified rate-strengthening friction law (Perfettini 2004). According to this model, ground deformation following a stress change (Δτ) will evolve as: u(t) = αV₀t + βV₀t_rlog [1 + d(exp (t/t_r) − 1)], where α and β are geometric constants, V₀ is the interseismic creep rate, t_r is the relaxation time and d = exp (Δτ/aσ), where a = ∂μ/∂log (V) gives the dependence of friction on sliding velocity, and σ is the effective normal stress acting on the block.

We find this model provides a good fit to the observations with the same relaxation time (t_r = 1.3 yr) and friction parameters for all the LOS time-series extracted over areas of afterslip (Fig. 5, black line), but the solutions are non-unique due to trade-offs between the geometric constants and friction parameters. The consistent relaxation time for all time-series suggests that the post-seismic deformation at Bam is compatible with creep on patches of the fault with relatively uniform friction properties. The following modelling builds on this suggestion by implementing a more geometrically realistic quasi-static model of afterslip driven by co-seismic stress changes.

5.2 Stress-driven afterslip models

Measurements of post-seismic deformation often record the effects of multiple superimposed deformation mechanisms including afterslip, viscoelastic relaxation and poroelastic rebound (Jónsson et al. 2003; Hsu et al. 2006; Jónsson 2008). We have shown that viscoelastic relaxation and deep afterslip have no resolvable effect on the short-wavelength surface deformation at Bam, hence can be ignored when modelling the ground deformation. We have also calculated the deformation resulting from poroelastic rebound by taking the difference between the co-seismic deformation for an elastic model calculated with a drained and undrained Poisson’s ratio (ν_d = 0.25 and ν_ud = 0.28, respectively) (Peltzer et al. 1998). Deformation caused by poroelastic relaxation is found to correspond to ∼1 cm LOS motion (Supporting Information Fig. S2), which is small compared to the observed transient LOS motion (∼3–4 cm). As a result, we can assume the InSAR observations around the fault tips, on a wavelength of 2–10 km, are dominated by shallow afterslip.

To model the surface deformation caused by stress-driven afterslip at Bam we perform a series of calculations using the quasi-static fault mechanics code RELAX (Barbot et al. 2009; Barbot & Fialko 2010a,b; Rousset et al. 2012), and compare our forward calculations to InSAR observations of ground deformation to determine the frictional properties of the fault zone.

The initial condition in the models is defined as the geometry and co-seismic slip distribution along the two faults taken from the inversion results of Funning et al. (2005), which are used to calculate the stress changes due to co-seismic slip. The magnitude of the stress changes are a function of the input slip model and the elastic moduli of the half-space (taken from the velocity models of Tatar et al.2005). Following the approach of Barbot et al. (2009) we adapt the model of Funning et al. (2005) by compacting the slip distribution to account for three effects: (1) the dynamic propagation of co-seismic slip into regions that otherwise would creep (Noda & Lapusta 2013), (2) the smoothing of the slip distribution performed in the co-seismic slip inversions, and (3) the possibility of afterslip overlapping with regions of reduced co-seismic slip (e.g. Jacobs et al.2002; Copley et al.2012). In addition, near-surface co-seismic slip under Bam city is removed, as this occurs beneath a region that is decorrelated in the co-seismic interferograms (Funning et al. 2005), and therefore is not robust. The slip distribution is compacted by removing areas in which slip is less than some factor (ζ) of the peak co-seismic slip, with the remaining slip distribution re-scaled to maintain the same co-seismic moment.

In our models the temporal evolution of afterslip, driven by the co-seismic stress changes, is governed by a generalized rate-strengthening friction relationship: V = V₀sinh (Δτ/aσ), where V is the sliding velocity, Δτ is the stress change, and V₀ and aσ are constitutive parameters that are varied to fit the observed time-series of ground deformation (Barbot et al. 2009). As the two faults we model at Bam are in close proximity and cut the same geology, both are assumed to have the same frictional properties.

Co-seismic stress changes alone are assumed to be responsible for driving afterslip, with the contribution of any long-term creep assumed to be small (as seen in Fialko et al.2005). Once the pattern of afterslip has been calculated, the resulting surface deformation is computed using Fourier domain Green’s functions (Barbot & Fialko 2010a), allowing comparison of the model with the InSAR time-series.

To find the best-fitting afterslip model and friction parameters we performed a grid search of forward calculations varying aσ, V₀ and ζ to find the solutions that minimize the reduced chi-square misfit between the InSAR time-series measurements and the corresponding stress-driven afterslip model predictions. We compute the misfit as:

\begin{equation} \chi _{\text{red}}^2 = \frac{1}{N-P} \sum _{i=1}^{n} \sum _{j=1}^k \left[ \frac{(\phi _{ij}^{\text{mod}} - \phi _{ij}^{\text{obs}})}{\sigma _{i}} \right]^{2}, \end{equation}

(1)

where N is the number of independent measurements, P is the number of variables used in the models, |$\phi _{ij}^{\text{mod}}$| is the modelled LOS displacement in the j’th pixel of the i’th interferogram, |$\phi _{ij}^{\text{obs}}$| is the same but for the observed interferograms, and σ_i is the estimate of the atmospheric noise standard deviation calculated in a non-deforming region of the ith interferogram. A |$\chi ^2_{\text{red}} =1$| indicates that the model fits the data to within the uncertainties of the observations, but without over-fitting.

5.3 Results of afterslip modelling

The best-fitting afterslip model is shown in Fig. 10 and corresponds to a |$\chi ^2_{\text{red}} = 1.17$| with the associated fitting parameters being aσ = 5.5 MPa, V₀ = 5 mm yr⁻¹ and ζ = 0.4 (i.e. only fault patches with co-seismic slip >1 m are used in the input slip model), although the models are relatively insensitive to variations in aσ (Fig. 10c). Other studies have found similar values for the fault friction parameters with aσ in the range 0.1–10 MPa (Hsu et al. 2006; Perfettini & Avouac 2007; Barbot et al. 2009; Wei et al. 2015) and V₀ being on the order of the interseismic loading rate, which at Bam is ∼1–3 mm yr⁻¹ (Walpersdorf et al. 2014).

Figure 10.

Results for the best-fitting stress-driven afterslip model on creeping faults. (a) shows snapshots of the observed post-seismic InSAR time-series at the end of the observation period, the model predictions for the best-fit stress-driven afterslip calculation, and the misfit between the two. Small black dots are the surface fractures, and the surface projection of the faults used in the calculations are shown in the modelled interferograms. The numbered dots represent locations in which a time-series of LOS displacements is extracted in Fig. 5, and point 5 is the location of the reference pixels. (b) shows the afterslip distribution that leads the surface displacements for (a). Arrows define the slip vectors showing motion of the hanging wall relative to the footwall and blank areas either represent the co-seismic slip asperity or regions of no afterslip. (c) shows the reduced chi-square between model and InSAR time-series for the grid search of friction parameters, with the best-fitting model in purple. The model shown has a geodetic moment release of 5.6 × 10¹⁷ Nm, 5 per cent of the co-seismic geodetic moment.

The calculated afterslip distribution after 5.7 yr extends outwards from the edge of the compacted co-seismic rupture on both faults and reaches a peak slip of around 20–30 cm (Fig. 10b), which is significantly smaller than the 60–90 cm of afterslip needed to fully relax co-seismic stress changes (see Supporting Information Fig. S3c). On the Arg-e-Bam Fault, slip is almost completely right-lateral and occurs across the depth range of 0–15 km, whilst on the Bam-Baravat Fault the direction of the calculated slip vector is highly variable and afterslip >10 cm is constrained only to the lower half of the fault.

Misfits in the spatial distribution of ground deformation between the models and observations are on the order of the atmospheric noise within the interferograms (Fig. 10a), apart from the narrow zone of negative LOS motion in both ascending and descending tracks directly over the Arg-e-Bam Fault. This negative LOS motion is caused by compaction of the top 1–2 km of the fault zone (Fielding et al. 2009), and is not modelled here as it is not associated with afterslip. This does not affect our estimates of the best-fitting model as we do not allow afterslip directly beneath the region of shallow fault zone contraction.

The temporal evolution of ground displacement is fit less well than the spatial pattern, as seen in the descending track of box 1 and the ascending track of box 3 (boxes shown in Fig. 4, time-series as red lines in Fig. 5). The LOS changes in box 3 are extracted over the eastern edge of Bam city and are both equally negative in the ascending and descending track, suggesting the ground surface is subsiding in this area. It may be that the tectonic signal has been contaminated by some anthropogenic effects (e.g. groundwater extraction). The misfits in box 1 show that the afterslip model predicts a much longer relaxation time than observed, and cannot fit the period of rapid deformation in the first year after the earthquake, which is a feature of the time-series in all areas of post-seismic transient deformation associated with afterslip. The best-fit model also predicts that afterslip will continue long after our observations end, which is not consistent with the InSAR time-series (Fig. 5).

The calculations described above cannot fit the observed post-seismic surface displacements without either (1) over-predicting the amplitude of surface deformation late in the post-seismic period, or (2) underpredicting the amount of early, rapid post-seismic deformation (see red lines, Fig. 5). This feature of the models suggests they overpredict the total amount of stress-driven afterslip following the 2003 earthquake. Calculations that simulate afterslip due to complete relaxation of the co-seismic stress changes confirm this proposition, revealing the resulting surface deformation should be 3–4 times larger than observed (see Supporting Information Figs S3a and c). In the following section we investigate whether this effect could be due to the elastic locking of the fault around the co-seismic rupture at Bam.

5.4 Locking on the faults at Bam

Afterslip on a creeping fault will continue until the stresses that drive creep have been relaxed, or are balanced by resistance to creep. As described in the previous section, calculations that simulate the complete relaxation of co-seismic stress changes at Bam via afterslip surrounding the co-seismic rupture significantly overpredict the surface deformation (Supporting Information Fig. S3c). This suggests either: (1) the calculated co-seismic stress changes that drive fault creep are too high, or (2) that creep is resisted by locked portions of the fault surface.

The co-seismic stress changes that drive afterslip in the models are controlled in part by the input slip distribution, which includes both the co-seismic slip inversion used and the compaction factor (ζ) applied. The geodetic moment of the input slip model would have to be as little as ∼10 per cent of the best-fitting co-seismic moment of Funning et al. (2005) in order for stress-driven afterslip on fully creeping faults to cause surface deformation consistent with the post-seismic InSAR measurements. Multiple independent slip inversions for the 2003 Bam earthquake (Fialko et al. 2005; Stramondo et al. 2005; Motagh et al. 2006) indicate the co-seismic moment can conceivably only be as low as 60 per cent of that suggested by Funning et al. (2005). Therefore uncertainties in the input co-seismic slip model cannot account for the over-prediction of afterslip following the Bam earthquake.

The elastic moduli of the upper crust will scale the co-seismic stress changes driving afterslip by making the region around the co-seismic rupture more or less compliant. For the total amount of surface deformation predicted in the stress-driven afterslip calculations described above to be consistent with the observed post-seismic deformation would require the elastic moduli of the crust at Bam to be only ∼10 per cent of the best-fit values determined from the seismic velocity models of Tatar et al. (2005). By combining the uncertainties from the velocity models with a realistic range in crustal densities (2500–2900 kgm⁻³), we find the elastic moduli are likely no lower than ∼75 per cent of the best-fit values of Tatar et al. (2005). Therefore uncertainties in the bulk elastic properties of the crust cannot account for the over-prediction of afterslip at Bam. Geodetic and seismological studies have also suggested that compliant fault damage zones can have elastic moduli some 50 per cent less than the adjacent undamaged country rock (Fialko et al. 2002; Cochran et al. 2009); even so, this still cannot reduce the elastic moduli enough to reconcile the model predictions with the InSAR observations.

Inelastic yielding of the near-surface during co-seismic rupture and slow interseismic deformation have also been suggested to reduce the stress changes available for driving post-seismic slip (Fialko et al. 2005; Kaneko & Fialko 2011). Although inelastic yielding models can account for some of the limited shallow post-seismic afterslip, it cannot account for the lack of deep (>5 km) post-seismic afterslip down-dip and along-strike from the co-seismic rupture, as the high confining pressures at depth reduce the inelastic response of the crust to the dynamic stress changes in the earthquake (Kaneko & Fialko 2011). Significant interseismic creep in the middle of the seismogenic layer is also unlikely as this would prevent the accumulation of elastic strain at the depths in which earthquakes nucleate. This model is therefore limited to possibly accounting for the lack of afterslip in the shallow fault zone only, although this suggestion remains unproven.

We present an alternative model to explain the limited post-seismic afterslip at Bam that is related to lateral and vertical variations in the frictional properties of the fault surface. Our models suggest that the amplitude and pattern of post-seismic afterslip may be limited by locked areas of the fault. This scenario is consistent with recent observations in which earthquakes that show significant afterslip-related moment release compared to co-seismic moment release (e.g. 2004 Parkfield (Barbot et al. 2009), 2003 Chengkung (Thomas et al. 2014)) are surrounded by creeping fault zones. In these cases, stress-driven afterslip can relax the majority of the co-seismic stress changes on the adjacent fault zone, and large amounts of afterslip occurs. On the contrary, if at Bam the adjacent fault zone remains locked and provided an elastic resistance to creep, it could possibly limit the amount of afterslip. Locking would also be consistent the lack of observed interseismic creep in the pre-earthquake InSAR measurements (Fialko et al. 2005).

To test this proposition we performed calculations in which the creeping area on both faults is restricted to small (<2 km) patches that are located around the edge of the compacted co-seismic rupture, with the rest of the fault zone remaining locked (see Supporting Information Fig. S4 for distribution of creeping patches). This geometry allows for the possibility that afterslip occurs on parts of the fault that also experienced co-seismic slip, similar to the patterns observed following the 1999 Hector Mine and 2006 Mozambique earthquakes (Jacobs 2002; Simons et al. 2002; Copley et al. 2012). The small creeping patches experience initial rapid afterslip following the earthquake, but eventually stop sliding as the relaxed driving stress becomes balanced by the elastic resistance from the adjacent locked zones. The result is a more non-linear temporal evolution of ground deformation, which is localized adjacent to the co-seismic slip patch (see blue line in Fig. 5 and Supporting Information Fig. S3b). The misfit between the predicted and observed deformation at the end of the post-seismic transient as a result of afterslip along two predominantly locked faults is significantly smaller than for calculations with faults that can creep over their whole surface area (Fig. 5, Supporting Information Figs S3b and c). Therefore the InSAR measurements are better explained by a model in which the faults at Bam remain predominantly locked, compared to one in which the fault zone is able to creep. The implications of potential locking on the Bam fault zone for seismic hazard and fault behaviour are discussed in Section 7.1.

The stress-driven afterslip models we have described here have allowed us to predict the post-seismic slip distribution on the fault surfaces at Bam. In the following section we use these predictions to test whether co-seismic stress perturbations on the Bam-Baravat Fault would lead to an afterslip distribution consistent with growing features of the overlying topography over multiple earthquake cycles.

6 GROWTH OF THE BAM-BARAVAT RIDGE AND LITHOSPHERIC STRESSES

Formation of short-wavelength topography, like the sharp step in elevation across the Bam-Baravat Ridge (Fig. 2d), has in some cases been found to be associated with creep along shallow faults, caused by changes in the local stress field due to co-seismic slip on adjacent faults (Nishimura et al. 2008; Copley & Reynolds 2014; Fattahi et al. 2015; Copley & Jolivet 2016). The InSAR observations at Bam indicate that parts of the Bam-Baravat Fault began creeping after the 2003 earthquake (see Fig. 4 and Fielding et al.2009), consistent with this mechanism of growth. The topography, and lack of any geomorphic evidence for long-term right-lateral motion along the front of the Bam-Baravat Ridge, suggests that the majority of this creep is manifest as thrusting. However, the stress-driven afterslip calculations presented above predict little or no shallow motion on the Bam-Baravat Fault, and a slip vector that has a variable rake between the surface and 4 km (Fig. 10b), which is inconsistent with the InSAR observations and would not lead to long-term growth of the Bam-Baravat Ridge through repeat earthquakes and the associated shallow post-seismic creep.

The direction of the slip vector on the Bam-Baravat Fault is parallel to the maximum shear stress resolved in the fault plane, therefore is a function of both the co-seismic stress changes and the far-field tectonic stress acting through the lithosphere. As co-seismic stress-driven afterslip alone is not consistent with forming the local topography, we investigate whether a component of far-field tectonic stress can modify the post-seismic slip distribution on the Bam-Baravat Fault to reconcile the afterslip models with both short-term post-seismic, and long-term geomorphic, observations.

6.1 Forward calculations for variable Δσ_xx

Reverse motion along the N–S trending Bam-Baravat Fault will occur if there is an E–W oriented compressive deviatoric stress (Δσ_xx) acting through the lithosphere. We simulate this effect by performing stress-driven afterslip calculations as described above, and by adding a range of Δσ_xx as a pre-stress acting on the faults. We then repeat the grid search to determine the best-fitting afterslip model and friction parameters in a |$\chi _{\rm red}^2$| sense. The upper bound on the value of Δσ_xx will then be defined by the increase in |$\chi _{\rm red}^2$|⁠.

The calculations are performed for the case of both fully creeping and partially locked faults. In the grid search the value of aσ is fixed to be the nominal best fit of 5.5 MPa, because the misfit is relatively insensitive to this parameter (Fig. 10c). We analyse the modelled average shallow rake direction on the top 4 km of the Bam-Baravat Fault to track the effect of increasing Δσ_xx on the rake of post-seismic slip.

6.2 Forward calculations for variable Δσ_xx: results

The forward calculations suggest that for values of Δσ_xx greater than 10 MPa, the misfit between the models and the InSAR observations is significant, irrespective of the friction properties (Fig. 11a). The large misfits are generated by two effects: (1) Δσ_xx resolves mainly as a normal stress on the near-vertical Arg-e-Bam Fault, leading to clamping and limited strike-slip afterslip in the shallow fault zone, and (2) an increase in slip along the Bam-Baravat Fault as Δσ_xx resolves mainly as a thrust-oriented shear stress that drives fault creep, over-predicting the observed surface uplift near the Bam-Baravat Ridge.

Figure 11.

The effect of an E–W orientated compressive deviatoric stress (Δσ_xx) on the fit between models of co-seismic stress relaxation and the observed surface displacements. (a) shows the change in the reduced chi-square of the best-fitting model to the observed InSAR time-series for a range of Δσ_xx. (b) shows the associated evolution in the rake along the top 4 km of the fault zone, which is mainly independent of the locking pattern on the fault. The pink shaded region corresponds to the proposed lower and upper bound on the E–W deviatoric stress on the basis of: (1) requiring reverse-faulting motion on the Bam-Baravat Fault, and (2) low misfits between model and observations. The change in rake due to the applied deviatoric stress is shown in (c) which is the post-seismic slip distribution 5.7 yr after the Bam earthquake with Δσ_xx = 0, and (d) is the same but for Δσ_xx = 10 MPa. As in Fig. 10 the arrows correspond to the motion of the hanging wall relative to the footwall.

A trade-off exists between increasing Δσ_xx enough such that the post-seismic slip vector along the Bam-Baravat Fault is predominantly thrust oriented, in agreement with the geomorphic evidence of motion on the fault, and maintaining low misfits between models and the InSAR observations. This trade-off can be used to place bounds on the magnitude of the deviatoric stresses acting across the faults at Bam. Slip on the shallow portion of the Bam-Baravat Fault becomes predominantly thrust oriented when Δσ_xx ≳ 2 MPa (Fig. 11b). An approximate upper bound corresponds to Δσ_xx ≲ 10 MPa, above which there is a significant deterioration in the fit between modelled and observed LOS displacements (Fig. 11b; in particular on the ascending track where the models begin to significantly overpredict the positive LOS displacements over the Bam-Baravat Ridge; see Supporting Information Fig. S5).

6.3 Mechanical effects of fault surface anisotropy

An alternative explanation for the misfit between predictions of afterslip and the growth of long-term topography is that slip along the Bam-Baravat Fault may not be free to occur along any rake direction, as assumed in our models. If the surface of the fault had some component of topography or fabric (Candela et al. 2009; Sagy & Brodsky 2009), such as aligned corrugations (e.g. Jackson & McKenzie 1999), then slip will occur parallel to this fabric. As a result, irrespective of the irregularity of the maximum shear stress resolved onto the fault plane, relative motion of the two fault surfaces will always be thrust oriented and consistently lead to the growth of the overlying topography. Some corrugations could have formed during a previous tectonic regime, which could explain the mismatch between the predicted orientation of the slip vector from the stress-driven afterslip models and the orientation of the corrugations that would be required to guide thrust-oriented slip. It is not possible to distinguish between the effects of fault surface anisotropy or tectonic loading with our current observations.

7 DISCUSSION

7.1 Limited afterslip and future seismic hazard at Bam

Limited post-seismic relaxation of the co-seismically stressed portion of the Arg-e-Bam Fault between 10–15 km depth, either through afterslip or viscous flow, indicates that significant elastic stresses remain stored in the fault zone. The lack of recorded historical seismicity at Bam (Ambraseys & Melville 2005), despite the city being inhabited for over 2000 yr (Berberian 2005), suggests that there has not been recent failure of the Arg-e-Bam Fault between the 2003 co-seismic rupture and the base of the seismogenic layer. The estimated GPS loading rate of 1–3 mm yr⁻¹ (Supporting Information Fig. S1) (Walpersdorf et al. 2014) would allow between 2 and 6 m of slip deficit to accumulate on the deep portion of the Arg-e-Bam Fault over this 2000 year time period, which is sufficient to critically stress the fault, as revealed by the 2.7 m of slip in the 2003 earthquake. Release of the accumulated elastic strain in the deep fault zone could occur in future earthquakes that rupture the bottom part of the seismogenic layer, or long-term aseismic creep that is currently invisible in the InSAR measurements.

The size of a single earthquake that could release the deep elastic strains caused by co- plus interseismic loading at Bam would be around M_w 6. Earthquakes that ruptured the top and bottom half of the seismogenic layer in separate events have been documented elsewhere (Berberian et al. 2001; Elliott et al. 2011) and can be explained by loading of the unruptured part of the fault following the earlier earthquake. Of these examples, the Qaidam earthquakes were separated by 10 months (Elliott et al. 2011), and the Golbaf-Sirch earthquakes by 17 yr (Berberian et al. 2001), showing that failure of the loaded half of the seismogenic layer can occur on timescales much shorter than the recurrence interval for faults in central and eastern Iran, which is generally >1000 yr (Fattahi et al. 2006; Walker et al. 2010).

Alternatively, long-term creep on the bottom of the Arg-e-Bam Fault could release the elastic strain aseismically. As < 50 cm of slip at depth may have occurred over the 5.7 yr of InSAR measurements, and assuming creep had an approximately constant rate, it would take ≳30 yr to release the full component of stored elastic strain. This would also require the deep portion of the fault zone to have a significantly different rheology to the shallow fault zone, as transient creep on the top 5 km of the Arg-e-Bam Fault has a short relaxation time of 1.3 yr.

Analysis of the short-wavelength InSAR measurements that are sensitive to deformation in the shallow crust supports the proposition that the fault zone remains predominantly locked along-strike of the Bam rupture. We find that the spatial and temporal evolution of post-seismic ground deformation at Bam cannot be simulated by stress-driven afterslip on a fully creeping fault zone, which implies a significant component of locking on both faults.

In summary, it is likely that significant elastic strains remain stored around the fault zone at Bam both along-strike and down-dip, which could be released in future earthquakes of M_w 6, therefore the future seismic hazard remains high.

7.2 Post-seismic afterslip on the co-seismic rupture

In order to fit the post-seismic InSAR measurements, our stress-driven afterslip models require that afterslip overlaps with areas of the fault that also experienced co-seismic slip (see Supporting Information Fig. S4). As the InSAR measurements are sensitive to slip in the top 5 km of the fault zone, the spatial overlap appears robust. Superimposed co-seismic and post-seismic fault motion suggests that if the creeping part of the fault is subject to dynamic stresses above a certain threshold (e.g. due to co-seismic slip on an adjacent fault patch), then the creeping part may break in earthquakes. Otherwise, it will respond to quasi-static stresses through slow afterslip and interseismic creep.

7.3 Post-seismic relaxation time and fault friction

Post-seismic deformation observed at Bam is consistent with afterslip along discrete fault surfaces. Similar observations have been made for other earthquakes in Iran including at Tabas-e-Golshan (Copley 2014; Zhou et al. 2016) and Sefidabeh (Copley & Reynolds 2014). However, both these earthquakes show unusually long-lived afterslip, with deformation still visible in post-seismic interferograms for at least 15–30 yr following the earthquakes. In most documented cases surface deformation due to afterslip decays on a time-scale of months to years (Savage et al. 2005). Copley (2014) proposed two possible explanations for the observation of long-lived post-seismic transients at Sefidabeh and Tabas-e-Golshan: (1) the low the noise levels in the InSAR data acquired over Iran allow small signals to be recognized that may have been missed elsewhere, therefore long relaxation times are potentially more visible, or (2) these faults, and others with long-lived afterslip, may have unusual frictional properties. The observation that the post-seismic deformation at Bam is complete in less than ∼6 yr, despite the good InSAR observations and low atmospheric noise, suggests that measurement bias is not to blame for the differences in the extent of post-seismic deformation. Instead it seems that contrasts in the relaxation time between Bam, Sefidabeh and Tabas-e-Golshan are a result of differing frictional properties of the fault surfaces.

7.4 Strength of the Bam-Baravat fault

The 2–10 MPa deviatoric stress we have estimated to be acting on the Bam-Baravat Fault is the E–W compressive stress, minus the lithostatic pressure, averaged over 0–4 km. This result supports accumulating geophysical evidence that the total deviatoric stresses in the elastic portion of the lithosphere are on the order of earthquake stress drops (∼1–30 MPa; Lamb 2006; Suppe 2007; Copley et al. 2009, 2011). If the fault was supporting shear stresses consistent with Byerlee’s laboratory-derived static friction range of 0.6–0.85 (Byerlee 1978), then the depth-averaged deviatoric stress acting on the faults at Bam, if critically stressed, would be ∼130–200 MPa. This is one to two orders of magnitude greater than the results of our calculations. We therefore suggest that the creeping patches of the Bam-Baravat Fault are weak compared to traditional views of fault strength, which is consistent with recent suggestions for both stick-slip and creeping faults (Lamb 2006; Suppe 2007; Jolivet et al. 2013; Craig et al. 2014).

8 CONCLUSIONS

We have analysed the post-seismic ground deformation following the 2003 Bam earthquake in Iran, and find the observations are consistent with co-seismic stress-driven afterslip on a predominantly locked fault zone interspersed with small creeping patches. Afterslip and viscoelastic deformation between 10 and 20 km depth that is ‘invisible’ to the InSAR measurements may relax <20 per cent of the co-seismic stress changes on the deep fault zone, whilst the shallow fault experiences limited afterslip and appears to remain mostly locked. This suggests that the majority of the stress changes following the 2003 Bam earthquake have been stored as elastic strains along-strike and down-dip of the co-seismic rupture. Interseismically accumulated elastic strain on the unruptured patches of the fault, and the co-seismically induced strains, could be released in a future M_w 6 earthquake, hence the risk of future seismicity remains high.

We also find that the post-seismic slip distribution resulting from co-seismic stress changes along the creeping Bam-Baravat Fault following the 2003 earthquake is inconsistent with forming the topography of the overlying Bam-Baravat Ridge. This finding is interpreted to reflect two possible mechanisms: (1) that a 2–10 MPa compressive deviatoric stress acts through the lithosphere, leading to thrust oriented creep that extends to the surface, consistent with forming the long-term topography and in agreement with the post-seismic InSAR data (which would suggest that the shear stresses supported by the creeping Bam-Baravat Fault are low and therefore the fault is weak), or (2) the Bam-Baravat Fault may have some component of surface fabric that controls the direction of the slip vector.

Acknowledgments

This work forms part of the NERC- and ESRC-funded project ‘Earthquakes without Frontiers’, and was partly supported by the NERC large grant ‘Looking into the Continents from Space’. SW was partly supported by the BGS. SW thanks James Jackson, Andy Howell and Christoph Grützner for comments on the manuscript. We thank the editor, Duncan Agnew, Roland Bürgmann and one anonymous reviewer for their comments on the manuscript. Figures were produced using the Generic Mapping Tools (GMT) software (Wessel et al. 2013).

SUPPORTING INFORMATION

Supplementary data are available at GJIRAS online.

Figure S1. Results of calculations to estimate the locking depth on the faults at Bam. The black dots are the northward component of the GPS observations with 1σ error bars taken from Walpersdorf et al. (2014), whilst the coloured curves represent the best-fit velocity profile across the faults for a given locking depth D (Savage & Burford 1973). This analysis suggests the GPS velocity gradient across Bam is too small, and the GPS spacing too coarse, to provide robust estimates of the locking depth.

Figure S2. Calculations for the surface deformation caused by poro-elastic rebound projected into the satellite LOS. The poro-elastic rebound response was computed by finding the difference between the co-seismic surface deformation caused by the Funning et al. (2005) slip model using a drained and undrained Poisson ratio (ν_d = 0.25 and ν_ud = 0.28, respectively). The resulting surface displacements were then projected into the satellite LOS. The black dots represent the co-seismic surface fractures from the 2003 Bam earthquake (Jackson et al.2006). Note the different colour scale used to plot the LOS motion in this figure compared to the figures in the main text.

Figure S3. Calculations for the surface deformation resulting from complete relaxation of co-seismic stress changes on the faults at Bam through afterslip, performed with variable degrees of locking. (a) shows the ground deformation from the descending (top) and ascending (bottom) InSAR tracks covering the whole observation period (January 2004 to October 2009). This is taken to represent the total ground deformation resulting from stress-driven afterslip, as the rate of LOS change after 5.7 years of post-seismic deformation is 0 mm yr⁻¹ within uncertainty (see Fig. 5 in the main text). (b) shows the simulated ground deformation due to the complete relaxation of co-seismic stress changes through afterslip along a partially locked fault. Afterslip in this calculation was restricted to the shallow fault zone around the edge of the co-seismic rupture, recreating the effect of locking. The corresponding temporal evolution of deformation is shown in Fig. 5 in the main text. (c) shows an equivalent calculation to (b), but in this case the faults are allowed to creep unrestricted in response to co-seismic stress changes. The calculation in (c) clearly shows that unrestricted stress-driven fault creep would lead to extensive afterslip that is not consistent with our InSAR observations. Both calculations in (b) and (c) are independent of the frictional properties applied to the faults, as friction only controls the temporal evolution of deformation, and not the finite deformation in response to a given stress change.

Figure S4. The co-seismic slip distribution of Funning et al. (2005), with the contour of the compacted co-seismic slip model used as the input to the stress-driven afterslip modelling, and the locations of the three creeping patches used in the locking calculations overlain on top. The distribution of the creeping patches overlaps with regions of significant co-seismic slip (∼1 m), suggesting that at Bam there may be post-seismic afterslip in regions that also moved co-seismically.

Figure S5. Results of the best-fitting afterslip distribution for a given component of E–W deviatoric stress acting across the faults. (a–c) show the slip distributions corresponding the the best-fit models for Δσ_xx = 0, 2 and 10 MPa, respectively and highlight the rotation of the afterslip vectors towards pure thrusting. (d) shows the misfit between the best-fit stress-driven afterslip model and the descending track InSAR observations, with Δσ_xx = 0 MPa, and after 5.7 years of post-seismic deformation. (e) shows the same for the ascending track. (f–g) are the same as (d–e) except these are the misfits between the model with Δσ_xx = 10 MPa.

Table S1. Compilation of geodetic inversions used in determining the empirical scaling between co- and post-seismic geodetic moment release.

Table S2. Table of small-baselines descending track interferograms formed in this study with the date format in YYYYMMDD.

Table S3. Table of small-baselines ascending track interferograms formed in this study with the date format in YYYYMMDD.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.