Shallow destructive earthquakes

Depths of earthquake occurrence and large slip distribution are critical for seismic hazard assessment. Numerous examples show that earthquakes with similar magnitudes, however, can result in significantly different ground shaking and damage. One of the critical factors is that whether the large slip was generated near the ground surface. In this article, we reviewed two aspects that are important on this regard, shallow slip deficit and nucleation depth. Understanding how shallow future earthquakes may nucleate in particular regions, such as shale gas fields, is critical for hazard assessment. Whether or not a strong earthquake may slip significantly at shallow depths (less than 3 km) plays crucial rules in seismic hazard preparation and should be further investigated by integrating high-resolution fault zone observations, dynamic rupture simulation, and fault zone properties. Moreover, precisely resolving shallow depth and slip distribution of earthquakes demands InSAR and/or other image data that can better capture the near-fault deformation to constrain the source parameters of earthquakes.


Introduction
Earthquakes occur in a variety of depths, from shallow crust to the base of the mantle transition zone. Destructive earthquakes, however, are mostly shallow, concentrating in the upper crust of continents and the shallow portions on the megathrust. Apparently, focal and centroid depths of earthquakes play critical roles in the intensity of ground shaking and thus directly determine the destructing power for earthquakes with similar magnitudes. Understanding what factors determined the depths of large slip is thus critical for seismic hazard assessment. Here we review two aspects from observations: shallow slip of large earthquakes and nucleation depths of moderate-sized earthquakes, and then highlight the significance of identifying and predicting such parameters for future hazard preparation of continental earthquakes.

Shallow slip deficit
For shallow earthquakes with magnitudes larger than 6, they are usually very destructive in producing significant damages to the infrastructure and society. However, a number of large continental earthquakes are observed to host the maximum coseismic slip at depths of 3−6 km ( Figure 1) (Fialko et al., 2005;Xu XH et al., 2016), leaving a shallow zone of slip deficit. Such shallow slip deficit (SSD) found in large earthquakes is independent of whether the ruptures had extended to the ground. For instance, the 26 December 2003 M W 6.5 Bam earthquake in southeastern Iran ruptured the fault at depths of 3−7 km with the maximum slip of ~2 m (Fialko et al., 2005). The rupture of this event did not break the surface. While the 1992 M W 7.3 Landers earthquake in southern California broke the ground but the maximum slip occurred in the middle of the seismogenic layer at 4−6 km ( Figure 1).

