Complex Rupture of an Immature Fault Zone: A Simultaneous Kinematic Model of the 2019 Ridgecrest, CA Earthquakes

The 4 July 2019 Mw6.4 and subsequent 6 July 2019 Mw7.1 Ridgecrest sequence earthquakes (CA, USA) ruptured orthogonal fault planes in a low slip rate (1 mm/year) dextral fault zone in the area linking the Eastern California Shear Zone and Walker Lane. This region accommodates nearly one fourth of plate boundary motion and has been proposed to be an incipient transform fault system that could eventually become the main tectonic boundary, replacing the San Andreas Fault. We investigate the rupture process of these events using a novel simultaneous kinematic slip method with joint inversion of high‐rate GNSS, strong motion, GNSS static offset, and Interferometric Synthetic Aperture Radar data. We model the Coulomb stress change to evaluate how the Mw6.4 earthquake may have affected the subsequent Mw7.1 event. Our findings suggest complex interactions between several fault structures, including dynamic and static triggering, and provide important context for regional seismic source characterization and hazard models.


Introduction
The Ridgecrest sequence earthquakes ruptured the Little Lake Fault Zone, a low slip-rate (1 mm/year) region dominated by northwest striking, dextral strike-slip faulting linking the Eastern California Shear Zone (ECSZ) and Walker Lane (Figure 1). The area is located on the western edge of the Basin and Range Province and bounded by Garlock Fault and Mojave Block to the south. To the west, the San Andreas Fault (SAF) is the main tectonic boundary between the Pacific and North American Plates. The SAF has a long history of stepping inland (Atwater & Stock, 1998), transferring the tectonic motion eastward over time. The Walker Lane/ECSZ currently accommodates 9-23% of the total relative plate motion between the Pacific and North American Plates (Dokka & Travis, 1990). It has been proposed that in the future, this diffuse region may coalesce into a major transform boundary and assume the bulk of the tectonic motion presently accommodated by the SAF (Faulds et al., 2005). The study region is therefore an example of an immature fault system in an early stage of development. Hence, the interaction between the 2019 4 July M w 6.4 and 6 July M w 7.1 earthquakes provides important context for the development of major tectonic boundaries and the future regional seismic hazard.
The 4 July M w 6.4 earthquake occurred on northeast trending faults, roughly 20 km north of the Garlock Fault. The surface trace suggests that rupture crossed a 1.3-km step over directly to the south of the epicenter (Figure 1). The subsequent 6 July M w 7.1 earthquake occurred 34 hr later and ruptured perpendicular to the former, along a northwest trending series of faults, arresting at the Garlock Fault. The rupture process for both events is spatially and temporally complex, including several discontinuous fault planes in addition to the delay between events. InSAR line-of-sight observations from Sentinel-1 ascending track 64 (look direction indicated by black arrow) overlain by assumed fault traces for the 4 July M w 6.4 and 6 July M w 7.1 earthquakes, in pink and cyan, respectively. Centroid moment tensor solutions for the two events (globalcmt.org) are shown in the same corresponding colors with lines connecting to their respective epicenter locations. The white arrow points to the location of a 1.3-km step over in the fault structure of the 4 July M w 6.4 event. The inset map (bottom right) shows the regional context of the study area in California (CA), south of Nevada (NV), and north of Baja California, Mexico (BC). The dextral San Andreas Fault (SAF) passes through the southwest corner of the mapped region, with the sinistral Garlock Fault (GF) transecting the study region. The brown shaded region denotes the approximate area of the diffuse Walker Lane (WL) dextral deformation zone.
The close spatial and temporal spacing of these two Ridgecrest earthquakes naturally raises questions about potential interaction between the two events. Earthquakes relieve stress accumulated by tectonic processes, yet cause a heterogeneous redistribution of regional stress that can increase stress on a local scale, pushing certain faults closer to seismic failure. Faults that are brought to failure in part due to the local stress change caused by another earthquake are said to be statically triggered; the initial earthquake either reduces the effective normal stress clamping the fault or increases the shear stress that promotes failure (e.g., King et al., 1994;Freed, 2004). Coulomb static stress change calculations have been used to explain earthquake sequences in eastern and southern California, including another Mojave-area sequence that began in April 1992 with the M w 6.1 Joshua Tree earthquake, followed two months later by the June 1992 M w 7.3 Landers and M w 6.3 Big Bear earthquakes, which were separated by just 3 hr (King & Cocco, 2001). Similarly, the 1999 M w 7.1 Hector Mine earthquake is thought to be triggered by the nearby 1992 M w 7.3 Landers earthquake, in a case of delayed static triggering due to a combination of coseismic and postseismic influences (e.g., Freed & Lin, 2001). In another triggering mechanism, rupture that initiates due to stress change from the passage of seismic waves associated with an earlier earthquake is said to be dynamically triggered. Near-field dynamic stress transfer has been observed to allow rupture to jump across fault segments separated by as much as 5 km (Freed, 2004). Immature fault systems may exhibit dynamic triggering more frequently than mature fault boundaries because they consist of families of disconnected structures, rather than having well-organized throughgoing faults (e.g., Gomberg, 1996).
The rupture kinematics in immature fault systems have been documented previously for earthquakes in the ECSZ and other tectonic environments. The 1999 M w 7.1 Hector Mine, California earthquake ruptured multiple fault planes, away from the main plate boundary, at a relatively low rupture velocity of 2.2 km/s (Ji et al., 2002). Similarly, the 1992 M w 7.3 Landers, California earthquake ruptured five distinct fault segments (Sieh et al., 1993) at an average velocity of 2.7 km/s (Peyrat et al., 2001). Elsewhere, the 2010 M w 7.2 El Mayor-Cucapah earthquake, at the California-Mexico border, ruptured multiple fault planes at an average velocity of 2.5 km/s (Wei et al., 2011). The 2012 M w 8.6 Wharton Basin sequence ruptured six young oceanic faults (Yue et al., 2012) at speeds between 1.5 and 2.5 km/s (e.g., Wei et al., 2013;Yue et al., 2012) and has contributed to the formation of a discrete plate boundary between the Indian and Australian plates (Hill et al., 2015). The 2016 M w 7.8 Kaikoura earthquake is particularly complex, rupturing more than 10 discrete faults, including both strike slip and thrust faulting motion, at low rupture velocities between 1.4 and 2.0 km/s (e.g., Cesca et al., 2017;Hollingsworth et al., 2017;Zhang et al., 2017). This low rupture speed is likely attributable to immature faults being less compliant, more geometrically complex, and rougher than their more developed counterparts (Perrin et al., 2016). In contrast, high rupture speeds, including sustained supershear rupture, are typically associated with well-developed faults devoid of splays or geometric complexities (e.g., Bouchon et al., 2010). For example, the 2018 M w 7.5 Palu, Indonesia earthquake reached supershear rupture in a well-developed fault damage zone with only large-scale fault bends (e.g., Bao et al., 2019;Socquet et al., 2019). The 2002 M w 7.9 Denali, Alaska and 2001 M w 7.8 Kunlun, Tibet earthquakes propagated mostly unilaterally along long, simple, established faults planes at high rupture speeds of 3.2 and 3.4 km/s, respectively (Ozacar & Beck, 2004).
Here we investigate the spatial and temporal distributions of slip during the M w 6.4 and M w 7.1 Ridgecrest earthquakes. We will show that the events exhibit characteristics representative of an immature fault in the process of coalescing to a simple geometry. We argue that both dynamic triggering-within faults that ruptured during the initial M w 6.4 earthquake-and static triggering of the M w 7.1 earthquake by the M w 6.4 rupture, likely occurred. We demonstrate a novel simultaneous kinematic slip inversion methodology that allows us to use Interferometric Synthetic Aperture Radar (InSAR) data that spans both earthquakes along with seismic (strong-motion accelerometer) and geodetic (high-rate Global Navigation Satellite Systems (HR-GNSS)) data sets that recorded each earthquake separately. We compare to recently published slip models that employed a variety of data types (InSAR: Barnhart et al., 2019;GNSS: Ross et al., 2019; GNSS/seismic: Liu et al., 2019). In our new approach, we use all of these different data types to solve for the kinematic slip model of both events concurrently, in a single inversion step.

