7.2.1.1. History and recent progress in coseismic InSAR studies

Before InSAR, it was extremely challenging to measure surface displacements caused by earthquakes. Earthquakes in a specific location are rare, so taking displacement measurements relied on serendipity – earthquakes could only be measured if precise survey data had been acquired prior to the event (e.g. Stein and Barrientos 1985) or if natural data, such as uplifted shorelines, could be interpreted (e.g. Plafker 1965). As such, displacements from around 20 earthquakes or fewer had been measured (Wright 2000). Data from these studies appeared to confirm the elastic rebound theory popularized by Reid (1910) but were in general too sparse to be used to determine details of earthquake mechanisms and slip distributions.

The widespread availability of SAR data acquired before and after earthquakes, along with the proliferation of measurements from global navigation satellite systems (such as the US Global Positioning System), means that spatially dense, precise observations of surface displacements are now available for almost all recent earthquakes where there is significant (1 cm or more) movement of the land surface. Earthquake deformation has been mapped for events ranging in size from small (M_w ∼ 5), shallow (< 10 km) events on continental faults (e.g. Lohman and Simons 2005) to on-land deformation from great earthquakes on subduction megathrusts (e.g. Kobayashi et al. 2011).

The ability of InSAR to capture coseismic deformation relies on pre- and post-event data existing. Mission managers at the ESA recognized this opportunity and implemented systematic background acquisition plans for the ERS and Envisat missions (Attema et al. 2000; Miranda et al. 2013) that ensured data were acquired regularly over tectonically active areas, although the 35-day revisit time and relatively wide orbital tube for ERS and Envisat meant that the method did not work for many earthquakes. Nevertheless, by the end of the Envisat mission in 2010, studies on 60 earthquakes using InSAR data had been published (Weston et al. 2012).