M W
Indeed, the centroid depths of global earthquakes larger than 5.5 show a cluster of ~15 km, despite their different tectonic settings and faulting mechanisms ( Figure 2a). Note that the catalog of Global Centroid Moment Tensor (GCMT) may lead to a bias in deeper focal depths of earthquakes, as they fixed the minimum depth of 12 km, probably the reason of the large peak at 15 km ( Figure 2a). However, the lack of relatively large earthquakes (e.g. >5) at depths shallower than 3 km is also reflected from other catalogs, in which moment tensor solutions of smaller magnitude earthquakes had been carefully resolved with regional velocity models using the Computer Programs in Seismology (Herrmann, 2013). For earthquakes in North America, most of them are distributed below ~7 km in the upper crust ( Figure 2b). In Yunnan, southwest China, centroid depths of earthquakes exhibit concentrations around 10 km (Figure 2c) (Xu Y et al., 2020). Although sporadic shallow earthquakes (e.g. less than 3 km) were reported in both catalogs, there is no doubt that accumulated coseismic slip occurs significantly below 5 km (Figure 2b and 2c).
Efforts have been made to quantitatively understand the source mechanism of the SSD. It has been suggested that the SSD was possibly compensated by shallow creep or afterslip on faults (Scholz, 1998). Although steady-state shallow creep was observed in certain fault segments, e.g. the southern section of the San Andreas Fault (Lyons and Sandwell, 2003) and a few faults in central California (Schulz et al., 1982), most faults have not been observed with shallow creep during interseismic period. Furthermore, geodetic observations indicate that afterslip was too low to account for the SSD generated in a few earthquakes (e.g. Fialko, 2004;Fialko et al., 2005). Therefore, neither shallow creep nor afterslip could be the source of the SSD.
Alternatively, the SSD has been attributed to shallow compliant fault zones that accommodate a distributed inelastic deformation within the uppermost few kilometres of the crust during the interseismic period (Fialko et al., 2005). The compliant fault zones may yield plastically during coseismic ruptures, leading to distributed off-fault inelastic deformation and thus compensating for the SSD. Such mechanism has been investigated by dynamic rupture simulations in both 2D (Kaneko and Fialko, 2011) and 3D models (Roten et al., 2017), in which the predicted degrees of the SSD were close to observations. Such  Fialko et al., (2005) and Kaneko and Fialko (2011). From Kaneko and Fialko (2011).  compliant fault zones have been documented by geodetic and seismic data, in particular along the fault zones ruptured during the 1992 Landers (Li YG et al., 1994;Peng ZG et al., 2003;Li HY et al., 2007) and 1999 Hector Mine earthquakes (Li YG et al., 2002), which both produced significant SSD (Fialko et al., 2005;Xu XH et al., 2016). Moreover, the compliant fault zones are seismically characterized as low velocity zones that had been found along a number of continental faults (Ben-Zion and Sammis, 2003;Cochran et al., 2009;Yang HF et al., 2011;Ma C et al., 2020;Tian FF et al., 2020;Yang HF et al., 2010;Yang HF et al., 2014;Yang HF, 2015;Yang HF et al., 2020a). However, the low velocity zones (LVZ) at shallow depths associated with seismogenic faults have been suggested to promote rupture extents in 3D dynamic rupture simulations, although the ruptures did not extend into the LVZ (Weng HH et al., 2016). If the rupture extends to the surface, slip at shallow depths becomes larger on the shallow portion of the fault that was surrounded by the compliant materials with smaller rigidity (Figure 3). Note that such conclusions were drawn based on the assumption of elastic response of the LVZ during coseismic ruptures (Weng HH et al., 2016;Chen X and Yang HF, 2020). When considering inelastic off-fault response in a 2D model, final slip on the fault was also promoted by the presence of the LVZ and the slip increased with the LVZ width (Duan BC, 2008). Therefore, it demands further investigations how the shallow compliant fault zone materials impact on rupture propagation and slip distribution. In particular, how to quantify the degree of plastic yielding of off-fault materials at shallow depths and how they may alter the slip distribution near surface are critical to precisely predict shallow slip distribution in future earthquakes along continental faults.
Furthermore, precisely resolving shallow slip behaviors requires extensive geodetic observations close to the surface rupture. In addition, the estimate of SSD may be discrepant among models with different geodetic datasets or other critical model inputs such as the fault geometry. For instance, the SSD of the 1992 Landers M W 7.3 earthquake was originally reported to be ~ 50% in the slip model ( Figure 1) that was inversed from the GPS data and InSAR data (Fialko, 2004). While a recent study (Xu XH et al., 2016), which revisits this event by further incorporating the optical imagery data and the SAR azimuth offsets, indicates a much smaller degree (20%) of the SSD (Figure 4a). Similarly, smaller SSDs are inferred for the Hector Mine M W 7.1 and the El Mayer-Cucapah M W 7.2 earthquakes after adding more geodetic data and more detailed processing into the inversions (Figure 4a) (Xu XH et al., 2016). Advances in the near-fault geodetic measurements are critical in precisely constraining the SSD.
In stark contrast, some earthquakes apparently do not exhibit such SSD and may significantly break the ground surface. For instance, the 3 August 2014 M W 6.2 Ludian earthquake ruptured the conjugate Zhaotong-Ludian fault system (Zhang GW et al., 2014), with the maximum slip of 1.6 m near the surface on the NNW-SSE oriented fault plane (Niu YF et al., 2020). Although slip distribution and the maximum slip differ in kinematic rupture models that were constrained from different data (Liu CL et al., 2014;Zhang Y et al., 2015), it is consistent that the maximum slip occurred near surface (Figure 4b), one of the reasons leading to severe damages during the 2014 Ludian earthquake. Furthermore, the 14 April 2010 M W 6.9 Yushu earthquake produced a maximum slip of 2.0 m at 3-km depth, with > 1.5 m slip near the ground surface on the Jiegu segment . Although details of slip distribution were slightly different, the characteristic of significant slip (~ 1 m) near the ground was consistent among different source models of the Yushu earthquake that were constrained from satellite images (Figure 4b) (Jiang GY et al., 2013;Tobita et al., 2011;Wang XW et al., 2014). Whether coseismic ruptures extend to the ground surface is critical for seismic hazard. Therefore, investigating what determined such behavior of earthquake ruptures is demanded with high priority.