Data
The kinematic slip inversion considers several distinct observational data sets. First, we select 12 strong motion and 12 GNSS sites based on proximity and azimuthal coverage of the two earthquakes ( Figure 1).

Geophysical Research Letters
We include static offsets estimated by the Nevada Geodetic Laboratory (geodesy.unr.edu) as well as high-rate (5 Hz) GNSS displacements from the Network of the Americas (unavco.org) postprocessed using precise point positioning  to achieve centimeter-level precision in the horizontal components and~3-cm precision in the vertical direction (Melgar et al., 2019). We low-pass filtered the HR-GNSS waveforms with a corner frequency of 0.5 Hz. We integrate strong-motion acceleration waveforms from the Southern California Seismic Network (scedc.caltech.edu) to velocity time series, and band-pass filter between 0.05 and 0.5 Hz. Finally, we use InSAR line-of-sight measurements from two acquisitions of the Sentinel-1 satellites operated by the European Space Agency. For ascending track 64 ( Figure 1) and descending track 71 ( Figure S1), the interferometric pairs span the dates 04-10 and 04-16 July 2019, respectively. The Sentinel-1 data processing techniques (Chen & Zebker, 2002;Lohman & Simons, 2005;Sandwell et al., 2016;Xu et al., 2016Xu et al., , 2017 are summarized in Text S1.

Kinematic Slip Inversion Methodology
We demonstrate a novel simultaneous kinematic slip inversion methodology that permits joint inversion of data types that observe just one of the two earthquakes as well as data types that include contribution from both events in a single observation. More specifically, we use static GNSS displacements, HR-GNSS displacement waveforms, and accelerometer-derived velocity waveforms, each of which observe the two earthquakes separately. We also use InSAR line-of-sight measurements that include deformation from both the 4 July M w 6.4 and the 6 July M w 7.1 events in each interferometric pair. We invert these two categories of data simultaneously to estimate a kinematic slip model of the two events in a single inversion step. We use the gradient of the interferometric phase ( Figure S2; Xu et al., 2020) to identify discontinuities in the deformation associated with surface faulting and define six fault planes-two related to the 4 July event and four related to the 6 July event (Figure 1)-that are assumed to be vertically dipping (corroborated by the moment tensor solutions). An alternate model has been proposed in which one of the northwest trending faults we have assigned to the 6 July event is instead associated with the 4 July event (e.g., Ross et al., 2019). We consider this geometry as well, in the supporting information. Because the surface traces have a complex curvilinear geometry, the fault planes are discretized into an irregular triangular mesh with subfaults with an average area of 2.9 km 2 ( Figure 2). The total slip is subject to Tikhonov minimum-norm regularization, and the regularization parameter chosen is consistent with a classic L-curve method and the expected total seismic moment ( Figure S3).
Our approach uses the multitime window method (Ide et al., 1996) to solve for the distribution of slip in time and space (implementation described in Melgar and Bock (2015)). We use five 50% overlapping triangleshaped source time functions of rise time 0.6 and 1.3 s for the M w 6.4 and M w 7.1 events, respectively. This allows variability in rupture speed across the faults. From the moment tensor solutions, the primary motion of the 4 July M w 6.4 event is left-lateral strike slip, while the 6 July M w 7.1 event motion is primarily rightlateral strike slip. We therefore constrain the rake of the slip along the M w 6.4 faults to be between −45°a nd 45°from north, and the slip along the M w 7.1 faults to be between 135°and 225°from north. We calculate Green's functions using the frequency-wave number technique (Zhu & Rivera, 2002). We assume a onedimensional velocity model from the Southern California Earthquake Center Mojave-area velocity model ( Figure S4) and a maximum rupture velocity of 2.0 km/s, which is roughly 0.55 times the shear wave velocity at the hypocentral depths of the earthquakes. We performed inversions at a range of rupture velocities between 1.4 and 3.0 km/s (see Figures 3a and S5) and found 2.0 km/s to be a mutually favored velocity for both the M w 6.4 and M w 7.1 events when considered independently; we prescribed that constraint to the simultaneous inversion as well. The surface trace shows that the 4 July M w 6.4 event ruptured two discrete fault planes across a 1.3-km step over. Rather than assume that the timing of rupture across the step over is continuous, we create a grid of potential rupture onset locations on the secondary fault to determine which location of secondary rupture initiation best fits the observational data (Figure 3a).
We prescribe weights such that each observation of the same data type has the same weight, yet there is variation between the different observation types. The final weighting was informed by both the data precision  Figure 2) after step over from event hypocenter (on fault segment F2). Total root-mean-square (RMS) difference for inversions of the 4 July event data only (GNSS and strong motion data, no InSAR) with prescribed maximum rupture velocities, Vrup, from 1.8 to 2.4 km/s and 25 potential secondary fault initiation locations, denoted by solid squares. The green square denotes the subfault on segment F1 closest to the event hypocenter (green star on segment F2). The white square denotes the case with minimum total RMS. (b) Snapshots of slip rate of the 4 July M w 6.4 event (fault segments F1 and F2) at six illustrative time frames after origin time (OT) assuming dynamic triggering at the location of minimum total RMS in (a). Event hypocenter is represented by the green star on fault segment F2, and secondary rupture initiation point is represented by yellow star on fault segment F1. and the L2 vector norm of the data type, to ensure that one data type does not overwhelm the inversion (see description in Melgar et al. (2016)). Text S2 describes these choices in detail and Table S1 has the final relative weights used in the inversion.
Simultaneous inversion of the two events requires careful organization of the Green's function matrix, G, within the inverse problem d = Gm, where d is the observational data set and m is the vector of model parameters (the slip on each subfault). In our formulation, d is a column vector of data ordered in the following way: 4 July static offsets, displacement waveforms, and velocity waveforms, 6 July static offsets, displacement waveforms, and velocity waveforms, InSAR line-of-sight observations from descending track 71, and finally, from ascending track 64. The G matrix must therefore properly relate each data set to the model parameters of the relevant earthquake(s). We build separate Green's function matrices for the 4 July-only and 6 July-only data sets (i.e., GNSS static offsets, HR-GNSS displacement waveforms, and accelerometer-derived velocity waveforms), and for the two InSAR data sets, which are affected by the model parameters of both events. We concatenate these matrices to form the full G matrix. An example for a single time window inversion is given by equation 1. For the multiwindow case, as we have implemented here, the G and d matrix rows are simply repeated n times, where n is the number of windows.

