Data‐Driven Two‐Fault Modeling of the Mw 6.0 2008 Wells, Nevada Earthquake Suggests a Listric Fault Rupture

Structural fault complexity at depth affects seismic hazard, earthquake physics, and regional tectonic behavior, but constraining such complexity is challenging. We present earthquake source models of the February 21, 2008, Mw 6.0 Wells event that occurred in the Basin and Range in the western USA, suggesting the rupture of both the shallow and deep parts of a listric fault. We use a large data set including 150 local seismic waveforms from the USArray combined with high‐quality Interferometric Synthetic Aperture Radar and teleseismic waveforms. Rather than imposing an a priori fault geometry in the source inversions, as is often done in the literature, we use a data‐driven approach whereby all the faulting parameters and number of faults are determined by the data alone. We find a two‐fault normal faulting solution comprising: (i) a shallow (centroid depth ∼4.6 km) sub‐event with Mw 5.3 and fault dip of ∼77°; and (ii) a deeper (centroid depth ∼8.8 km), larger Mw 6.0 sub‐event on a fault with shallower dip angle (∼41°). Our preferred two‐fault model is consistent with aftershocks and with the tectonics of the region. The local USArray waveforms used in the modeling are key to detect the rupture of both shallow and deep parts of the possible listric fault. The lack of such dense and uniform coverage of earthquakes in other regions on Earth may explain why the full seismic rupture of listric faults may have gone undetected in the past. Thus, earthquake slip on whole listric faults may be more common than previously thought.