Nucleation depths of moderate-sized earthquakes
In addition to the SSD produced by large continental earthquakes, one important factor is how shallow earthquakes nucleate, because even a moderate-sized earthquake at very shallow depth can be quite damaging. Numerous reports show that large earthquakes do not nucleate at shallow depths. Based on spontaneous rupture models, Das and Scholz (1983) have suggested that ruptures initiating at greater depths, an environment that is anticipated to have larger stress drop, can propagate over the entire fault plane. In contrast, the ruptures starting from shallow depths where frictional strength and stress drop are smaller are inhibited due to the insufficient strain energy (Das and Scholz, 1983). Although frictional strength may not necessarily increase with depth from dynamic rupture modeling results with near-field constraints (Weng HH and Yang HF, 2018;Yao SL and Yang HF, 2020), it is indeed rare to observe very shallow focal depths of large earthquakes.
It is understood from low-velocity rock sliding experiments that rock samples under conditions representing depths shallower than 2−3 km generally exhibit velocity strengthening (VS), i.e. frictional resistance increases with slip rate, unfavorable for fault acceleration (Blanpied et al., 1998;He CR et al., 2007). Therefore, these areas are anticipated with aseismic slip while in contrast earthquakes should only occur in velocity weakening (VW) conditions, which correspond to upper crust (Scholz, 1998). It is, however, recognized that the VS zones may be broken due to a large rupture nucleated in nearby VW regions (Kaneko et al., 2010;Yang HF et al., 2012). But ruptures would not nucleate within the VS region by themselves.
However, reports have shown that moderate-to-large magnitude earthquakes do occur at very shallow depths, some of which became quite destructive. By analysis of depth phases such as sP and Rg of foreshocks and aftershocks of the 1968 M S 6.8 Meckering, west Australia, earthquake, Langston (1987) suggested that the foreshocks were located mostly at less than 2 km depth and most aftershocks occurred within 1 km of the surface, consistent with a previous study of the mainshock based on teleseismic waveforms. The locations of foreshocks and aftershocks implied that the mainshock rupture initiated at depth of 2 km and propagated downward to 7 km, which was inferred from the maximum depth of aftershocks (Langston, 1987). Such shallow nucleation depth and downward propagation of rupture marked the 1968 M S 6.8 Meckering event a very rare earthquake, plausibly a result of high strength in the cold shield crust.
Such shallow earthquakes were not only reported in intraplate shield, but also along major block boundaries. In the northeastern end of the rupture zone of the 12 May 2008 M W 7.9 Wenchuan earthquake, a M S 5.7 aftershock occurred at Qingchuan on 24 July 2008. Based on waveform modeling of sPL phases of its aftershocks, they were located at depth less than 3 km (Luo Y et al., 2010). In addition, the well-developed Rayleigh wave of the M S 5.7 earthquake recorded at a station within 15 km of epicentral distance implied a focal depth of no more than 3 km (Luo Y et al., 2010). Such shallow earthquakes were attributed to stress increase from the deep coseismic slip during the 2008 M W 7.9 Wenchuan earthquake. More recently, on 11 August 2016, a M W 4.1 thrust earthquake occurred in the Eastern Sichuan Basin and led to three injuries, in a region with historically low seismicity. Seismic waveform modeling constrained the centroid depth of this event to be 1 km, which was consistent with ground deformation that was captured by InSAR data. This event was interpreted to be induced by removal of 10 m thick surface rock layer during the construction (Qian YY et al., 2019).
Moreover, a shallow earthquake with M W 4.3 occurred on 25 February 2019, in the shale gas field of Rongxian, Sichuan. This earthquake was constrained by InSAR and seismic data to have a centroid depth of 1 km on the Molin fault ( Figure 5), one of the major reasons for this earthquake to be very destructive (Yang HF et al., 2020b). Two fatalities and 12 injuries were caused by the M W 4.3 Rongxian earthquake, despite its relatively small magnitude. Analysis of spatial and temporal correlation between the earthquake and nearby hydraulic fracturing wells suggests a possible causal relationship. But it remains enigmatic why such a shallow fault was activated to slip coseismically (Yang HF et al., 2020b;Wang MM et al., 2020). The host rocks of the Molin fault are marlstone and limestone. According to results of laboratory frictional experiments on sedimentary rock samples (mudstone, sandstone, limestone) from the Sichuan Basin, all samples under conditions representing shallow depths (< 2 km) exhibit VS behavior (Verberne et al., 2010), and thus are anticipated to slip aseismically. The shallow 2019 February M W 4.3 earthquake apparently challenges such traditional view and raises a significant question of how to evaluate potential seismic hazard and risks for hydraulic fracturing near very shallow faults. It is also demanded to further investigate the probability of coseismic slip on very shallow faults.
Accompanying the occurrence of shallow damaging earthquakes, challenges arise in resolving earthquake source parameters, especially in focal depths, which are important in elucidating mechanisms of why such earthquakes occurred. For instance, the focal depth of the 2019 February M W 4.3 Rongxian earthquake was suggested to be ~2 km (Lei XL et al., 2020;Yi GX et al., 2020), slightly greater than the results in Yang HF et al., (2020b). In general, such small difference would not matter much for earthquakes occurring below 5 km. However, the precise depth estimation becomes critical for shallow earthquakes, e.g. the 2019 Rongxian earthquake, because the Molin fault was found terminating at a depth of ~1.5 km (Figure 5b) (Wang MM et al., 2020). Focal depths with large uncertainties may lead to problematic interpretation of source faults in the region.
In addition, the focal depths of the two largest earthquakes in the Weiyuan Shale Gas Field (WSGF) are in debate, probably due to the limited station coverage. On September 8, 2019, a M W 5.0 earthquake occurred in the east cluster to Weiyuan. Using the crust 1.0 velocity model, Wang MM et al., (2020) constrained the focal depth at 5 km, close to the result (4.5 km) using another velocity model that was derived from relocated local  earthquakes (Yi GX et al., 2020). However, two other studies resolved its focal depth at 2.5 km (Lei XL et al., 2020) and 2.9 km (Sheng et al., 2020), respectively, which were slightly shallower than the injection depth (3.0 km). Similarly, the focal depths of the 18 December 2019 M W 4.9 earthquake differed from 4 km (Yi GX et al., 2020), 3 km (Sheng et al., 2020), to 2.5 km (Lei XL et al., 2020). Resolving earthquake source parameters in high resolution is necessary to infer whether and how these earthquakes were induced by fracking activities in the WSGF.