Coulomb Stress Change Methodology
We evaluate the Coulomb stress change caused by the 4 July M w 6.4 earthquake on the faults associated with the M w 7.1 earthquake nearly 34 hr later, using our best fitting slip model of the M w 6.4 event. We use the formulation of Lin and Stein (2004) and Toda et al. (2005) for a homogeneous half-space and assume a rigidity of 36.6 GPa, a Poisson's ratio of 0.25, and a coefficient of friction of 0.4 (see Text S3 for details).

Kinematic Slip Inversion
The inversion estimates the seismic moment of the 4 and 6 July events to be 4.37 × 10 18 Nm (M w 6.36) and 4.41 × 10 19 Nm (M w 7.03), respectively. Comparisons between the observed and modeled data are available in Figures S1 (InSAR), S6 (GNSS static offsets), and S7 (HR-GNSS and seismic waveforms). The InSAR residuals ( Figure S1) can be improved with higher weighting (Figure S8), at the expense of larger and rougher surface slip. Moreover, the InSAR observations likely include postseismic deformation on the order of these residuals ( Figure S9), and thus, we have selected a final model that allows for smoother slip patterns. The respective rupture durations are 12.2 and 23.8 s (Figure 2). The maximum moment rate of the first earthquake occurs about 60% into the rupture duration, while the larger event reaches maximum moment rate early on, at only~25% of the total rupture duration. An alternate model in which fault segment F4 (see Figure 2) is assigned to the 4 July event rather than the 6 July event shows similar features and we do not distinguish which fault model is preferred ( Figure S10).
The 4 July M w 6.4 event began on a small fault segment, jumping a step over to continue onto a larger fault plane to the southwest (Figure 1). We evaluated this earthquake independently of the M w 7.1 mainshock (excluding the InSAR data, which spanned both events) in order to identify how the rupture transitioned to the secondary fault in time and space. We considered the case of dynamic triggering of the secondary fault from the passing seismic waves generated from slip on the initial fault. We constructed a grid of potential rupture initiation locations along the secondary fault (Figure 3a), including an additional node at the subfault on the secondary plane closest in distance to the event hypocenter, implicitly considering the case for continuous rupture, without any triggering. We constrain the rupture onset times such that rupture begins on the secondary fault at the time when the S wave from the event hypocenter passes the subfault of interest. We calculate the overall root-mean-square difference between the observations and model for each of these potential initiation locations. While our temporal resolution does not allow us to rule out continuous rupture entirely, we find that the data prefer that rupture on the secondary fault initiates at roughly 15-km depth, 6 km from the fault edge nearest the event hypocenter, and with a maximum rupture velocity of 2.0 km/s (Figure 3a). This dynamic triggering of slip on the secondary fault plane results in dual slip pulses migrating across the two fault planes (Figure 3b). Due to our limited spatial and temporal resolution, we have tested only the case of dynamic triggering from the passing S wave, but note that triggering from an earlier or later seismic phase would likely result in a shift in the preferred location of rupture initiation across the step over. Our inversion uses longer-period data; thus, we cannot discern this clearly. The fault segments associated with the 6 July M w 7.1 event are all connected, without any step overs in the fault geometry; thus, we assume continuous rupture from the hypocenter for this event.
The main asperity of the 4 July M w 6.4 event is centered around 7-km depth toward the center of the secondary fault plane, F1 (Figure 2). Slip reaches~1 m in this location; however, the maximum slip amplitude (1.1 m) occurs in the shallow part of the initial fault plane, F2. Our model suggests surface offsets up to 0.8 m resulting from this first event. The maximum slip associated with the M w 7.1 event reaches~4 km in the main asperity, which is located in the area between the M w 7.1 hypocenter and the cross-cutting M w 6.4 fault planes (see Figure 1). Rupture appears to be concentrated in two lobes on the throughgoing M w 7.1 fault (F3), with slip reaching >20-km depth. The shorter fault segments associated with the M w 7.1 mimic the patterns on the adjacent segment of the throughgoing fault plane. In order to determine whether the features exhibited in the inversion are robust, we performed a jackknife test, which shows that the locations of main asperities (Figures 2 and S11) are well constrained, but our model is less certain of the shallowest slip, slip at the fault edges, and below the main asperities. Details are available in Text S4 and Figure S11.
Although there are slight differences in fault geometry and parameterization, our slip model shares many features of recently published static and kinematic inversions, including the two lobes of slip and a maximum slip of~4 m for the 6 July event and just over 1 m of maximum slip for the 4 July event (Barnhart et al., 2019;Liu et al., 2019).