The launches of the EU Copernicus Sentinel-1 C-band SAR satellites, Sentinel-1A (launched in 2014) and Sentinel-1B (launched in 2016), have radically increased the volume and quality of SAR data available for InSAR analysis. Together, the satellite constellation can acquire SAR data with a six-day revisit, and a small orbital tube means that interferograms can be formed between all acquisitions (Torres et al. 2012). Funning and Garcia (2019) examined Sentinel-1A data for 35 “potentially detectable” earthquakes that occurred between April 2015 and December 2016, finding that around 70% of these had coherent signals. However, they found that 31% of the events, mostly in vegetated low latitude regions, could not be detected with the available C-band data from Sentinel-1A. Since that study, Sentinel-1 has increased its duty cycle and acquires data globally every 6 or 12 days, depending on the location, including on both ascending and descending tracks in the tectonic belts. Several automatic systems have been built that provide rapid response interferometric data to the community for analysis (e.g. Barnhart et al. 2019; Lazeck`y et al. 2020).

C-band data (∼6 cm wavelength) can suffer from poor coherence in vegetated areas. Further improvements can be made in these cases using data from longer-wavelength L-band SAR missions (∼20 cm wavelength), such as from PALSAR on the Japanese ALOS-2 satellite (Rosenqvist et al. 2014). The launch of the joint US–India L-band SAR mission NISAR, currently scheduled for 2022 (Rosen et al. 2017), will ensure L-band SAR data are acquired globally for all ascending and descending passes with a 12-day revisit time. The EU Copernicus program is also planning an operational L-band SAR mission to complement the C-band data from Sentinel-1 (Elliott 2020).

7.2.1.2. What have we learned about faults from InSAR studies of coseismic deformation

For individual earthquakes on land, the detailed maps of surface deformation provided by InSAR typically provide more accurate locations and depths than are available from global teleseismic seismology, often revealing significant mislocations in the seismic catalogues (Weston et al. 2011). For shallow events larger than M_w ∼ 6, it is usually possible to derive detailed slip distributions from InSAR, again with the geodetic solutions resolving finer spatial details than is possible from global seismology. The large number of dense geodetic earthquake deformation observations has led to a significant increase in our knowledge of earthquake behavior in several key areas. Here, we summarize some of the major advances in our understanding about faults and earthquakes that have been made possible through InSAR observations.

Fault complexity: a key area where our understanding of earthquakes has changed because of the availability of precise, dense geodetic data from InSAR is earthquake complexity. Earthquakes with M_w less than around 6 can typically be modeled with uniform slip on a single rectangular fault plane in an elastic half-space (Okada 1985), particularly if they do not rupture the surface. However, larger, shallow earthquakes usually require more complex models of variable slip on segmented faults or multiple fault strands (e.g. Elliott et al. 2016). Several recent earthquakes, including the 2010 M_w7.1 Darfield (New Zealand) earthquake (Elliott et al. 2012), the 2010 M_w7.2 El Mayor Cucapeh (Mexico) earthquake (Huang et al. 2017) and the 2019 M_w6.4/7.1 Ridgecrest (USA) sequence (Xu et al. 2020), have required significant complexity to match the geodetic observations, with slip occurring on multiple fault strands.

The type example of a complex, multi-fault rupture is the 2016 M_w7.8 Kaikoura (New Zealand) earthquake, which ruptured more than 170 km of faults in the northeast part of the South Island of New Zealand (Figure 7.1; Hamling et al. 2017). In total, more than 20 crustal faults had surface ruptures, and these include faults that were previously mapped as well as previously unknown structures (Litchfield et al. 2018). Both crustal faults and the deeper subduction interface may have failed simultaneously (Furlong and Herman 2017). The earthquake defied conventional wisdom about the ability of fault segmentation to arrest/control rupture; the rupture appeared to jump across fault stepovers as high as 15–20 km (Hamling et al. 2017). Without the detailed geodetic imaging work from InSAR and optical correlation, it would have been extremely challenging for field teams to map the complexity of this earthquake – scarps from surface ruptures erode very rapidly in New Zealand. The new understanding of the potential for complex, multi-fault ruptures is now being incorporated in the latest generation of seismic hazard models.

Surface versus deep fault slip during earthquakes: another area where dense geodetic imaging data have changed our understanding of earthquakes is the relationship between fault slip seen at the surface in an earthquake and the slip that occurs at depth. Earthquakes on “blind” thrust faults, where minimal slip reaches the surface, were known before InSAR (e.g. Yielding et al. 1981; Ekström et al. 1992), but we now know that for all types of earthquake, the surface slip is in fact only a partial guide to the slip that occurs at depth. For example, modeling of the 1995 M_w6.4 Dinar earthquake, a normal-faulting earthquake in southwest Turkey, requires 1.2–1.4 m of slip as shallow as 2 km beneath the surface, but the maximum observed surface rupture magnitude was only 0.25 m (Fukahata and Wright 2008). Similarly, in the 2003 M_w6.5 Bam earthquake (strike-slip, Iran), models require 2–3 m of slip at ∼ 4 km depth, but the maximum surface slip was again 0.25 cm. This observation is not limited to moderate earthquakes – Kaneko and Fialko (2011) compared the slip profiles of the Bam earthquake and four other strike-slip faults with M_w > 7.1, showing that they all had a shallow slip deficit as high as 60%. Although more detailed analysis of near-field date by Xu et al. (2016) reduced the magnitude of the discrepancy for some of these earthquakes, it is clear that the surface slip often only provides a minimum estimate for the fault slip that occurred in the subsurface.

The observation of a discrepancy between surface observations and slip at depth has important implications for studies of paleoseismology – magnitudes estimated from paleo-offset data and empirical relationships between slip and magnitude (Wells and Coppersmith 1994) are likely to be lower bounds of the true magnitude of earthquakes. The causes and implications of the shallow slip deficit (SSD) have also been investigated. There are several plausible explanations. The first is that fault slip branches near the surface onto multiple strands, and without precise and detailed near-field data, this can map into an apparent deficit (e.g. Xu et al. 2016). A second possibility is that friction in the shallowest portion of the fault is velocity-strengthening – some shallow slip then occurs as afterslip or inter-seismic creep. Hussain et al. (2016) investigated this possibility for part of the North Anatolian Fault that ruptured in the 1999 Izmit/Kocaeli earthquake, finding that without additional periods of rapid shallow slip, this mechanism could only partially explain the observed shallow slip deficit. Finally, it is possible that inelastic distributed deformation, for example due to folding, accommodates some of the strain in the shallowest portion of the fault. It seems likely that these mechanisms may operate concurrently in some cases. High-resolution near-field displacement data from optical image correlation is invaluable for investigating these mechanisms (see discussion later in this chapter).

Dynamic triggering: although geodetic imaging techniques like InSAR provide a snapshot of the total deformation that occurs between two image acquisitions, the dense spatial coverage of the observations has revealed unexpected deformation patterns in a few instances that were initially missed by seismology. One of the best observations of dynamic triggering is from Pagli et al. (2003), who observed deformation on the Reykjanes Peninsula in Iceland that was consistent with an undetected M_w5.8 strike-slip earthquake. They showed that the earthquake was dynamically triggered by a M_w6.6 earthquake around 100 km away. The triggered earthquake was not initially detected by seismic networks because the seismic waveform was partially hidden by the waveforms from the larger event.

Another beautiful example of dynamic triggering was presented by Nissen et al. (2016), who investigated an M_w7.1 reverse-faulting earthquake that occurred in Pakistan in 1997. They found that the earthquake was in fact not a single event but a doublet, with two clear areas of uplift associated with separate faults. Joint analysis of the InSAR with seismology showed that the initial rupture was followed just 19 s later by rupture of a second reverse fault, located ∼ 50 km away. Again they invoke dynamic triggering to explain the occurrence of this doublet. Dynamic triggering is also likely to have played a part in the Kaikoura (New Zealand) earthquake (Figure 7.1) discussed earlier. The potential for dynamic triggering is an additional factor that needs to be accounted for in future seismic hazard models and scenarios.

Structural control of earthquake rupture: although some earthquakes, such as Kaikoura, appear to propagate through structural barriers (Hamling et al. 2017), InSAR has documented other examples where pre-existing structure does appear to play an important role in controlling earthquake rupture. Several good examples of this have been observed in the India–Asia collision zone. The 2001 M_w7.8 Kokoxili earthquake ruptured more than 400 km of the Kunlun Fault in northern Tibet, with slip reaching as high as 8 m in places. Lasserre et al. 2005 showed that areas of significant moment release corresponded to the first-order along-strike structural segmentation of the fault, with very low slip observed and modeled on some segment boundaries. Elliott et al. (2011) showed that depth segmentation can also act as a barrier to rupture – they investigated two M_w6.3 reverse-faulting earthquakes that occurred in the North Qaidam Thrust system (Tibet) at almost the same location but separated in time by nearly 300 days. They showed that the first of the two earthquakes ruptured the lower part of the seismogenic crust, with the later event rupturing a near continuation of the same fault plane in the upper seismogenic crust. They speculated that the gap between the two earthquakes was due to a structural barrier caused by the intersection of a reverse fault dipping in the opposite direction.

Structural control also appears to have played an important role in the 2015 M_w7.8 Gorkha (Nepal) earthquake, which partially ruptured the locked shallow Himalayan megathrust (Grandin et al. 2015; Elliott et al. 2016a). Qiu et al. (2016) used a structural geology model to show that the coseismic slip was confined to a flat region of the megathrust surrounded by two steeper ramps. They also ran dynamic simulations that showed that their fault geometry was sufficient to confine some (but not all) ruptures to the lower or upper part of the seismogenic fault; some earthquakes in their simulations were larger and ruptured the entire locked megathrust. Importantly for Nepal, the shallow, locked portion of the megathrust remains locked and has not moved following the earthquake (Ingleby et al. 2020) – the hazard remains very high there.

7.2.2.1. History and recent progress in inter-seismic InSAR studies

Mapping inter-seismic strain accumulation with geodetic methods has been a long-standing goal of the geodetic community as it has the potential to strongly inform seismic hazard models. If Reid’s 1910 elastic rebound hypothesis is correct then measuring where strain is accumulating, and at what rate, should allow us to quantify the long-term rates of fault movement, and hence, the rate of seismicity. Elliott et al. (2016) showed that in the western US and in eastern Turkey, areas with dense GNSS coverage and high tectonic strain rates, there is a good correlation between seismicity rate and geodetic strain rate, and data from the GNSS-derived Global Strain Rate Model (Kreemer et al. 2014) have been used to make global seismicity forecasts that compare well with forecasts based on instrumental seismicity (Bird and Kreemer 2015). However, dense GNSS data are not available globally. Here, we review progress in using InSAR to measure inter-seismic strain and show that it has the potential to significantly improve strain rate models built from GNSS data alone.

Using InSAR to measure inter-seismic strain is challenging compared to using InSAR for earthquakes – the inter-seismic signal in any one short-interval interferogram is much smaller than the noise, particularly the noise arising from spatial and temporal variations in tropospheric delays. The first studies to clearly show the signature of inter-seismic strain in InSAR results were able to reduce the noise in individual interferograms by stacking a relatively small number of datasets from ERS-1 and ERS-2 with multi-year time spans (Peltzer et al. 2001; Wright et al. 2001b). With the ESA operating Envisat with a more systematic background acquisition strategy in tectonic areas, alongside increasing volumes of data from systems such as ALOS, increasing computer power, and improved algorithms (e.g. Doin et al. 2011; Agram et al. 2013; Morishita et al. 2020), an increasing number of studies have been able to map inter-seismic tectonic deformation over ever larger and/or more challenging areas (e.g. Grandin et al. 2012; Wang and Wright 2012; Tong et al. 2013; Cakir et al. 2014; Pagli et al. 2014; Hussain et al. 2018; Wang et al. 2019).

The launch of the Sentinel-1 constellation has again provided a radical increase in the volume of SAR data acquired over tectonic areas (Elliott et al. 2015). The shorter revisit time increases coherence and means that average crustal velocities can be obtained with improved accuracy. Several groups have developed tools capable of mass processing Sentinel-1 data (e.g. Bekaert et al. 2019; Foumelis et al. 2019; Lazeck`y et al. 2020), with the ESA’s Geophysical Exploitation Platform (Foumelis et al. 2019) and COMET’s LiCSAR system (Lazeck`y et al. 2020) focused on mass producing interferometric products for the continental seismic belts. Weiss et al. (2020) presented the largest regional analysis of Sentinel-1 to date, processing Sentinel-1 data using LiCSAR and LiCSBAS (Morishita et al. 2020) over 800,000 km² of Anatolia. They produced average velocities for both ascending and descending tracks and used these in combination with GNSS data to provide east–west and vertical velocities and to estimate strain (Figure 7.3). The results reveal both large-scale features, such as the North Anatolian Fault, as well as localized deformation, for example due to fault creep or water extraction. Even in Anatolia, where GNSS coverage is good, the inclusion of InSAR data in the inversion of strain rate sharpens the strain rate model considerably (compare the North Anatolian Fault in panels (c) and (d) of Figure 7.3). The work in Anatolia is being expanded at present to produce complete velocity fields with this approach for the entire Alpine–Himalayan Belt and East African Rift.

7.2.2.2. What have we learned about faults from InSAR studies of inter-seismic deformation

Geodetic determination of fault slip rates: for individual faults, InSAR has regularly been used to obtain estimates of their present-day slip rate (e.g. Peltzer et al. 2001; Wright et al. 2001b; Cavalié et al. 2008; Wang et al. 2009; Karimzadeh et al. 2013; Fattahi and Amelung 2016). To do so requires a few assumptions – the rheology of the seismogenic layer is typically assumed to be elastic such that strain observed in the crust away from the fault is eventually be released by slip on the fault, and the rate of strain accumulation observed geodetically is typically assumed to be representative of the long-term strain rate. If we make those assumptions, geodetic slip rates can be estimated, alongside parameters such as “locking depth” that are indicative of the seismogenic layer thickness (Wright et al. 2013). Many of the geodetic slip rates, for example on the North Anatolian and San Andreas Faults, agree with estimates from geology/paleoseismology, but in a few cases, the geodetic slip rates have disagreed with the prevailing wisdom. For example, Wright et al. (2004) used InSAR to show that the faults in western Tibet were generally accumulating strain at a slower rate than had been estimated in some geological models (e.g. Tapponnier et al. 2001). More recent geological studies have yielded Quaternary rates that are in good agreement with the geodetic rates (e.g. Cowgill et al. 2009).

Beyond the rates and styles of deformation at individual fault systems, geodetic measurements of inter-seismic strain at faults have led to new insights into fault behavior in a number of key areas. We discuss these below.

Focused strain throughout the earthquake cycle for most faults: inter-seismic deformation, particularly at strike-slip faults, is often modeled using a simple elastic rheology, in which the deep part of a fault slips at a steady rate whereas the shallow part is locked. For strike-slip faults, a screw dislocation formulation can be used (Savage and Burford 1973); other fault geometries can be modeled with a backslip formulation (Savage 1983). These elastic dislocation models give rise to relatively focused areas of inter-seismic strain around active faults; in the standard screw dislocation model, 70% of the motion for a strike-slip fault, for example, occurs within a distance of two locking depths from the fault.

Wright et al. (2013) noted that these simple elastic dislocation models were ubiquitous in studies of inter-seismic deformation; their review found 187 examples in the literature. The fact that focused strain is seen at most faults suggests that it occurs at all phases of the earthquake cycle. This argument is supported by observations of late-stage inter-seismic deformation at several major strike-slip faults that later failed in earthquakes, including the North Anatolian Fault at Izmit (Turkey), the San Andreas Fault at Parkfield (California), the Denali Fault in Alaska and the Manyi Fault in Tibet (Elliott et al. 2016). This is a somewhat surprising result – if the viscosity of the lower crust is as low as some studies of post-seismic deformation have suggested (see below), then we would not expect focused inter-seismic strain late in the earthquake cycle (e.g. Hetland and Hager 2006).

There are a few exceptions, however, where faults do not appear to have clear focused strain associated with them during the inter-seismic period. A clear recent example is from Wang et al. (2019), who used ERS and Envisat data to map inter-seismic tectonic strain across south-central Tibet. Of the five extensional grabens in their study region, only two had clear focused extensional strain in the inter-seismic geodetic velocity field.

Temporal invariance of inter-seismic strain rate: typical inter-event times for major earthquakes on large faults are hundreds to thousands of years, so modern geodetic measurements are inevitably only available for a small portion of this period. To assess whether strain rates vary throughout this time period, indirect methods are required. Meade et al. (2013) compared short-term geodetic estimates of slip rate for 15 continental strike slip faults with long-term geological estimates of slip rate, estimated over multiple earthquake cycles, finding very good agreement overall between these rates. They argued, persuasively, that because the geodetic rates sample random short intervals within the earthquake cycle, the geodetic rates must be fairly constant throughout the earthquake cycle, apart from during a relatively short post-seismic period.

Hussain et al. (2018) took a different approach to the same problem, examining present-day strain rates along the entire North Anatolian Fault from InSAR alongside estimates of strain from GPS measurements before and after the 1999 Izmit earthquakes. As major earthquakes have occurred at different times in the past, these spatial samples of strain rate can be thought of as sampling different time periods of the earthquake cycle, assuming the entire fault behaves in a similar way. They found that strain rates were approximately constant apart from an approximately 10-year period following the 1999 earthquakes (Figure 7.4). This observation requires that the relaxation time for the majority of the mid to lower crust be as long as the inter-event time (Figure 7.4; Yamasaki et al. 2014; Hussain et al. 2018).

The fact that, for large strike-slip faults at least, the inter-seismic strain rate outside the post-seismic period appears to be fairly time invariant is good news for hopes to use geodetic observations to help constrain seismic hazard models. The primary input for most modern probabilistic seismic hazard models are instrumental and historical records of previous earthquakes. These can be biased, particularly in areas of moderate to low strain, because of the relatively short observation period compared to the inter-event time (see e.g. Stein et al. 2012). If short-term geodetic strain rates really are largely representative of long-term rates of fault activity, then they have the potential to be an important dataset for constraining earthquake hazard (e.g. Bird and Kreemer 2015).

7.2.3.1. History and recent progress in the study of aseismic deformation transients

As geodetic data from InSAR and GNSS have shown improvements in accuracy and sampling period, it has become increasingly evident that aseismic deformation transients are an important phenomenon in fault zones. They can occur as periods of accelerated post-seismic deformation following earthquakes, when the rate of deformation is higher than the inter-seismic period (e.g. Nur and Mavko 1974), or as “slow earthquakes”, where the aseismic deformation transients are not obviously associated with the earthquake cycle (e.g. Sacks et al. 1978).

The first post-seismic studies using InSAR came just a few years after the first coseismic studies by Massonnet et al. (1996) and Peltzer et al. (1996). Although these authors studied the same earthquake, their observations were limited to just a few post-seismic interferograms and they invoked different mechanisms to explain the data. A number of studies in the 2000s presented coarsely sampled time series showing how post-seismic deformation evolved through time (e.g. Fialko 2004; Ryder et al. 2007; Biggs et al. 2009). Prior to the launch of Sentinel-1 in 2014, there were around 23 earthquakes with geodetic observations of post-seismic deformation, of which more than half had been detected by InSAR. The latest generation of InSAR missions, including Sentinel-1, allow the spatial and temporal history of post-seismic deformation to be mapped with unprecedented detail (e.g. Figure 7.5; Floyd et al. 2016; Ingleby et al. 2020; Wang and Bürgmann 2020).

Aseismic deformation has also been observed with InSAR. Some of the earliest InSAR observations of aseismic deformation were for creeping sections of major faults, including the Hayward and San Andreas faults in California (Bürgmann et al. 2000; Johanson and Bürgmann 2005) and the North Anatolian Fault in Turkey (Cakir et al. 2005). With more recent data it has been possible to map the spatial and temporal distribution of creep, showing that it often accrues in aseismic transient slip events (“slow earthquakes”). Furuya and Satyabala (2008) were the first to present InSAR observations of a slow earthquake, for the Chaman fault system in Afghanistan – lasting around a year and with an equivalent moment magnitude of 5.0. By using time-series analysis, several authors have now been able to map the spatial and temporal distribution of creep events at continental strike-slip faults (Rousset et al. 2016; Pousse Beltran et al. 2016; Khoshmanesh and Shirzaei 2018). Cavalié et al. (2013) presented the first observations from InSAR of a slow slip event on a subduction zone, at Guerrero, Mexico, but in general, the deformation caused by slow slip events at subduction zones has been too small to be reliably detected using InSAR.

7.2.3.2. What have we learned about faults from InSAR studies of aseismic deformation transients

Post-seismic deformation is now commonly observed: for most moderate to large earthquakes, post-seismic deformation is now routinely observed – it is now recognized as an important part of the earthquake deformation cycle. However, multiple mechanisms have been invoked to explain the observations (Wright et al. 2013), notably frictional afterslip (e.g. Floyd et al. 2016), viscoelastic relaxation (e.g. Ryder et al. 2007) and poroelastic rebound (e.g. Jónsson et al. 2003). Often, different mechanisms are invoked for the same earthquake from geodetic datasets with different time or spatial scales (Wright et al. 2013). For the 1992 Landers earthquake, post-seismic deformation has been attributed to afterslip and fault zone contraction (Massonnet et al. 1996), poroelastic deformation (Peltzer et al. 1996), deep afterslip (Savage and Svarc 1997), a combination of poroelastic deformation and deep afterslip (Fialko 2004), linear Maxwell viscoelastic relaxation of the lower crust (Deng et al. 1998) and power-law viscoelastic relaxation of the mantle (Freed and Bürgmann 2004).

Multiple processes contribute to post-seismic deformation: when looking at individual earthquakes with detailed time series of geodetic data, such as was possible for the 2014 M_w 6.1 South Napa earthquake by combining GNSS with Sentinel-1 InSAR, the spatial and temporal distribution of deformation appears as somewhat complex (Figure 7.5; Floyd et al. 2016); afterslip appears to start at different times in different locations and have different temporal decay patterns (compare the slip histories in panels (a) to (c) of Figure 7.5). To try to make sense of the complexity, Ingleby and Wright (2017) compiled time histories of post-seismic velocity for all observations of post-seismic deformation from continental earthquakes (151 velocity observations from 34 events). They showed that the maximum post-seismic velocities are inversely proportional to the time since the earthquake, like Omori’s law for aftershocks (Utsu et al. 1995). This requires frictional afterslip or a highly nonlinear viscoelastic deformation of a shear zone (Ingleby and Wright 2017). However, the maximum post-seismic velocities studied by Ingleby and Wright (2017) are likely to be dominated by shallow fault zone processes, such as afterslip. Studies of large earthquakes show that a combination of frictional afterslip and deeper viscoelastic deformation is required (e.g. Weiss et al. 2019). To separate the contributions of these mechanisms requires long time periods of precise observations over large spatial scales.

A wide range of complex aseismic phenomena occur: in addition to post-seismic deformation, InSAR has helped identify a range of complex and varied aseismic creep behavior, including triggered slip, fault creep and slow earthquakes. Triggered aseismic slip of nearby faults has been observed following a number of major earthquakes (e.g. Wright et al. 2001a; Fielding et al. 2004; Fujiwara et al. 2016). Xu et al. (2020) presented a compelling analysis of 169 fractures away from the main faults that deformed during the 2019 Ridgecrest earthquake sequence (M_w 6.4, 7.1, California). They showed that fractures that experienced coseismic stress changes that were in the same sense as the regional stress field had narrow deformation zones and were likely to have failed due to permanent inelastic deformation, such as fault slip. Faults that were stressed coseismically in the opposite sense to the regional stress had wider deformation zones, with deformation likely to be the result of enhanced strain in a compliant fault zone.

Detailed studies of time series of aseismic creep are starting to reveal complexities in the temporal behavior of creep. For example, Rousset et al. (2016) used a dense time series of high-resolution InSAR data to document a pulse of fault creep that took place over around a month in 2013. The pulse had a magnitude of 15 mm and was the only significant fault motion to occur in their observation window (July 2013 to April 2014). Similarly, Khoshmanesh et al. (2015) showed that creep at the Parkfield section of the San Andreas Fault was not steady in time, instead occurring in pulses separated by one to three years. Khoshmanesh and Shirzaei (2018) suggest that the pulsing is caused by variations in normal stress that result from pore fluids trapped in hydraulically-isolated parts of the fault.

Aseismic behavior is related to lithology: spatial maps of aseismic fault creep have been powerful in illuminating the close relationship between aseismic creep and fault zone lithology. In Taiwan, Thomas et al. (2014) showed that creep occurred on the Longitudinal Valley Fault in Taiwan in areas where fault zone gouge was dominated by a particular lithological unit, derived from unlithified subduction zone mélange, that favored pressure solution creep. In Turkey, Cetin et al. (2014) showed that creep on the North Anatolian Fault is fastest in areas where the lithology has low frictional strength, including andesitic-basaltic, limestone and serpentine units caught up in the fault zone. The relationship between aseismic deformation and lithology has also been observed in post-seismic data – regions of shallow afterslip that occurred following the 2014 South Napa earthquake correlated with younger sediments and alluvium; Floyd et al. (2016) argued that the complex variations in spatial and temporal behavior (Figure 7.5) may be a function of different clay contents, mineralogy or pore pressures within these units.

As the quality of InSAR time series improve globally, we expect to see further advances in our understanding of the aseismic behavior of faults.

7.3.2.1. Optical image correlation: general workflow

Measuring the horizontal surface displacement field from raw optical images requires three main steps: (1) a pre-processing step during which the pre- and post- earthquake images are co-registered and orthorectified, (2) the correlation of the pre- and post-orthoimages and (3) the post-processing of the output surface displacement maps. Some optical datasets are provided in a processed format (with radiometric and geometric corrections applied), e.g. Landsat 8, and Sentinel-2. In such cases, we may proceed directly to the correlation without any pre-processing; however, any artifacts arising from the pre-processing stages cannot be easily corrected from the processed imagery (for example, DEMs used by space agencies to orthorectify images are often lower resolution than the images, which can lead to related topographic noise in the final displacement map). If datasets are provided without geometric correction, the user must pre-process the data prior to correlation. For more details about the pre-processing and the correlation steps, we refer the reader to Chapter 2 of this book as well as to the review papers by Avouac and Leprince (2015) and Hollingsworth et al. (2020).

Basic post-processing to minimize the influence of noise on the output surface displacement fields typically includes outlier removal, based either on displacements with unrealistic values, low signal-to-noise ratio, or values which differ significantly from their neighbors (e.g. Zinke et al. 2019), and various spatial smoothing operations (e.g. low-pass, median or non-local mean filtering). Non-local mean filtering, which is a spatial averaging filter with weights applied based on the similarity with each target pixel, is particularly appropriate for denoising displacements containing sharp discontinuities, such as those produced by surface-rupturing earthquakes, because it preserves sharp edges, and thus preserves fine details of the surface rupture (Buades et al. 2005). More specific corrections can be applied depending on the case, such as the removal of ramps, stripes, or undulations due to misregistration errors, charged-couple device misalignment or unknown variations of the satellite attitude. More details about these errors and their corrections can be found in section 7.3.2.3.

Several software packages are available to process optical images. The most widely used by the earthquake science community include: COSI-Corr (Co-registration of Optically Sensed Images and Correlation (Leprince et al. 2007b)), MicMac (Multi Images Correspondances par Méthodes Automatiques de Corrélation (Rosu et al. 2015)) and ASP (Ames Stereo Pipeline (Beyer et al. 2018)).

7.3.2.2. Measuring the 3D displacement field of earthquakes from stereo optical images

With the increasing number of stereo optical images available, it is now possible to measure the vertical component of the displacement in addition to the horizontal component. Three different methods can be used; all of them require pre- and post-stereo pairs of optical images³ (with enough stereoscopic information to isolate the pre- and post-earthquake topographic signal from the tectonic signal).

DEM differencing: the most straightforward technique is to (1) compute pre- and post-earthquake DEMs from pre- and post-stereo optical images, (2) orthorectify the pre- and post-images using the relevant DEMs, (3) correlate the orthorectified images to obtain the horizontal displacement field, (4) resample the pre-DEM to incorporate the horizontal displacement and (5) compute the vertical difference between the two DEMs. If the horizontal displacement is not accounted for when differencing DEMs, the resulting vertical displacement field will be overprinted by noise correlated with misregistered topography (e.g. Oskin et al. 2012). If the amount of horizontal displacement is small enough, or in areas of smooth topography, then the pre- and post-DEMs can be simply subtracted without accounting for the horizontal displacements (Kuo et al. 2019). This technique has been successfully applied to measure the 3D displacement field of several earthquakes, such as the Norcia, Hualien and Balochistan earthquakes (Barnhart et al. 2019a; Kuo et al. 2019; Delorme et al. 2020).

3D triangulation: in the second approach, the 3D location of each point before and after the earthquake is triangulated from displacement maps computed for the stereo pre- and post-event images (Avouac and Leprince 2015; Zinke et al. 2019). Let us consider two stereo pairs of optical images covering the same area: one taken before the earthquake and the second one after (Figure 7.6). One of the pre-earthquake images is defined as the primary (i.e. reference) image, and three disparity maps are computed for the primary image and the three secondary images (Figure 7.6). As stereo images are taken with off-nadir incidence angles, the pre-earthquake disparity map contains only topographic information, while the two pre-/post-earthquake disparity maps contain a combination of stereoscopic effects and ground displacements from the earthquake (Figure 7.6). Using a ray-tracing approach, the x,y and z positions at each grid location in the pre-event correlation map can be retrieved using the pre-event correlation map and the camera positions of the pre-stereo images. Similarly, the x,y and z positions of these same points after the event can be retrieved using the post-event camera positions, and the pre-/post-correlation maps. The 3D displacements are then obtained by simply differentiating the pre- and post-3D positions. An example of computing a 3D displacement field is given in section 7.3.3.2.

Point-cloud registration: finally, the third approach is to use point-cloud registration algorithms (e.g. the iterative closest point and fast global matching algorithms) on pre- and post-earthquake DEMs computed from stereo optical images. This method allows more expressive transforms to be retrieved, albeit with limited subpixel precision, thereby revealing, in the case of an affine transform, local changes in rotation, scale and shear, in addition to rigid 3D translations (e.g. Nissen et al. 2014). The advantage of using this technique instead of image correlation-based approaches is that the similarity metric is based on surface shape/morphology rather than surface reflectance. Therefore, optical images with vastly different reflectance properties (e.g. agricultural settings, regions with strong seasonal changes in vegetation and snow) are likely to fail with traditional correlation. However, if the pre- and post-event DEMs can be well-resolved for each epoch and precisely co-registered, point-cloud registration techniques can be used to retrieve the 3D ground deformation. Point-cloud registration is commonly used to measure the 3D surface deformation with LiDAR datasets (e.g. Nissen et al. 2014; Scott et al. 2018), which are generally high resolution (negating the need for subpixel performance), very high quality and well-resolved at high frequencies. This is not always the case with DEMs derived from stereo imagery, although with increasing resolution, image quality and availability of satellite and aerial datasets, and tools for processing them, this technique is finding increasing applications with optical datasets (e.g. Barnhart et al. 2019a).

7.3.2.3. Limitations and sources of error

Error and bias in surface displacement fields obtained with optical image correlation have a variety of sources. First, errors can arise from imperfect modeling of the imaging system. For example, systematic misregistration can produce a ramp over the entire displacement field. Moreover, unrecorded variations of the satellite attitude, such as pitch oscillations, can lead to across-track undulations in the displacement field (e.g. Dominguez et al. 2003; Leprince et al. 2007a), while misalignment of the charge-coupled devices (CCDs) of pushbroom sensors can cause along-track stripes (e.g. Leprince et al. 2007b). These artifacts are very common (Avouac and Leprince 2015) and can be mitigated using detrending and destriping tools (e.g. Ayoub et al. 2017). When using film-based images, unaccounted for distortions of the film (e.g. thermomechanical warping) and scanning artifacts can also bias the correlation maps (e.g. Michel and Avouac 2006; Ayoub et al. 2009). Bias affecting the entire displacement field (e.g. orbital error and film distortion) make the measurement of the absolute displacement difficult. Indeed when such errors affect the correlation, only the measured offsets across the fault are reliable whereas the absolute values of displacements are not. These errors can be mitigated by setting the far-field displacements to zero (e.g. Vallage 2015). However, when working with high-resolution images, the image coverage is too tight, so there are no zones of zero displacement. In that case, a possible solution is to work with images at different scales, anchoring high-resolution correlation with lower-resolution correlation but larger coverage. If GPS data are available, they can also be used to estimate a correction to apply to the whole optical displacement field (e.g. Dominguez et al. 2003). Topographic artifacts can appear when images are orthorectified with a lower resolution DEM, if the same DEM is used to orthorectify both pre- and post-images despite a change in topography during the earthquake, and if images are not well registered to the topography (e.g. Ayoub et al. 2009). Finally, any changes of the ground surface between the pre- and post-event images will impact the resulting displacement field. The correlation will fail in areas of drastic change, such as those due to urbanization, farming, snow, clouds or slope failure, while less drastic changes will likely be expressed as coherent noise in the correlation map. For example, if the pre- and post-images are not taken at the same time of the year, shifts in the positions of shadows between the pre- and post-acquisitions will be resolved by the correlator, leading to coherent noise in the surface displacements. This bias can be avoided by using pre- and post-images taken at the same epoch of the year (i.e. with the same illumination conditions). Alternatively, shadowing bias can also be characterized and corrected using time-series analysis (Lacroix et al. 2019; Hollingsworth et al. 2020). Indeed, because the position of the sun (i.e. azimuth and elevation) at the time of image acquisition changes throughout the year, the orientation and length of the shadows will change locally, depending on the topography. Therefore, the displacement bias induced by the shadows will be different depending on the time of the year at which the pre- and post-images were taken. By correlating all possible combinations of images taken throughout the year, the redundancy of the measurements can be used to invert a robust displacement time series for each pixel. As the time series obtained is the sum of various signals (e.g. tectonic signal + illumination bias + other noise), the isolation and removal of the illumination component can be achieved using, for example, blind source separation approaches (e.g. independent component analysis) or decomposition into known basis functions (if their functional form is known).

7.3.2.4. Perspectives and future developments

One of the main limitations of the optical correlation technique is the difficulty of measuring small displacements (cm to mm) produced by small earthquakes or by transient processes such as post-seismic deformation or inter-seismic creep (especially if the displacement signal at the surface has a long spatial wavelength). In the near future, optical images with increasing spatial and temporal resolutions will be increasingly available, along with improved correlation schemes. This large archive of optical data will thus allow measurement of displacements of increasingly smaller amplitude and with higher temporal resolution. In particular, we will be able to make use of the large archive of optical images to construct dense time series of optical displacements. Indeed, the correlation of every possible combination of optical image available between two dates will provide a set of displacement measurements with high redundancy that can be inverted for a robust time series. Time-series analysis is already used to measure small displacements using InSAR data and has started to be used with optical data to monitor slow-moving landslides and glaciers (which typically experience larger displacements compared to those produced throughout the earthquake cycle). With the increase in temporal and spatial resolution of optical imagery, we might hope that such processing methods will enable the measurements of transient deformation relating to the seismic cycle, such as post-seismic deformation or inter-seismic creep. Furthermore, machine-learning approaches applied to optical geodesy will help to improve correlator performance (faster, increased robustness and better accuracy) and may be used to detect and correct sources of noise coming from the imaging system (CCD and jitter), thereby increasing our ability to robustly measure small displacements.

7.3.3.1. Horizontal surface displacement of the 1992 Mw 7.3 Landers earthquake from high-resolution aerial images

As a first application example, we describe the horizontal displacement field of the 1992 Mw 7.2 Landers, California, earthquake measured from the correlation of aerial images. The pre-earthquake dataset was composed of ∼40 National Aerial Photography Program (NAPP) overlapping photographs of 1 m resolution acquired in July 1989 (i.e. two years before the earthquake). The post-earthquake dataset was composed of nearly 400 photographs of 20 cm resolution taken two days after the earthquake. The post-earthquake dataset was acquired specifically to map the rupture trace of the Landers earthquake and was composed of several segments of overlapping photos that follow the surface rupture (Figure 7.7). Due to the very high resolution of the images, the spatial coverage is very narrow (image width of ∼2 km, Figure 7.7). The aerial photos of each dataset were acquired with a frame camera, meaning that the incidence angle increases radially outwards from the optical center. Partially overlapping images therefore have different incidence angles and can be used to extract a DEM for the overlapping region. Therefore, for each dataset, DEMs were extracted and used to orthorectify the images. Each segment of the post-earthquake dataset was processed independently. The only metadata available for the post-images were the approximate latitude/longitude coordinates of the image center and the focal length of the camera. Therefore, using the scale of the images (estimated from several GCPs), the camera elevation was estimated and the camera information (position and orientation) refined using the bundle adjustment tool in ASP to ensure that the overlapping images were consistent with each other (i.e. to make sure that the x,y and z coordinates of features seen in several overlapping images were identical). The DEM was then extracted and registered to a reference DEM (used as ground truth) using the iterative closest point (ICP) alignment tool in ASP (“pc_align”). We used a 5 m DEM covering the entire study zone (extracted from stereo Worldview satellite images). The alignment step with the reference DEM ensured that each segment was well-registered to the same external reference and therefore consistent with neighboring segments. Once a satisfactory alignment was achieved (the alignment could be evaluated by simply differencing the aerial and reference DEMs), the aerial DEM was used to orthorectify the aerial photos. Finally, the orthoimages belonging to a given segment were mosaiced together.

The same processing chain was used for the pre-earthquake images except that all images were processed simultaneously (they were all overlapping). Both the pre- and post-images were orthorectified at 1 m resolution to match the lower resolution of the pre-earthquake images. Then, the pre- and post-orthorectified images were correlated using the phase correlator of COSI-Corr with an initial window size of 64×64 pixels, a final window size of 32×32 pixels and a step size of 8 pixels. This resulted in a displacement field with ground sampling every 8 m and truly independent displacements every 32 m (although the windowing function used in the correlator resulted in nearly independent measurements every ∼16 m). Finally, measurements with a low signal-to-noise ratio (< 0.97) and outliers with displacement values higher than ± 10 m were discarded and the EW and NS displacement fields were smoothed using the non-local-mean filter implemented in COSI-Corr (Buades et al. 2005).

The surface displacement field is shown in Figure 7.7. The processing chain applied (bundle adjustment and alignment with a reference DEM) allowed a seamless horizontal displacement map to be obtained (Figure 7.7). The high-resolution displacement field allowed detection of secondary ruptures, which are clearly visible as sharp discontinuities yet were only partially mapped in the field (Figure 7.7, zoom 1). From these measurements, we can also characterize a wide zone of complex deformation that occurred north of the Kickapoo step-over (Figure 7.7, zoom 2).

7.3.3.2. Three-dimensional surface displacement of the 2016 Mw 7.8 Kaikoura earthquake from Worldview satellite images

On November 13, 2016, an Mw 7.8 earthquake struck the southeast coast of the South Island of New Zealand, 60 km southwest of the city of Kaikoura. This earthquake produced the most complex surface rupture ever observed, while the source mechanism indicated about equal thrust slip and right lateral strike-slip motion. Field investigations measured surface strike-slip up to 12 m and vertical slip up to 6.5 m. Considering the significant role that thrust faulting played in this earthquake, computation of both horizontal and vertical surface displacements was essential to sufficiently characterize the kinematics of this earthquake. Using stereoscopic Worldview 1–3 images of approximately 0.5 m resolution, Zinke et al. (2019) measured the 3D displacement field of the Kaikoura earthquake using the method describe in section 7.3.2.2 (3D triangulation approach). Before proceeding to the correlation, they first corrected sub-pixel artifacts in the raw images (from CDD array mis-alignments) using ASP. Then, for each overlapping pair of preand post-stereo images, the pre-earthquake image with an incidence angle closest to nadir was selected as the primary (to maximize the correlation with the three other images). Then, the three secondary images were coregistered to the master pre-images, and the four images were orthorectified using a coarse DEM (8 m) to remove the long-wavelength topographic distortions (thus allowing the remaining high-frequency stereoscopic distortions to be retrieved with a relatively small 32×32 correlation window). Then, the primary image was correlated with the three secondary images to compute disparity maps that contained information on both the short-wavelength (<8 m) topographic distortions and the tectonic displacement. Finally, they used a ray tracing approach between the disparity maps and satellite position to obtain the pre- and post-event 3D positions of each point of the ground. The 3D displacement was then obtained by simply differentiating the pre- and post-3D positions of each point of the ground (Zinke et al. 2019). ICP alignment was then used to align each correlation map with its neighbor, and the full displacement map was then aligned to either an independent correlation map (e.g. derived from correlation of Landsat 8 images) or GPS data.

The 3D displacement field of the Kaikoura earthquake is shown in Figure 7.8 and highlights a very geometrically and kinematically complex rupture composed of 19 faults, including one fault previously unmapped structure, the Snowflake Spur fault (Zinke et al. 2019). The surface slip amplitude is highly variable, with significantly larger slip in the northern region, where the slip is subparallel to the plate motion direction, than in the southern part, where the slip is subperpendicular to the plate motion direction (Figure 7.8).

7.3.3.3. Horizontal surface displacement field of the 1979 Mw 7.1 Khuli-Boniabad, Iran, earthquake from KH-9 spy satellite images

The optical image correlation technique allows us to take advantage of the substantial archive of historical aerial and satellite images to document the surface displacement fields of past earthquakes. The oldest events that have currently been studied with optical data occurred in 1975–1984 (Krafla seismo-volcanic rifting crisis crisis (Hollingsworth et al. 2012, 2013)), 1978 (Tabas earthquake (Zhou et al. 2016)) and 1979 (Khuli-Boniabad earthquake (Marchandon et al. 2018a)); nevertheless, any earthquake can be studied for which pre-imagery exists, potentially extending as far back as the beginning of aerial photography in the early 20th century. Correlation of aerial photos therefore provides a powerful method to considerably increase the database of documented surface ruptures, which are important for fault mechanics and seismic hazard research (as in the SURE database). For this purpose, the declassified KH-9 US spy satellite images are particularly useful. The KH-9 spy satellite was active between 1973 and 1980 and acquired ∼29,000 images globally with a footprint of 250×125 km and a resolution of 6–9 m. The KH-9 negatives are scanned by the USGS at 7 µm (3,600 dpi) using a high-performance photogrammetric film scanner, and the digital copies can be purchased from the USGS Earth Explorer website for USD 30. Although KH-9 image metadata have not been declassified, KH-9 images were acquired with a camera similar to the Large Format Camera (LFC) used by NASA in 1984, which featured a 23×46 cm frame camera with a focal length of 30.5 cm. For the purposes of orthorectification, they can be treated as very high altitude aerial photos (Surazakov and Aizen 2010; Hollingsworth et al. 2012). In the example presented in this section, the surface displacement field of the 1979 Mw 7.1 Khuli-Boniabad, Iran, earthquake was measured using a KH-9 pre-earthquake image and a 10 m Sentinel-2 post-earthquake image (Marchandon et al. 2018a; Figure 7.9). A Sentinel-2 image taken in the same season (i.e. January) as the KH-9 image was used to limit the illumination bias in the correlation (Hollingsworth et al. 2017; see also section 7.3.2.3). Using the COSI-Corr software package, the KH-9 image was coregistered and orthorectified using the post-earthquake Sentinel-2 image as reference and a 0.5 m resolution Pléiades-derived DEM of the studied area. A very long-wavelength trend was removed from both components of the EW and NS displacement fields by subtracting a second-order polynomial using the detrending tool of COSI-Corr. Such distortions are commonly seen when correlating film-based images and are usually due to thermo-mechanical warping of the photographic film (Michel and Avouac 2006). Nevertheless, the sharp discontinuity associated with the rupture was unaffected by long-wavelength denoising.

The left-lateral surface rupture trace is clearly visible on the EW displacement field (Figure 7.9). The 50-km long rupture trace trends ∼N70°, before bending abruptly N160° (west of point F) at its western termination. Surface displacements in the western part of the rupture (west of point F) are typically higher than in the eastern part of the rupture (between points B and C). This is shown by the two example profiles P1 and P2, which reveal fault offsets of ∼ 4.92 ± 0.73 m and ∼ 2.24 ± 0.26 m.

7.3.4.1. Characterization of off-fault deformation

During an earthquake, a fraction of the coseismic strain is released off the main fault, either as localized slip on minor secondary faults or as inelastic deformation expressed as warping, granular flow or rigid-block rotation (Figure 7.10; Nelson and Jones 1987; Shelef and Oskin 2010). Given its nature, coseismic distributed deformation is difficult to observe and quantify during field surveys. Indeed, field offset measurements are typically made within 10 m of the primary fault strand and thus capture only the on-fault component of the deformation. Some field studies have also attempted to capture the coseismic distributed deformation. For example, Rockwell et al. (2002) used linear features crossing the fault (e.g. tree alignment, fences, walls or canals) to measure deformation accommodated off the main fault by inelastic processes during the 1999 ˙Izmit and Düzce earthquakes. They found that up to 40% of the total slip was accommodated off the primary fault plane over lateral distances up to 100 m. The same approach can be taken with visual inspection of high-resolution optical images (e.g. Zinke et al. 2014; Gold et al. 2015). However, estimation of off-fault deformation in this way can be time-consuming and interpretive and may only be achieved in areas where long linear features cross the fault.

By capturing the pattern of near-field deformation, optical correlation allows measurement of the total shear accommodated across the fault at high spatial resolution. Optical displacement fields can thus be used to characterize the off-fault deformation (e.g. magnitude, spatial distribution and width) in recent surface-rupturing earthquakes, and to evaluate which parameters control its characteristics (e.g. Zinke et al. 2014, 2019; Milliner et al. 2015; Vallage et al. 2015). The total surface slip accommodated across the fault zone, along with the width over which it is accommodated, are usually measured from fault-perpendicular stacked profiles spaced regularly along the fault strike (e.g. Zinke et al. 2014; Milliner et al. 2015). On each side of the displacement profile, a linear regression is fitted from the far-field to the point where the displacement profile shows an inflection near the fault (care must be taken to account for any near-fault artifacts introduced by smoothing operations, as well as correlated noise). Then, the linear regressions are extrapolated to the fault trace and the relative difference between the extrapolated regressions defines the offset across the entire fault zone, while the distance between the inflection points located each side of the fault defines the fault zone width (Figure 7.10). The amount of off-fault deformation is then quantified by subtracting the on-fault slip (e.g. offsets measured in the field) from the total slip measured from optical displacement maps (i.e. total slip across the entire fault zone), and expressed as a percentage of the total slip. The method relies on the assumption that the differences between optical and on-fault offsets reflect only the degree of optical frequency division (OFD) and not a reduction of localized slip on the primary fault plane in the uppermost tens of meters of the crust. A comparison between the along-strike offsets measured from optical displacement maps and those measured within a few meters of the main fault strand (by inspection of satellite imagery) of the 2013 Mw 7.7 Balochistan earthquake (Zinke et al. 2014) is shown in Figure 7.11.

The average amount of off-fault deformation inferred from high-resolution near-fault displacement fields for various earthquakes ranges from 20% to 60%, while the mean fault-zone width ranges from ∼100 to ∼400 m (Zinke et al. 2014; Gold et al. 2015; Milliner et al. 2015; Kuo et al. 2019; Scott et al. 2018; Zinke et al. 2019). It should be noted, however, that the data and method used to estimate OFD can provide significantly different results for a given earthquake (e.g. 45% vs. 28% for the Balochistan earthquake in Zinke et al. (2014) and Gold et al. (2015), respectively). The parameters controlling the amount of OFD from one rupture to another are not completely understood, although fault zone maturity is thought to play a significant role: with increasing fault maturity (e.g. increasing cumulative offset, see Manighetti et al. (2007)), the fault becomes geometrically smoother (Brodsky et al. 2011), while developing wider and more enhanced damage zones (Faulkner et al. 2011) and more established (and mechanically weak) gouge zones, ultimately leading to increased localization of slip during an earthquake (Dolan and Haravitch 2014; Socquet et al. 2019). Cumulative long-term OFD has also been estimated using long-term markers of deformation; for example, Salyards et al. (1992) measured the paleomagnetic rotation of sediments at a site on the San Andreas fault and found that as much as 78% of the deformation was distributed over a fault zone width of ∼45 m. Shelef and Oskin (2010) used geologic features (i.e. dikes and inactive faults) offset by active strike-slip faults in the Eastern California Shear Zone to conclude that distributed displacement accommodates up to 25% of the total displacement within 100–200 m of the studied faults; such values are also consistent with coseismic estimates from a structurally mature section of the 2016 Kaikoura earthquake (Zinke et al. 2019), with occurring within ∼100 m of the fault zone.

Along the length of a surface rupture, the amount of OFD can vary significantly, with OFD ranging from 0% (deformation entirely localized at the primary fault core) to 100% (completely distributed deformation). However, on average, the deformation is more localized at structurally simple sections of the fault and more distributed at structurally complex ones. For example, for the 2013 Mw 7.7 Balochistan earthquake, the percentage of OFD was 26% for structurally simple segments of the rupture, and 56% for structurally complex segments (Zinke et al. 2014). Similarly, the mean amount of OFD for the 1992 Mw 7.2 Landers earthquake was 15% in structurally simple parts of the fault and 50% in structurally complex parts of the fault (Milliner et al. 2015).

The type of near-fault materials also seems to control the amount of OFD. For the Balochistan earthquake, Zinke et al. (2014) found a general tendency for the deformation to be more distributed in younger and thicker sediments, although the trend was weak. Similarly, for the Hector-Mine earthquake, Milliner et al. (2016a) found larger OFD where the rupture propagated through sediments compared with bedrock ( ∼40% vs. 20%). This observation is not systematic however. For the Landers earthquake, Milliner et al. (2015) found similar and large ranges of OFD in areas of sediment–sediment and basement–basement interfaces (mean of and respectively) and lower OFD (mean of 20%) in areas with a strong bimaterial interface (i.e. where the fault juxtaposes bedrock against sediment). Lower amounts of OFD on bimaterial interfaces has also been observed for the 2018 Mw 7.5 Palu earthquake (Socquet et al. 2019) and likely results from enhanced localization along the strong contrast in material rigidity (Andrews and Ben-Zion 1997).

Finally, the amount of OFD can be impacted by the local topography. Weaker sediments tend to dominate the surface geology in flat and low-lying areas (e.g. flood plains), while steep mountainous regions often feature strong bedrock materials; therefore, we may expect to see more OFD where fault ruptures cross flat, sediment-filled areas, and less OFD in mountainous settings. Nevertheless, Zinke et al. (2019) found the opposite (albeit relatively weak) relationship along the Kekerengu segment of the 2016 Mw 7.8 Kaikorua earthquake, with more OFD in the steep mountains compared to the low-slope coastal regions. One explanation proposed by Zinke et al. (2019) is topographic amplification of SH seismic waves, which could in turn enhance the damage within the fault zone, especially where it intersects the surface (e.g. Boore 1972; Ma 2008). Delorme et al. (2020) also noted the influence of topography on surface slip following the 2016 Norcia earthquake (normal), which featured locally enhanced (i.e. very shallow) slip on the primary fault plane along steep sections of the fault. Hollingsworth et al. (2019) also noted that low slope regions generally featured more OFD compared with steep slope regions for the Norcia event, with the steepest slopes featuring locally enhanced slip (shallow slip surplus). Gravitational effects have been postulated as a contributing factor for locally enhancing the slip where the fault breaks the surface (Delorme et al. 2020).

High-resolution optical displacement fields can also be used to infer surface strain tensors (e.g. Scott et al. 2018; Barnhart et al. 2020). The zone of inelastic failure can be readily observed by delineating the areas where the threshold between elastic and inelastic strain is exceeded (a ratio of 0.5% is commonly assumed for surface materials; for example, Scott et al. (2018). In a recent study, Barnhart et al. (2020) inferred the surface-strain tensors of the 2019 Mw 7.1 Ridgecrest, California, earthquake and showed that the zone of inelastic deformation inferred from the shear strain tensor (where the 0.5% threshold was exceeded) was narrow (width of ∼25 m) and did not correlate with the width of off-fault deformation inferred from fault-perpendicular profiles ( ∼63 m). Furthermore, they showed that the width inferred from the strain tensors was consistent along the rupture compared to the OFD width, which was much more variable. They also showed that no correlation exists between the inelastic failure width and the near-fault materials. According to the authors, the inelastic failure width is influenced by the fault maturity (increasing width with maturity) and by the rupture speed (narrower width with increasing rupture speed).

Understanding the degree to which the deformation is localized within the fault zones has several important implications. First, if the deformation is distributed within hundreds of meters from the main fault trace, fault slip rates deduced from measurements of offset geological markers made within tens of meters of the main fault trace could be significantly underestimated (Dolan and Haravitch 2014). If real, this underestimation has several consequences. First, the estimation of seismic hazard using probabilistic approaches incorporating geological fault slip rates could lead to a systematic underestimation of the real hazard. Second, any interpretation of the differences observed between geological slip rate and geodetic loading rate would be biased. The degree to which the deformation is localized can also influence the frequency content and magnitude of strong ground motions during an earthquake (Dolan and Haravitch 2014). Better understanding the link between off-fault deformation and strong ground motion could therefore help to make more accurate predictions of strong ground motions in future earthquakes and thus better assess seismic hazard. Finally, partitioning of slip into damage zones has major implications for how we interpret information recorded in paleoseismic trenches, as offset markers only capture a component of the slip at depth, and may thus lead to an incorrect estimate of past earthquake magnitudes, which also form an integral part of seismic hazard assessment.

High-resolution optical maps of coseismic surface displacements can also be combined with dynamic rupture modeling incorporating off-fault damage to determine the most likely rupture scenario for an earthquake. Indeed, by comparing the pattern of off-fault damage observed in the optical displacement maps with the one simulated by the dynamic rupture modeling, Klinger et al. (2018) succeeded in determining the most likely rupture scenario at the Jordan–Kekerengu–Papatea triple junction during the 2016 Mw 7.8 Kaikoura earthquake.

7.3.4.2. Improvement of fault slip model

Over the last 40 years, fault slip distribution models are routinely inferred each time a significant earthquake occurs. As the estimation of fault slip from surface data is a fundamentally ill-posed problem, slip models for a given earthquake can be very different depending on the type and quantity of data, the inversion method and the assumed fault geometry and Earth model. In the last three decades, the development of space-based geodetic techniques, in particular InSAR, has given access to increasingly complete descriptions of earthquake surface displacements, thereby increasing the robustness of slip models. However, InSAR does not routinely provide near-fault measurements for large surface-rupturing earthquakes. By improving the near-fault coverage, optical correlation data can increase the reliability of fault slip models (e.g. Dominguez et al. 2003; Xu et al. 2016; Gombert et al. 2017; Marchandon et al. 2018b; Delorme et al. 2020; Lauer et al. 2020). For example, by re-estimating the slip model of the 1997 Mw 7.2 Zirkuh earthquake, Marchandon et al. (2018b) have shown that the addition of near-field optical data in the modeling changes the obtained fault geometry, the slip amplitude and the slip distribution significantly compared to the slip model based on InSAR data only.

Furthermore, near-fault optical data have been used to better understand the origin of the shallow slip deficit (SSD) commonly seen in fault slip models for large continental strike-slip earthquakes (Fialko et al. 2005). By reassessing the fault slip model of three large strike-slip earthquakes and incorporating near-field displacements from optical correlation data, Xu et al. (2016) have shown that a significant part of the previously inferred SSD for each event was an artifact due to a lack of near-fault coverage. The impact of the lack of near-fault coverage for SSD has also been shown by Lauer et al. (2020) for the Balochistan earthquake and by Marchandon et al. (2021) using synthetic tests. However, because optical data are typically much noisier compared to InSAR, noise-related model uncertainties should be systematically estimated when using optical data in fault slip inversions in order to avoid any misinterpretations of the inferred slip models, in particular concerning the distribution of slip with depth (Marchandon et al. 2021).

7.3.4.3. Impact of the fault structural properties on the coseismic slip distribution

It has long been noticed that fault complexities can impact rupture propagation during an earthquake (King and Nabelek 1985). This impact is reflected in the coseismic slip distribution of earthquakes and shows that fault complexities such as step-overs or bends correlate with (1) zones of low slip amplitude (e.g. Lasserre et al. 2005) and (2) the start and endpoint of ruptures (e.g. King and Nabelek 1985). By providing the geometry of the surface rupture and the distribution of surface fault slip, optical correlation results have helped to increase the number of examples that illustrate the impact of fault geometry on coseismic slip (e.g. Klinger et al. 2006; Milliner et al. 2015; Vallage et al. 2015). For example, analysis of the surface coseismic slip distributions of the Khuli-Boniabad and Zirkuh earthquakes obtained from optical correlation have shown that the ruptures of these earthquakes are characterized by several segments, with the inter-segments located in zones of large-scale structural complexity of the Dasht-e-Bayaz and Abiz faults (Marchandon et al. 2018a, 2018b).

Furthermore, high-resolution optical measurements allow a more in-depth analysis of the surface slip variability. In particular, the frequency content of surface coseismic slip distributions can be analyzed by computing the power spectral density of surface slip (e.g. Milliner et al. 2016; Bruhat et al. 2020). Power spectral density (PSD) is the measure of the power content of a signal versus frequency and can be used to characterize fractal properties of a signal. By computing the PSD of the surface slip distribution of various earthquakes, it has been shown that the surface slip distribution is fractal (Milliner et al. 2016; Bruhat et al. 2020) (i.e. a power law describes the relationship between the slip amplitude and the spatial wavelength). The fractal nature of surface slip echoes the fractal nature of fault surfaces. Indeed, although faults are typically treated as planes, faults are in reality irregular surfaces at all scales (e.g. Candela et al. 2012). Dynamic rupture simulations have shown that these geometric irregularities, i.e. roughness, induce local stress perturbations that affect the rupture prorogation (e.g. rapid accelerations and decelerations of the rupture front, variability in moment release and super-shear transitions) and induce slip heterogeneities (e.g. Dieterich and Smith 2009). By analyzing the coseismic slip variability of real earthquakes, it could thus be, in principle, possible to investigate the roughness properties of the faults with which they occurred. Bruhat et al. (2020) made a first attempt towards this goal. Using a catalog containing thousands of dynamic rupture simulations on a self-similar fractal fault (H=1) with varying roughness, they investigated the relationship between the fault roughness and the spectral characteristics of the resulting slip distribution. They found that the slip distribution deviates from self-similarity (Hurst coefficient decreases) with increasing fault roughness. They also analyzed the spectral characteristics of various real earthquakes and found that the real slip distributions were significantly more self-affine (lower Hurst coefficient) compared to the simulated slip distributions, thus making it hard to draw any conclusions about the fault roughness with which these real earthquakes occurred.