Discussion and future directions
Although most large earthquakes nucleate in the middle or near the base of the seismogenic zone and leave the SSD, there are exceptional cases of causing large slip at very shallow depth and thus leading to severe ground shaking and damage, as discussed in above examples of the 2010 M W 6.1 Yushu and 2014 M W 6.9 Ludian earthquakes. Evaluating whether or not large shallow slip may occur along continental faults is thus necessary for future seismic hazard assessment. Indeed, earthquake rupture propagation is a complex nonlinear process that is governed by a number of factors. Heterogeneities of fault zone properties, including fault zone materials (e.g. Duan BC, 2008;Harris and Day, 1997;Weng HH et al., 2016), fault geometry (e.g. Yang HF et al., 2013;Yu HY et al., 2018), stress distribution and frictional properties (e.g. Weng HH et al., 2015;Weng HH and Yang HF, 2018;Yao SL and Yang HF, 2020), play critical roles in slip distribution. Although stress condition on seismogenic faults remains elusive, various lines of evidence indicate that stress distribution on faults is heterogeneous. As a result, ruptures initiating at different positions on the fault plane may result in distinctly different rupture extent and slip distribution , making it more challenging to predict potential shallow slip in future earthquakes. Integrating seismological observations on fault zones, frictional properties of rock samples representing shallow faults, and dynamic mechanism of rupture propagation is desired to better tackle such a challenging problem.
On the other hand, accurately resolving focal depths of shallow earthquakes is not only critical for elucidating tectonic setting and potential inducing mechanisms of faulting, but also holds important implications for earthquake emergency response. As shown in the case of 8 September 2019 M W 5.0 earthquake in Weiyuan, it is controversial to resolve the centroid depth using conventional seismological techniques, whose resolutions are subject to station coverage, frequency band of waveforms, and imperfect velocity models. Indeed, the residual difference in waveform-based seismological solutions such as from the gCAP method (Zhu LP and Ben-Zion, 2013) may not be that distinct for shallow earthquakes. For instance, the misfit values of focal depths at 1 km and 2 km for the 2019 Febuary M W 4.3 Rongxian earthquake are quite close (Yang HF et al., 2020b). Even results of forward waveform modeling show subtle changes in synthetic waveforms at most stations, including the closest one that was 13.5 km in distance ( Figure 6). Therefore, resolving accurate depths of shallow earthquakes demands additional constraints.
Detailed analysis of depth phases such as sPL and sP is very helpful to constrain the earthquake depth (e.g. Langston, 1987;Luo Y et al., 2010), but is usually subject to station coverage and data quality. In contrast, InSAR data has been proven very effective to constrain the depths of shallow moderate-magnitude earthquakes (e.g. Qian YY et al., 2019;Yang HF et al., 2020b). Since the application on the 1992 Landers earthquake (Massonnet et al., 1994), InSAR data have been extensively used to constrain co-, inter-, and post-seismic deformation of numerous large continental earthquakes. In addition to focal depths, it has also been well noted that constraining slip distribution jointly from InSAR and seismic data can significantly improve the resolution and robustness in kinematic source models (Yue H et al., 2020). Incorporating InSAR data to constrain source parameters of earthquakes is very helpful to precisely resolve shallow depth and slip distribution of earthquakes.
Although moderate earthquakes usually produced little damage, they may become destructive if they occur at shallow depths, as evidenced by recent earthquakes in the Sichuan Basin (Qian YY et al., 2019;Yang HF et al., 2020b;Sheng et al., 2020). Even though M4 earthquakes may not occur that often, earthquakes with smaller magnitudes (e.g. 3) may be frequent and can cause damage to underground infrastructure, such as underground gas storage Jiang GY et al., 2020), CO 2 sequestration (Zoback and Gorelick, 2012), and horizontal wells to conduct hydraulic fracturing during shale gas development. In addition to small earthquakes, aseismic slip may also lead to damage to the infrastructure. For instance, casing distortion in the Shale Gas Fields at Sichuan Basin is severe (Xie J, 2018), and it is unclear whether the distortion was caused by earthquakes or aseismic slip. Understanding the fault slip at shallow depths in a full spectrum spanning from long-term and episodic creep, slow slip events, and earthquakes is necessary not only for advancing our understanding of faulting mechanics and earthquake hazard preparation, but also holds critical keys to prevent underground infrastructure from potential failure and damage. 3 Rongxian earthquake at 1 (red) and 2 km (blue) depths, respectively. Station names, waveform cross correlation coefficient, and epicentral distance (in km) are marked for each station.