Coulomb Stress Change
The slip during the M w 6.4 event results in a reorganization of stress in the surrounding region. Our analysis shows a decrease in Coulomb stress on the main throughgoing fault of the M w 7.1 event in the area that is cross-cut by the M w 6.4 faults (Figure 4, fault F3). Similarly, the western splay at the south end of the M w 7.1 rupture (Figure 4, fault F6), located just south of the M w 6.4 fault, experienced a decrease of Coulomb stress in response to the M w 6.4 rupture. There is, however, a localized area of increased Coulomb stress (<0.1 MPa) at the location of rupture initiation of the M w 7.1 event. Aftershock studies demonstrate that small stress increases on this scale (<0.1-0.3 MPa) are sufficient to bring critically stressed faults to failure (e.g., Freed, 2004). The calculated Coulomb stress change is directly affected by the result of the M w 6.4 slip model, which is in turn affected by our choices of data weighting. To ensure that the increased Coulomb stress at the location of the M w 7.1 event hypocenter is a robust feature, we repeat this analysis using slip models in which we vary the weight assigned to the InSAR data. We find that the increased Coulomb stress at this location is a persistent feature ( Figure S12). In the alternate case where subfault F4 is assumed to have ruptured during the 4 July event, the resulting Coulomb stress along segment F3 shows that the 6 July hypocenter is at the border between areas of increased and decreased Coulomb stress ( Figure S10f).