with depth, where brittle faulting in the upper crust generates less rotation than fault creep in the lower part of the crust (Jackson & McKenzie, 1983). However, it is not clear whether earthquakes rupture whole listric fault structures or only parts of them, notably their deeper, shallower-dipping components. Some seismic studies based on first motion analysis (e.g., Doser & Smith, 1989;Jackson & White, 1989) and on down-dip sub-event seismic waveform inversions (e.g., Braunmiller & Nabelek, 1996) suggested that listric normal faults do not rupture seismically along their full down-dip widths. On the other hand, some seismic and geodetic studies suggested the full rupture of listric faults (e.g., Cakir & Akoglu, 2008;Senturk et al., 2019). However, a listric fault configuration was imposed in those studies as a priori information without verifying whether the data used actually require such a configuration. Knowing whether earthquakes rupture whole or only parts of listric faults has important seismic hazard implications, since simulations suggest that fault listricity could have an important effect on peak ground velocities (Passone & Mai, 2017). As listricity increases, the peak ground velocity in the hanging wall increases and enhances constructive interference effects. The latter can double seismic wave amplitudes compared to those generated by a planar fault (Passone & Mai, 2017).
On 21 February 2008, a Mw 6.0 normal faulting earthquake occurred in Wells, Nevada, in the Basin and Range tectonic province in the western USA . This intraplate event occurred on an unmapped fault segment, filling the gap between the Town Creek section of the Clover Hills fault and the Eastern Snake Mountains fault further North (Henry & Colgan, 2011;Ramelli & dePolo, 2011, Figure 1). The event had a maximum Modified Mercalli Intensity of VIII at Wells (dePolo & Pecoraro, 2011) and caused significant damage in the town of Wells, where seven houses and buildings were destroyed, more than 100 buildings were damaged dePolo & Pecoraro, 2011) and three people were injured (dePolo & Lotspeich, 2011). Previous studies of this event used seismic and InSAR data to determine a normal faulting mechanism with fault strike between 19° and 35°, fault dip in the range 33°-41° and moment magnitude Mw 5.8-6.2 (Dreger et al., 2011;Mendoza & Hartzell, 2009;Nealy et al., 2017). However, existing source models of the event show discrepancies, notably regarding slip at depth (see Table 1 with solutions from previous studies). Mendoza & Hartzell (2009) used local seismic data to find that the earthquake slipped mainly between 6.5 and 10.8 km depth. On the other hand, Dreger et al. (2011) found shallower slip when using InSAR data (at ∼5 km depth). Moreover, Nealy et al. (2017) reported that including InSAR data in the source inversions led to shallower slip (∼7.4 km) than when using a fault geometry inferred from the aftershocks, with slip at ∼11.5 km depth. Figure 1 shows the area around the Wells earthquake together with the aftershocks located by Smith et al. (2011). While previous studies considered that the Wells earthquake occurred on a single planar fault, the aftershock distribution of Smith et al., 2011 suggests a listric fault configuration (see the Movie S1 in the supporting information). In addition, the GCMT and USGS moment tensor solutions of the event estimate non-double-couple components of ∼17%-45% (Table 1), which suggests source complexity.
The Wells earthquake occurred right in the center of the region where the USArray (Transportable Array) was deployed at the time of the event, providing excellent dense seismic coverage. Moreover, it was also captured by multiple tracks of very high-quality InSAR data from the Envisat satellite. Here, we present models of the Wells earthquake built from simultaneous inversions of seismic and InSAR data whereby all faulting parameters and the number of faults are determined by the data alone, without imposing any a priori fault geometry. We show that the data require the seismic rupture of two fault planes suggesting a listric fault and discuss its implications.  Henry and Colgan (2011) and from Ramelli and dePolo (2011). Bottom: Surface projection of the two fault planes determined in this study (red and blue rectangles). The circles are the relocated aftershocks of Smith et al. (2011) color-coded according to their depth. The topography is based on the Shuttle Radar Topography Mission Global one arc second elevation model (NASA JPL, 2013

InSAR Data
The Envisat satellite had a repeat interval of 35 days between radar images acquired on the same orbital track. The ASAR sensor operated in the C-band with a wavelength of 5.6 cm. The ascending image pair (track 220, with the acquisition on November 14, 2007 andApril 02, 2008) and descending image pair (track 399, acquisition on August 13, 2007 andApril 13, 2008) are processed with the software package InSAR Scientific Computing Environment (ISCE, Agram et al., 2013). Two additional InSAR tracks were available from Envisat (tracks 127 and 492), but the data were either very noisy (track 127) or covered a too long time interval (track 492), and hence we did not use them. In any case, the minor differences in incidence angle permitted by these additional tracks would be unlikely to significantly improve our resolution of the deformation field for the Wells earthquake (e.g., Wright et al., 2004). The 1 arc-second Shuttle Radar Topography Mission (SRTM) digital elevation model (Farr et al., 2007;NASA JPL, 2013) is used to remove the topographic phase and to geocode the InSAR data. The phase is unwrapped with Statistical-cost Network-flow Approach for Phase Unwrapping(SNAPHU; Chen and Zebker, 2001). In order to reduce the computational cost of the analysis, the Line Of Sight (LOS) displacements are downsampled using a quadtree decomposition (e.g., Funning et al., 2005;Jonsson, 2002;Weston et al., 2014). Consequently, the downsampled interferogram for the ascending track consists of 633 data points and of 599 data points for the descending track ( Figure 2). In order to remove orbital errors and subsequent long-wavelength phase ramps in the interferograms, three additional parameters for each geodetic data set are solved for in the inversion process and removed from the data: a translation (static offset) in the radar line-of-sight, as well as two linear gradients in the north and east directions. As explained in the Results section below, we found an excellent agreement between the observed InSAR data and the predictions from the co-seismic deformation. We would expect that any unmodelled postseismic deformation that would be present in the InSAR data would lead to systematic residuals in the same place in the ascending and descending interferograms, which do not appear in our results. GPS data coverage is limited in the area, but there is at least one station triplet shown by Hammond et al. (2014) that has a co-seismic offset due to the Wells earthquake, but no nonlinear FRIETSCH ET AL.  Dreger et al. (2011) based on local seismic data; "DS Local BMT": distributed slip solution obtained from local seismic data imposing the fault geometry obtained from a body wave moment tensor inversion (Dreger et al., 2011); "DS EGF": distributed slip solution obtained from an Empirical Green's Function approach (Mendoza & Hartzell, 2009); "DS InSAR, aftershock": distributed slip solution obtained from InSAR data, with the fault geometry from aftershock relocations (Nealy et al.,2017); "DS seismic, aftershock": distributed slip solution from local seismic data (Dreger et al., 2011) based on the fault geometry inferred from aftershocks of Smith et al. (2011); "DS InSAR": distributed slip solution obtained from InSAR data (Nealy et al., 2017). Abbreviations: InSAR, Interferometric Synthetic Aperture Radar. On that basis, we have no reason to expect any detectable postseismic deformation in the InSAR data and thus there are no additional corrections accounting for the contribution of post-seismic deformation to the co-seismic InSAR LOS-displacements.

Seismic Data
We use three-component seismic waveforms bandpass filtered in the range ∼17-180 s, depending on the data type. The filtering parameters are chosen to achieve a good signal-to-noise ratio without substantial loss of information. We select the station with the best signal-to-noise ratio for each 5° azimuthal interval in order to get a balanced azimuthal distribution of the stations ( Figure 3). Table S1 in the supporting information shows all the relevant data selection and processing parameters, notably the epicentral distance range and the filtering parameters used. Figure 4 shows illustrative examples of local and teleseismic waveforms (all waveforms are shown in Figures S3-S5). The epicentral distance range of the teleseismic data (see Figure 3) was selected to avoid sampling the complex lowermost mantle with body waves and also to avoid the interference of multiple orbit wave trains. We choose time windows that fully capture each seismic phase. While the windowing of the body waves is based on phase arrivals estimated by the travel time calculator TauP (Crotwell et al., 1999), the maximum of the envelope function is used to select the time window for the teleseismic surface waves.

Source Inversion Method
We Contribution to L2-norm Misfit 10 -4 component, following the approach of Tape and Tape (2015). The v parameter is calculated directly from the moment tensor eigenvalues and is similar to other measures, such as lune longitude γ or the classical ϵ parameter (Ekström et al., 2012). While source studies in the literature typically impose a priori fault geometries in the inversions, our approach consists of letting the data alone determine the optimal fault parameters and number of faults. Thus, the search for the source parameters for each fault is totally free within reasonable bounds (Table S2). Our method uses the Powell algorithm (Powell, 1964) as its local optimization scheme in an inversion with multiple Monte Carlo restarts in order to find the global minimum of the L 2 -norm misfit. The misfit function m 2 is based on the difference between the theoretical seismograms t S and the observed seismograms d S , as well as the difference between the theoretical geodetic t G and observed geodetic d G displacements with weighting factors for the seismic and geodetic components given by α S and α G , respectively: The number of restarts of the algorithm depends on the complexity of the source model. Having extensively experimented with the number of restarts and iterations, we found that 200 to 1,000 Monte Carlo restarts are enough for a single planar fault model to converge, while a two-fault model needs at least 100,000 iterations.
The InSAR displacements are modeled using elastic dislocation theory for the displacement of a rectangular fault in a homogeneous half-space (Okada, 1985) assuming a constant shear modulus μ = 30 GPa and FRIETSCH ET AL.
10.1029/2020JB020263 6 of 15 corresponding to a pure double-couple source model. Following Weston et al. (2014), the seismic centroid location is fixed to the center of the fault plane obtained from single fault inversions using only InSAR data.
Using the centroid location obtained from the inversions of InSAR data, synthetic seismograms and seismic sensitivity kernels for source parameters are computed with the highly accurate spectral element method, using the SPECFEM3D_GLOBE package (Komatitsch & Tromp, 2002). The global 3-D whole mantle shear wave model SGLOBE-rani (Chang et al., 2015) is employed combined with the CRUST2.0 (Bassin et al., 2000) crustal model. Despite being relatively expensive computationally, given the linear relationship between the moment tensor and the seismic waveforms, sensitivity kernels are pre-computed for the centroid location obtained from the InSAR data and stored before being used in the source inversions.
Model uncertainties are estimated using a heuristic approach, by defining error bounds based on the source models that fit the data within 20% of the minimum misfit value. This misfit threshold value was chosen such that the corresponding models all fit the data reasonably well. First, the different seismic and geodetic data sets are inverted separately to explore their sensitivity to the different source parameters. Then, the weighting factors (α S and α G ) in the joint data inversions are adjusted to give a similar importance to each data set (see Frietsch   The multiple fault inversion scheme successively adds subfaults to the modeling as required by the data. Every source inversion starts by considering a single fault. The number of subfaults is then increased one by one until the improvement in data misfit is smaller than 5% or if the moment magnitude of the additional sub-event becomes smaller than the previous one by one magnitude unit or more. The 5% threshold in the data misfit criteria is consistent with previous seismic studies that used similar thresholds taking data noise into consideration (e.g., Debayle & Ricard, 2012;Parisi & Ferreira, 2016). After performing extensive inversions, we found that once three-fault models were considered, they led to marginal improvements in the data fit (smaller than 5%) compared to two-fault models. Hence, in this study we take the conservative approach of focusing on simpler one-and twofault models.

Results
Table 2 and Figure 1 show the results from single fault and two-fault source inversions of the Wells earthquake. As expected, all the solutions show predominantly normal faulting mechanisms with moment magnitude Mw 6.0. The two-fault solution from the joint data inversion shows a change in fault dip angle from ∼77° at shallow depth (∼5 km) to ∼41° at greater depth (∼9 km). It also shows that the earthquake started with a deeper, larger magnitude (Mw ∼ 6.0) sub-event and extended upwards, with the second shallower, FRIETSCH ET AL. Single fault solutions are represented by a blue beach ball, while two-fault solutions are shown with red beach balls. The fault parameters shown are the fault strike, dip, rake, the double-couple Percentage ("DC"), the CLVD v component (Tape & Tape, 2015), the time lag to the GCMT time, the seismic moment, the moment magnitude, the L 2 -norm misfit (Equation 1; see details in the main text) and the corresponding stereonet plot. We compute the non-double-couple (ϵ) percentage of the solutions from the eigenvalues λ i of the moment tensor using the same approach as e.g. in the GCMT catalog: ϵ = 100 × min(|λ i |)/max(|λ i |) (Ekström et al., 2012), with the double-couple component being DC = 100 − ϵ. We report solutions obtained with various subsets of data whereby "Local" refers to local seismic data (with Epicentral distance smaller than 10°), "Tele" refers to teleseismic data (with Epicentral distance between 30° and 140°), "InSAR" refers to the use of InSAR data alone and "Joint" refers to joint data inversions combining the local, teleseismic and InSAR data. Abbreviations: InSAR, Interferometric Synthetic Aperture Radar.  (2011), which suggested an earthquake nucleation deeper than the event's aftershocks (Figure 1, Movie S1). Apart from the teleseismic source inversions, the differences between the retrieved fault strike values for the two sub-events do not exceed ∼10°. This is within the expected error range for fault strike estimations from source inversions (e.g., Weston et al., 2011Weston et al., , 2012, suggesting that the two sub-events ruptured the same curved fault. Our source models provide a very good fit to the data, with low waveform misfits (Figures 2, 4 and S3-S5 in the supporting information). The final part of the waveforms occasionally show slightly poorer data fits, which is probably due to inaccuracies in the crustal model used in the forward modeling. When going from a single fault model to a two-fault configuration, the local seismic data show the most substantial improvement in data fit of ∼15% (Table 2), well above the 5% data misfit threshold explained in Section 3. When visually comparing the waveform fits of the single and two-fault models it is not very easy to assess the improvements in data fit mostly due to the large magnitude difference between the two sub-events, which leads to the signal of the smaller magnitude sub-fault being partly covered by the signal of the larger magnitude sub-fault. Nevertheless, we can see clear misfit improvements for example for the vertical components of stations L10A, L11A, and K10A in Figure 4. This is consistent with synthetic inversion tests that showed that InSAR and teleseismic data (i.e., without local data) cannot distinguish two normal sub-faults with down-dip segmentation where one of the sub-events is an order of magnitude smaller than the other (Frietsch et al., 2019). Indeed, the two-fault solutions from teleseismic and InSAR inversions presented in Table 2 only show some slight improvements in data fit compared to the one-fault solutions and show some unlikely features. This includes large differences in fault strike between the two sub-events (∼127°) for the teleseismic inversions and a difference of ∼40° in rake for the InSAR-only solutions. In addition, one of the sub-events obtained from the inversion of teleseismic data alone has a strong, non-double-couple component. This is unlikely since it would suggest non-shear deformation, and it must instead be due to instabilities in the inversions as found, for example, by Frietsch et al. (2018) in inversions of other moderate magnitude events. Moreover, the InSAR data inversion predicts right lateral oblique slip for the upper section and left-lateral oblique slip for the lower section of the fault, which is also improbable. This illustrates the difficulties of the InSAR and teleseismic data in resolving the second, lower magnitude sub-event. Although it may sound counter-intuitive that the InSAR data are not able to discriminate the shallow sub-event, Frietsch et al. (2019) showed that this is due to the event's normal faulting geometry, whose displacements are hard to constrain with InSAR when the event's magnitude is relatively low. This would not be a problem if the earthquake had e.g. a strike-slip faulting mechanism, as the fault segments would be spatially separated in map view, rather than superposed.
In order to test the significance of the differences in data misfit between the single and two-fault models obtained from the joint data inversions, we performed a F-test by computing: where χ 2 is the misfit function, n is the number of data points used in the inversions, and p and r are the numbers of model parameters for the two models which are compared (p > r). Under the null hypothesis that the two-fault model does not provide a significantly better fit than the single fault model, the null hypothesis is rejected if F is greater than the critical value of the F-distribution for a given desired false-rejection probability P (Bevington & Robinson, 2002). Since errors in the data and in the modeling used in our inversions are essentially unknown, here, we use the misfit function χ 2 based on the L 2 -norm misfit in Table 2 following, for example, Forsyth (1975) and Ferreira et al. (2019). We find that the probability that the improved fit of the two-fault model compared to the single fault model occurs by chance is less than 1% (i.e., the two models differ at the P = 0.01 level). Thus, the F-tests support the significance of the improved data fit of the two-fault model. Moreover, we also tested our source model against independent data, that is, data not used in the inversions, such as , for example, station ELK, which is a near nodal station at an epicentral distance of ∼0.5°. We found that when considering the source parameter uncertainties, the corresponding synthetic waveforms fit ELK's data well (see Figure S12). Figures 5 and 6 show parameter trade-off plots for the 40 best-fitting solutions from our joint source inversions and Table S3 shows the estimated uncertainties from our inversions. As expected, the source parameters are better constrained for the larger magnitude sub-event ( Figure 5) than for the lower mag-FRIETSCH ET AL.
10.1029/2020JB020263 9 of 15 nitude sub-event ( Figure 6), which is also consistent with the likely artificial lower double-couple component of this latter event obtained from the inversion of the local seismic data alone ( Table 2). The trade-off plots from inversions based on teleseismic or InSAR data individually (Figures S8-S11) show less tightly clustered solutions than when using local seismic data and/or the joint data sets (Figures 5, 6, and S6-S7). Moreover, the uncertainties obtained from the individual data inversions are larger than from the joint data inversions. This highlights the limited power of teleseismic and InSAR data at resolving the two sub-events and the importance of combining complementary data sets. Local seismic data ( Figures S6-S7) have a stronger sensitivity to the time shift than teleseismic data ( Figures S8-S9). In addition, the inversions based solely on teleseismic data lead to well-known trade-offs between the seismic moment and the fault dip as well as between the seismic moment and the fault rake. The good sensitivity to the location is a strength of the InSAR data set; thus, the events' latitude and longitude are well constrained (Figures S10-S11). Finally, we note that as part of our source inversion tests we also performed source inversions imposing a pure double couple constraint and found that the solutions from the joint data inversions were very similar to those reported in this study. While both previous studies of Smith et al. (2011) andNealy et al. (2017) mention the possibility that the 2008 Wells earthquake may have occurred on a listric fault, neither study tested this possibility. The twofault model we present here shows slip both at shallow (3.2 -6.0 km depth) and greater depths (6.8-10.7 km depth), resolving the depth discrepancy of previous studies e.g., between Dreger et al. (2011) andNealy et al. (2017). Amongst other techniques, these studies used distributed slip inversions, which can be highly non-unique (e.g., López-Comino et al., 2015). The multiple fault inversion used in this study is an attractive, data-driven and simpler alternative to robustly constrain earthquake source parameters. When considering the two-fault solution obtained in our study, the surface projection of the shallow sub-fault fills a gap between the Town Creek Flat Fault (TCFF, Figure 1) and the Eastern Snake Mountain Fault (ESMF), while the surface projection of the deep sub-fault is further away into the Snake mountains, which is less geologically plausible. Moreover, Figure 1 shows that, as expected, the aftershocks surround the two imaged sub-faults, corresponding to regions of low slip (e.g., Woessner et al., 2006). This adds further confidence to our twofault solution. We note that for internal consistency of our inversion algorithm, future work will focus on the use of a layered half-space in the modeling of the geodetic data based on the same local layered structure in the source region as in the 3-D Earth models used in the seismic modeling. Nevertheless, this should not change the main results of the joint source inversions because previous studies showed that using layered models versus homogeneous half-spaces in the geodetic modeling affects mainly the source depth and by no more than about 30% (e.g., Hearn & Bürgmann, 2005;Lohman, 2002;Marshall et al., 1991;Savage, 1987).

Discussion
Listric faults are common in the Basin and Range and have been imaged with seismic reflection (Anderson et al., 1983). However, up to now it was not entirely clear whether a whole listric fault could rupture (e.g., Braunmiller & Nabelek, 1996;Jackson & White, 1989). For example, Stein and Barrientos (1985) and Ward and Barrientos (1986)  earthquake, which also occurred in the Basin and Range, ∼340 km away from Wells. Both studies found that a planar fault was in better agreement with their leveling data set. Other studies focusing on other regions also reported the rupture of entire listric faults, notably in Turkey (e.g., Cakir & Akoglu, 2008;Senturk et al., 2019). However, these studies imposed a listric fault configuration a priori in their source inversions and did not verify whether the data actually require such configuration compared to simpler fault systems.
In contrast, in this study we let the fault geometry and number of faults be freely determined by the data FRIETSCH ET AL. We note that the discrepancy of ∼4 km between the location of the mainshock hypocenter and the modeled fault is within depth errors in source inversions (∼10 km when using teleseismic data and ∼6 km when using regional data; e.g., Weston et al., 2012). Moreover, the observed differences between the hypocenter and our deeper fault plane may also be due to the fact that we model the source's centroid, not the hypocenter, and to the use of a homogeneous elastic half-space in the InSAR data modeling (see the Discussion section). We also note that the best-fitting one-fault solution is very similar to that of the largest magnitude, deeper sub-event in the two-fault solution shown in blue. InSAR, Interferometric Synthetic Aperture Radar.
alone. Hence, to the best of our knowledge, this is the first time that free source inversions demonstrate that geophysical data require a rupture of both the shallow and deep parts of a two-fault system suggesting a listric normal fault.
The long duration of the local seismic signals observed in this study (Figure 3) was also noticed by Biasi and Smith (2011). While this may be partly caused by thick sediments in the Town Creek Flats (Biasi & Smith, 2011), we find that this is linked to a lower magnitude, shallower sub-event rupturing ∼15 s after the deeper, main sub-event. We suggest that the second sub-event is due to static triggering and also argue that an alternative explanation, such as unmodelled structural effects, is highly unlikely, as it would be very unusual for structural effects to lead to the exact same azimuthal dependency of the waveforms as in the radiation pattern of the second sub-event. Although the shallower sub-event shows stronger parameter trade-offs (Figures 5 and 6) and larger uncertainties (Table S3) than the deeper, larger magnitude sub-event, their parameter uncertainties obtained from the joint data inversions are not excessive and are comparable to those reported for global earthquakes (Weston et al., 2011(Weston et al., , 2012. The two-fault solution suggesting a listric fault found in this study for the 2008 Wells earthquake is a good example of source complexity at depth, for which there were no previous clues from surface observations. We show that both deep and shallow parts of possible listric faults can rupture in a single earthquake, which could have important implications for seismic hazard assessment in the region (Passone & Mai, 2017).
We note that our solutions consist of a relatively simple two-fault system and that ideally many planar subfaults could be considered in our inversions to define a realistic listric, curved fault. However, as mentioned previously, even when considering three-fault configurations, the data fit only improved marginally, showing that current data sets cannot resolve such detailed faulting configurations. Future higher frequency analysis such as, for example using back projection (e.g., Ishii et al., 2005) could potentially bring further insights into the event's source process. Nevertheless, the two-fault solution found by our data-driven approach along with comparisons with aftershocks and the region's tectonics provide strong support for the rupture of both the shallow and deep parts of a likely listric fault in the region. More generally, this study forms a useful basis to study the complexity of other moderate magnitude events, which remains relatively under-studied.

Conclusions
In this study, we used a unique data set of local seismic waveforms from the USArray combined with high-quality InSAR data and teleseismic waveforms to investigate the complexity of the Mw 6.0, February 21, 2008 Wells, Nevada earthquake. The preferred two-fault solution obtained suggests a listric fault comprising a shallow (centroid depth ∼4.6 km), sub-event with Mw 5.3 and fault dip of ∼77°, and a deeper (centroid depth ∼8.8 km), larger Mw 6.0 sub-event on a subfault with shallower fault dip of 41°. The rupture initiated on the deeper subfault and propagated updip, with an interval of ∼15 s between the two sub-events. This is in agreement with the aftershock distribution of Smith et al. (2011). Listric faults are common in areas of extension such as the Basin and Range and thus the two-fault solution obtained in this study is consistent with the tectonics of the region. We show that ruptures can occur in both the shallow and deep parts of possible listric faults, which has key implications for seismic hazard analyses (Passone & Mai, 2017). The local USArray waveforms used in this study were key to unraveling that both shallow and deep parts of the possible listric fault slipped in the earthquake. The absence of such excellent, dense network coverage in other regions may explain why the rupture of whole listric faults has not been detected previously. They may be more common than previously thought.

Data Availability Statement
The seismic data used in this study are freely available in the IRIS repository (http://ds.iris.edu/ds/nodes/ dmc/data/) and the InSAR Envisat ASAR data are freely provided by the European Space Agency (ESA) and by the WInSAR/UNAVCO SAR archive: https://earth.esa.int/web/guest/missions/esa-operational-eo-missions/envisat/instruments/asar; https://www.unavco.org/data/sar/sar.html.