Discussion
Complex, multifault, and slow-velocity rupture appear to be common features of earthquakes in immature fault systems transitioning into major tectonic boundaries. The Ridgecrest sequence seems to display this behavior as well. Our kinematic slip inversion prefers a slow rupture velocity of 2.0 km/s (~55% of shear wave speed at the hypocenter depth), consistent with previous studies noted in the introduction that suggested that young faults tend to rupture slower than mature faults. We observe evidence of dynamic triggering across a 1.3-km step over during the 4 July M w 6.4 event. Interestingly, the model prefers a case where the rupture across the step over begins near the center of the second fault, rather than at the edge nearest the first activated fault.
From calculation of the static Coulomb stress change imposed by the 4 July M w 6.4 event, we find a heterogeneous static stress change on the fault planes associated with the subsequent 6 July M w 7.1 event, wherein the hypocenter of the M w 7.1 event is located directly in a localized region of increased Coulomb stress. The stress perturbation is small (<0.1 MPa), but consistent with previous studies that demonstrate small stress changes can be sufficient to initiate rupture on critically stressed faults. We note that while we have not addressed the M w 6.4 aftershock sequence here, it is likely that those events also have an effect on the stress change in the area of the M w 7.1 nucleation. Similar analysis that includes the effect of the M w 6.4 aftershock sequence agrees with our interpretation that the M w 7.1 hypocenter is in a region of increased stress due to the previous events (Barnhart et al., 2019).
While mature faults have generally simpler geometries, incipient shear zones are often more geometrically complex, segmented, and involve multiple fault planes with variable orientations that may not be optimally oriented for failure (Crider & Peacock, 2004). Independently, these fault strands are small and incapable of hosting large earthquakes characteristic of major plate boundary faults. Through time, these immature fault strands tend to localize deformation and ultimately link, forming throughgoing fault structures capable of hosting large-magnitude earthquakes (Manighetti et al., 2007;Perrin et al., 2016;Thomas et al., 2013;Wesnousky, 1988). The large magnitude, tectonic setting, slow rupture velocity, and structural complexity of the Ridgecrest sequence suggest that these events are a manifestation of the increasing contribution of the Walker Lane and ECSZ to the overall plate boundary deformation. The earthquake interactions presented here have important implications for earthquake hazard models, as small immature faults may not be capable of hosting large-magnitude earthquakes, yet cascading rupture through multiple fault strands can accommodate moderate-to large-magnitude earthquakes in areas where such events are unexpected.

Conclusions
We apply an innovative kinematic slip inversion approach that allows simultaneous inversion of InSAR data that spans the two largest Ridgecrest sequence earthquakes alongside data sets that recorded each earthquake separately (strong motion accelerometer, HR-GNSS, and GNSS static offset) in a single inversion step. Our kinematic model is indicative of dynamic triggering across a step over associated with the 4 July M w 6.4 event, and we demonstrate that the M w 6.4 event may have statically triggered the M w 7.1 event by increasing the Coulomb stress near the hypocenter of the M w 7.1 earthquake. Our analysis of the kinematic rupture processes of the Ridgecrest sequence supports the concept that the Eastern California Shear Zone and Walker Lane comprise an immature, incipient fault zone, to which the San Andreas Fault is transferring an increasing amount of the overall plate motion between the Pacific and North American Plates.