Aseismic Fault Slip During a Shallow Normal-Faulting Seismic Swarm Constrained Using a Physically Informed Geodetic Inversion Method

Improved imaging of the spatio-temporal growth of fault slip is crucial for understanding the driving mechanisms of earthquakes and faulting. This is especially critical to properly evaluate the evolution of seismic swarms and earthquake precursory phenomena. Fault slip inversion is an ill-posed problem and hence regularization is required to obtain stable and interpretable solutions. An analysis of compiled finite fault slip models shows that slip distributions can be approximated with a generic elliptical shape, particularly well for M ≤ 7.5 events. Therefore, we introduce a new physically informed regularization to constrain the spatial pattern of slip distribution. Our approach adapts a crack model derived from mechanical laboratory experiments and allows for complex slipping patterns by stacking multiple cracks. The new inversion method successfully recovered different simulated time-dependent patterns of slip propagation, that is, crack-like and pulse-like ruptures, directly using wrapped satellite radar interferometry (InSAR) phase observations. We find that the new method reduces model parameter space, and favors simpler interpretable spatio-temporal fault slip distributions. We apply the proposed method to the 2011 March–September normal-faulting seismic swarm at Hawthorne (Nevada, USA), by computing ENVISAT and RADARSAT-2 interferograms to estimate the spatio-temporal evolution of fault slip distribution. The results show that (a) aseismic slip might play a significant role during the initial stage and (b) this shallow seismic swarm had slip rates value decomposition (TSVD) to reject model space basis vectors associated with small singular values. Instead of regularizing the volume variation itself, they minimized the volume change rate, to avoid large discontinuities. Here, we improve previous methods by (a) regularizing the fault slip distribution using a prescribed parameterization derived from a laboratory-based crack model and (b) introducing a statistically optimal truncation criterion that allows to automatically separate signal and noise in the spatio-temporal fault slip distributions. We demonstrated the validity of this approach using synthetic experiments and comparing it against a compilation of published slip distribution models. Finally, we applied the new proposed methodology to the 2011 Hawthorne seismic swarm (Nevada, USA). The 2011 Hawthorne seismic swarm is located at the central Walker Lane, which accommodates the Pacific-North American transform plate motion by oblique-normal faults and block rotations (Thatcher et al., 1999; Wesnousky, 2005). The 2011 Hawthorne swarm consists of 10 M4+ events, and the largest earthquake among them is an M4.6 event (Smith et al., 2011; Zha et al., 2019); a recent study using satellite images reveals clear surface deformation signals before the M4.6 event, and the geodetic moment is much higher than the seismic moment, indicating that aseismic slip dominates the fault behavior (Jiang & González, 2021). By applying our proposed methodology, we retrieved the fault-slip spatio-temporal evolution, and explored the interactions between the fault slip and the seismicity. results slip larger coseismic ruptures of fault fault slip lower

slow-slip triggers and/or precedes the occurrence of a seismically dynamic rupture. Within the scope of increasing our capacity to distinguish between the earthquake nucleation models, a promising venue is to increase our ability to image how fault slip evolves in space and time. Although fault slip evolution is not necessarily the only cause of seismicity migrating, improvements in this direction may provide crucial data to examine hypotheses for earthquake nucleation mechanisms.
Fault slip imaging improvements are particularly desirable to estimate (seismic and aseismic) slip propagation parameters, such as slip rate, and gain deeper insights into the physics controlling regular earthquakes and slow-slip phenomena. Regular earthquakes are known to show peak and average slip rates of the order of 1 and 0.1 m/s (Takenaka & Fujii, 2008). While slow-slip phenomena show much lower slip rates, for example, slow slip events (SSEs), fault creep, or slip related to fluid injection. For example, in the case of SSEs in subduction zones, the peak slip rates vary around 0.1 ∼ 3 cm/day (Bletery & Nocquet, 2020;Ozawa et al., 2019;Radiguet et al., 2011;Rousset et al., 2019). In the case of the episodic creep event, the slip rates in continental faults are 0.5 ∼ 3 cm/year (Hussain et al., 2016;Jolivet et al., 2012;Schmidt et al., 2005;Scott et al., 2020). In fluid injection experiments, the slip rates have been observed to be much higher, up to 4 × 10 −3 mm/s (35 cm/day) (Guglielmi et al., 2015).
Hence, to evaluate (seismic and aseismic) fault slip characteristics, a better description of how fault slip propagates in space and time is necessary. Including complex propagation patterns of fault slip such as pulse-like and crack-like ruptures (Lambert et al., 2021;Marone & Richardson, 2006). Such patterns have been observed during regular earthquakes but are also associated with slow-slip phenomena: with slow slip transients migrating further away from where they started along strike (or dip) or remain stationary through time. Observations of some SSEs and "Episodic Tremor and Slip" (ETS) show pulse-like rupture characteristics with elongated slipping areas, for example, the Cascadia subduction zone (Michel et al., 2019), and with along strike migration speeds of ∼10 km/ day (Rousset et al., 2019;Wech et al., 2009). In contrast, slip propagation of meter-scale fluid injection experiments indicates stationary patterns: Bhattacharya and Viesca (2019) proposed a model in which the slip grows as an expanding ellipse, with the injection point as the slipping center. The latter phenomenon is also found in some SSEs on subduction zones, for example, the deep Manawatu and Kaimanawa SSEs on the Hikurangi subduction zone (Wallace, 2020). Here, we aim to improve fault slip mapping in space and time to contribute to the advancement of the study of fault slip processes using, yet underutilized, satellite InSAR observations.
In this research, we developed a new method to interpret directly wrapped phase InSAR observations to estimate the spatio-temporal fault slip, in particular, in the context of a favorable tectonic setting, continental seismic swarms (e.g., small-amplitude surface deformation signals and/or phase discontinuities due to surface ruptures). InSAR has been used to map surface displacements with high spatial resolution and subsequently model fault slip. But so far, it is more common to estimate static slip distributions than jointly invert for the time series of slip evolution (Floyd et al., 2016;Ingleby et al., 2020). The problem of retrieving time series of source parameters from nonsimultaneous and temporally overlapped multi-sensor observations is ill-posed; however, the oscillations of the solution caused by the rank deficiency of this problem can be reduced by applying regularization or temporal filtering (Samsonov & D'Oreye, 2012). Grandin et al. (2010) introduced a temporal smoothing scheme as an additional constraint to retrieve the time series of magma volume changes. Additionally, González et al. (2013) used truncated singular value decomposition (TSVD) to reject model space basis vectors associated with small singular values. Instead of regularizing the volume variation itself, they minimized the volume change rate, to avoid large discontinuities. Here, we improve previous methods by (a) regularizing the fault slip distribution using a prescribed parameterization derived from a laboratory-based crack model and (b) introducing a statistically optimal truncation criterion that allows to automatically separate signal and noise in the spatio-temporal fault slip distributions. We demonstrated the validity of this approach using synthetic experiments and comparing it against a compilation of published slip distribution models. Finally, we applied the new proposed methodology to the 2011 Hawthorne seismic swarm (Nevada, USA). The 2011 Hawthorne seismic swarm is located at the central Walker Lane, which accommodates the Pacific-North American transform plate motion by oblique-normal faults and block rotations (Thatcher et al., 1999;Wesnousky, 2005). The 2011 Hawthorne swarm consists of 10 M4+ events, and the largest earthquake among them is an M4.6 event (Smith et al., 2011;Zha et al., 2019); a recent study using satellite images reveals clear surface deformation signals before the M4.6 event, and the geodetic moment is much higher than the seismic moment, indicating that aseismic slip dominates the fault behavior (Jiang & González, 2021). By applying our proposed methodology, we retrieved the fault-slip spatio-temporal evolution, and explored the interactions between the fault slip and the seismicity.

Static Fault Slip Models
Slip inversions with kinematic models are ill-posed problems in which the solution is nonunique and unstable, and unphysical slip distributions can be estimated by least squares algorithms, that is, extremely rough oscillatory slip distributions. Harris and Segall (1987) introduced Laplacian smoothing as the regularization scheme. This minimizes the second derivative of slip and can prevent cases with large stress drops. Du et al. (1992) plotted a trade-off curve for misfit as a function of slip roughness, and manually picked a smoothing factor within the inflection point of the curve to find an optimal balance between data fit and model roughness. Matthews and Segall (1993) determined the optimal smoothing factor in the trade-off curve objectively by implementing the cross-validation method. Much later, Fukahata and Wright (2008) and Fukuda and Johnson (2008) introduced the Bayesian approach, Akaike's Bayesian Information Criterion (ABIC), to solve the slip distribution. While Fukahata and Wright (2008) emphasized the significance of fault geometry as a nonlinear constraint, Fukuda and Johnson (2008) overcame the deficiencies of ABIC with positivity constraints, and then applied the adapted ABIC to simultaneously estimate the slip distribution and smoothing parameter objectively in a Bayesian framework. Fukuda and Johnson (2010) then devised a mixed linear-non-linear Bayesian inverse formulation and extended their work for the joint slip and geometry inversion. In response, Minson et al. (2013) argued that the nonphysical regularization scheme (i.e., Laplacian smoothing) is unnecessary, and developed a fully Bayesian approach to sample all possible families of models compatible with the observations, via a parallel computing framework. Ragon et al. (2018) further extended the work of Minson et al. (2013) and accounted for the uncertainty in fault geometry. Instead of Laplacian regularization, Amey et al. (2018) developed an inversion package slipBERI, and incorporated self-similarity, characterizing the seismic slip distribution in real earthquakes, as a prior assumption within the Bayesian inversion of earthquake slip.
All the previous methods are based on kinematic models that do not take into account the relationship between stress and slip in the fault. Alternatively, dynamic source models satisfy physical constraints on the propagation of shear fractures on Earth, but few dynamic source models are considered to constrain the slip inversions. As an alternative, Di Carli et al. (2010) proposed using elliptical patches to describe the slip distribution in the kinematic and dynamic inversion of near-field strong motion data at low frequencies. Soon afterward, Sun et al. (2011) put forward a mechanical slip inversion, imposing a uniform stress drop on the fault plane. The resulting slip distribution is inherently smooth, so the smoothing norm and the smoothing factor are unnecessary. Tridon et al. (2016) assumed a circular stress patch in volcano research, inverting the displacement for shear and normal stresses simultaneously, along with the fault geometry.
In this study, we apply a new methodology named Geodetic fault-slip Inversion using a physics-based Crack Model (GICMo) (Jiang et al., 2022). In this method, we take advantage of a one-dimensional analytical crack model proposed by Ke et al. (2020). The model was theoretically and experimentally validated in self-contained ruptures within a 3-m-long saw-cut granite fault. This new crack model features non-singular (finite) peak stresses at the rupture tip. In Jiang et al. (2022), we expanded the one-dimensional model into two dimensions to produce elliptical fault slip shapes/patches. We assume that one of the focal points of the ellipse is the crack center (with the maximum slip) and the elliptical perimeter to be the crack tip. Therefore, the slip distribution on the fault plane is controlled by a very compact and reduced set of parameters. The geodetic-inverted fault slip infers that it is possible that the crack center can be located at the rupture center, for example, the 2009 L'Aquila earthquake (Walters et al., 2009). To adapt to this possibility, we relax the constraint that the maximum slip should coincide with the crack center location, and allow it to move along the x axis inside the ellipse. Hence, our crack model contains only eight parameters as demonstrated by Equation 1 and where s is the slip distribution; x 0 , y 0 are the locations of the crack center; a and e are the semi-major axis and eccentricity of the ellipse; α is the ratio controlling the location of the crack center along x axis: the crack center is located at the ellipse center, left/right vertices when α = 0, −1/1; λ is the ratio controlling the displacement transition from the center to the edge of the elliptical crack; d max is the maximum slip; θ is the rake angle.
In the GICMo method, once the crack model parameters are provided, the slips for all fault patches are then determined based on the two-dimensional crack model discussed above. Then, the fault slip distribution is 10.1029/2021JB022621 4 of 21 forward modeled to estimate surface displacement. Following Jiang and González (2020), a misfit function is constructed based on the wrapped phase residuals and the weighting matrix. The misfit function is then regarded as the likelihood function fed into the Bayesian process to retrieve the posterior distribution of crack model parameters. In the Bayesian process, the Markov chain Monte Carlo algorithm is adopted as the probability sampling approach based on the Metropolis-Hasting rule.
Here, we design a synthetic static slip to compare the performance of our method, GICMo, and a state-of-the-art method, slipBERI (Amey et al., 2018). The geodetic inversion package, slipBERI, solves for fault slip with GNSS and unwrapped InSAR phases in a Bayesian approach using von Karman regularization, and simultaneously solves for a hyperparameter that controls the degree of regularization. A normal fault with pure down-dip slip is simulated as the synthetic fault model. To imitate the slipping patterns observed in the published finite-source rupture models SRCMOD (Mai & Thingbaijam, 2014) (e.g., Bennett et al., 1995;Elliott et al., 2010;Ichinose et al., 2003), the inner region is a square area with a larger displacement, and the outer region is an annulus area with a smaller displacement ( Figure 2). Due to the difference in the ingestion data, the synthetic phases are unwrapped phases for slipBERI and wrapped phases for GICMo. The displacement phase is forward calculated based on the synthetic fault slip distribution and the dislocation model. To increase its resemblance to reality, decorrelation and atmosphere noises are simulated and added, whose amplitudes are 10% of 2π for wrapped phase cases or the peak amplitude of the deformation phase for unwrapped phase cases, which is based on the signal-to-noise ratio from a real interferogram in Section 4 (RS2-20110322-20110415). The simulated noise-plus-deformation interferogram is resampled with a quadtree algorithm within the downsampled unwrapped and wrapped phases (Bagnardi & Hooper, 2018;Jiang & González, 2020). In addition, the covariance matrix is estimated based on the phase in the far-field. Finally, the downsampled phases and covariance matrix are fed into slipBERI and GICMo to retrieve the slip distributions. Figures 2b-2d show the modeled slip distribution inverted by GICMo and slipBERI, and Figure  S1 in Supporting Information S1 shows the modeled phase and phase residuals. The conclusions are listed below.
1. Both GICMo and slipBERI provide the first-order accuracy of the slip distribution, including the locations of the crack center and the magnitude of the slip peak. 2. We interpolate the slip distribution onto a 0.5 × 0.5 km patch mesh, and calculate the root-mean-square error (RMSE) of the slip distribution compared with the synthetic slip distribution. We find that the RMSEs are 1.5 cm for the one-ellipse model, 2 cm for the von Karman smoothing model, and 3 cm for the Laplacian smoothing model, which are approximately similar. However, the great advantage is that the parameters to be solved in GICMo are independent of the fault mesh discretization, and the number of parameters is 30 times less in this case than 201 in slipBERI for this case.

Bayesian Inversion of Fault Slip Time Series Using a Physics-Based Crack Model (Time-GICMo)
The temporal evolution of fault slip is critical to understanding the driving mechanism of slow slip. It is difficult to find one SSE where one interferogram can coincidentally capture the beginning and the ending of the activity. Instead, a common scenario is that the slip increment is captured by interferograms. In this section, we develop We consider a system of N increments of fault slip (Δs n ∈ [Δs 1 , …, Δs N ] between dates and ) based on the nonlinear inversion estimation from the corresponding wrapped interferogram, and the raw images of interferograms are acquired at M unique dates (t ∈ [t 1 , …, t M ]). The aim is to solve for the temporal evolution of fault slips (s ∈ [s 1 , …, s M ]) for each date. We assume that the slip rate between adjacent dates (v m ∈ [v 1 , …, v M−1 ]) is constant, so the slip increment Δs n can be expressed by the sum of fault slip increment between adjacent dates, . The linear expression for N increments of fault slip is shown in Equation 2, as illustrated by González et al. (2013) (2)

of 21
where P is the observation vector, Q is the unknown vector, and B is the designed matrix. Considering there are N increments of fault slip, the matrix dimension is (N × 1) for P, (N × (M − 1)) for B, and ((M − 1) × 1) for Q. Then, we decompose matrix B by using the SVD methods, where U is an orthogonal matrix with columns that are the basis vectors of the data space (N × N), V is an orthogonal matrix with columns that are the basis vectors spanning the singular values of the model ((M − 1) × (M − 1)), and S is a diagonal matrix of the singular values ((N × (M − 1)) × 1). A solution for this problem can be obtained as follows, If rank (B) < m, the solution obtained using the SVD technique may contain numerical instabilities when there are small singular values. In this case, a more stable solution can be achieved using the TSVD method (Aster et al., 2019), which rejects model space basis vectors associated with small singular values, up to a certain threshold. As an improvement upon González et al. (2013), we apply an optimal hard threshold for singular values truncation proposed by Gavish and Donoho (2014). Gavish and Donoho (2014) proposed that the optimal hard threshold for singular value is 4∕ √ 3 of the median singular value. This criterion is empirically proven to be the best hard thresholding, independent of model size, noise level, or true rank of the low-rank model. This improvement allows us to define the degree of regularization based on an objective criterion, which generates a parsimonious low-rank model solution in the presence of noisy data. Note that in order to retrieve a realistic solution, a non-negative constraint is added in solving for slip rate vector Q implemented by using MATLAB function lsqnonneg (https://uk.mathworks.com/help/optim/ug/lsqnonneg.html). It is physically appropriate because slip along faults rarely re-rupture backwards (Hicks et al., 2020).

Time-Dependent Fault Slip Inversion Experiments
In this section, we describe two experiments to investigate if this method can retrieve pulse-and crack-like rupture propagation patterns in space and time. We tested the performance of the inversion method to recover fault slip evolution from each of the two-ellipse models.
The first synthetic case aims to explore the inversion with overlapping ruptures (Figure 3). Several recent studies have suggested spatial overlap between coseismic slip and afterslip (Barnhart et al., 2016;Bedford et al., 2013;Bürgmann et al., 2002;Johnson et al., 2012;Pritchard & Simons, 2006;Salman et al., 2017;Tsang et al., 2016). A series of overlapping elliptical cracks are simulated in Figure 3a and forward inversion is performed to calculate the surface displacement due to the slip increment between adjacent cracks. We aimed to compare the results based on various geodetic inversion algorithms: (a) the one-ellipse model, as described in Section 2.1, (b) a von Karman regularization algorithm (Amey et al., 2018), and (c) the two-ellipse model with different crack centers. Inversion results are shown in Figures 3b-3d, and the modeled phase and residuals are shown in Figures S2 and  S3 in Supporting Information S1. The main conclusions are as follows: 1. The RMSEs of the fault slip residual is the lowest in results based on the two-ellipse model with different centers. The triangle patch size in the crack model is ∼0.84 km, and the rectangle patch size in slipBERI is 1.5 km. In this way, we interpolated the modeled slip distributions to grid points with 1.17 km spacing, and then calculated the RMSE of the fault slip residual. In each case, the RMSE of slip residuals based on the two-ellipse model with different centers (Figure 3d) are the smallest, and the average RMSE for the one-ellipse model, the von Karman smoothing model and the two-ellipse model are 0.9, 1.6, and 0.6 cm. 2. The two-ellipse model is superior to the one-ellipse model in the F-test for the residual of the interferometric phase. The two-ellipse model has more free parameters, leading to an inherent improvement in the data fit.
To objectively compare the model performances, we use the F-ratio statistic to test the significance of the decrease of residuals between models (Stein & Gordon, 1984). The statistical test checks if the empirical F-ratio (F emp ) is larger than the theoretical (F theory ). In this case, the comparison of the one-ellipse model and two-ellipse model leads to F emp = 72.8 ≫ F theory = 2.6.
10.1029/2021JB022621 7 of 21 The second synthetic case aims to explore the inversion with the containing ruptures ( Figure 4). A growing rupture has been widely observed and studied in fluid injection experiments (Bhattacharya & Viesca, 2019;Cappa et al., 2019;Guglielmi et al., 2015). The rupture center is located at the injection point, and the radius of the slipping zone grows at a rate up to 10 −6 m/s. A set containing elliptical ruptures is simulated in Figure 4a and a forward inversion facilitates the surface displacement calculation. We aimed to retrieve the slip increments from the observed interferometric phase with various methods described above (one-ellipse model, von Karman smoothing model, and two-ellipse model). On noticing that the slip distribution is not well resolved by the two-ellipse model with different centers, we added another constraint to the two-ellipse model so that both cracks share the same center. The inversion results are shown in Figures 4b-4e, and the modeled phase and residuals are shown in Figures S4 and S5 in Supporting Information S1. The main conclusions are as follows: 1. The average RMSE of slip residuals based on various inversion models (one-ellipse model, von Karman smoothing model, two-ellipse model with different centers, and one center) are 1.3, 1.3, 1.0, and 0.8 cm. The one-ellipse model failed because the slip increment in containing ruptures no longer could be described by one complete crack. Indeed, slipBERI showed better performance because it inferred the region with the slip peak. The two-ellipse model with different centers is even better but was not well resolved, for example, the slip increment from t 1 to t 2 (second image in Figure 4c). Therefore, the two-ellipse model with the same center is the most appropriate for reconstructing the cracks' locations, sizes, and maximum slips. 2. In the F-test of the interferometric phase residuals, the two-ellipse model with the same center is superior to the two-ellipse model with different centers, and the one-ellipse model is the least useful.

Regional Tectonics and Seismicity
We apply our algorithm to the 2011 Hawthorne seismic swarm, which occurred on the central Walker Lane ( Figure 5). The Walker Lane is a 500 km-long and 100 km-wide deformation region consisting of N-NW right-lateral shear and extension (Wesnousky, 2005). It is located between the northwest translating Sierra Nevada microplate and the westward extending Basin and Range Province. The Walker Lane accommodates 20% ∼ 25% of the current relative motion (50 mm/year) between the Pacific and North American plates (Argus & Gordon, 1991;Faulds & Henry, 2008). The central Walker Lane accommodates the deformation budget of ∼8 mm/year between the Basin and Range province and the central Sierra Nevada (Bormann et al., 2016). The distributed dextral shear in central Walker Lane is accommodated by oblique-normal faults, block rotations, and partitioning of oblique deformation between sub-parallel normal and strike-slip faults. The total long-term strain rate is 51 nanostrain/year extension directed N77°W and 38 nanostrain/year contraction directed N13°E (Kreemer et al., 2014), much higher than the central Basin and Range (Kreemer et al., 2009).
Being a geologically young and developing fault system, the Walker Lane shows high levels of seismicity over the instrument period, including >10 M6+ earthquakes in the last century. Since 2000, the Walker Lane was struck by a few seismic sequences with some accompanied by aseismic slip evidence. For example, for the 2008 Mogul earthquake sequence, geodetic observation and modeling indicated significant aseismic slip (Bell et al., 2012), and the migration speed of the largest foreshock cluster is consistent with aseismic slip (Ruhl et al., 2016); for the 2014 Virginia City Swarm, migration rate of small earthquakes was consistent with rates observed elsewhere associated with pore fluid diffusion and aseismic creep (Hatch et al., 2020).  The 2011 Hawthorne seismic swarm lasted from March to September and consisted of 10 M4+ earthquakes according to the U.S. Geological Survey (USGS) hypocenter catalog (https://earthquake.usgs.gov/earthquakes/ search/). This sequence occurred in the footwall block of the Wassuk Range segment at the central Walker Lane (Faulds & Henry, 2008), and this segment experiences a significant extension of 1.5 ± 0.3 mm/year (Hammond & Thatcher, 2007). Early moment tensor solutions show the shallow depths in this sequence (Smith et al., 2011), and further hypocenter relocation together with the focal mechanisms of the M4+ events consistently reveal a W-NW-dipping normal fault zone with centroid depths between 2 and 4 km (Zha et al., 2019). The 2011 Hawthorne sequence is close to the Aurora-Bodie volcano (Lange & Carmichael, 1996), but no volcanic signature was observed in near-source seismograms, which infers this sequence is not likely related to the magmatic activity (Smith et al., 2011;Zha et al., 2019). In this research, we were able to identify three stages with respect to the timing of the most energetic event (M4.6) occurred: an initial stage (pre-M4.6 stage) from 15 March to 17 April, a shorter period around the most energetic stage (co-M4.6 stage), and the post-energetic stage (post-M4.6 stage) lasting until 17 September.

Multi-Satellite Geodetic Data Sets
We processed ENVISAT and RADARSAT-2 data and generated eight SAR interferograms to quantify surface displacements ( Figure 6). SAR images were acquired between February and September 2011 from two tracks: one ascending track from the Canadian Space Agency RADARSAT-2 satellite, look angle of 35° and heading angle of 350°; and another descending track from the European Space Agency (ESA) ENVISAT satellite, track 343, look angle of 35° and heading angle of −166°. Interferograms were processed in two-pass differential mode, using a 30 m resolution digital elevation model (DEM) derived from the Shuttle Radar Topography Mission. ENVISAT-ASAR data were processed using Doris software (Kampes et al., 2003) and ISCE software, RADARSAT-2 data using GAMMA software (Werner, 2000). Overall, we obtained 8 short baseline differential interferograms. The computed interferograms have temporal separations ranging from 24 to 120 days. Considering the dominant extensional mechanism and N-S fault striking in this region, the preferred movement direction of the ground displacement is E-W. Consequently, the satellite flight direction favors surface displacement observations in this normal faulting system.
Interestingly, two ascending RADARSAT-2 interferograms during the pre-M4.6 stage indicated clear surface displacement signals (Figures 6d and 6a), ∼4 cm away from satellite line-of-sight motion. In interferograms covering the co-M4.6 stage, it is notable that surface displacement signals were larger in magnitude and located further north with respect to the pre-M4.6 stage (Figures 6b-6f). During the early post-M4.6 stage, surface displacements were detected along a very narrow spatial band with clear phase discontinuities, suggesting surface ruptures (Figure 6g). For one interferogram covering the late post-M4.6 stage (Figure 6h), the phase was dominated by atmospheric noise and no clear deformation signal was detected. Analysis of interferograms suggests that fault slip may have occurred along a fault system with a two-plane geometry, which is consistent with the finding from early moment tensor solutions (Smith et al., 2011).
Note that the two ascending RADARSAT-2 interferograms provide a unique opportunity to look into the preseismic slip, which is not available in other reported cases due to the data limitation. For example, for the 2008 Mogul earthquake swarm, Bell et al. (2012) measured the surface deformation covering the whole earthquake swarm using InSAR and they found that the modeled cumulative geodetic moment is ∼2 times the cumulative seismic moment, indicating a significant portion of aseismic slip. However, they cannot separate the preseismic deforma-  tion signal because there is no available interferogram covering the preseismic stage only. In addition, the GPS observations covering the 2008 Mogul earthquake swarm cannot constrain the preseismic slip well due to the low signal-to-noise ratio in GPS solutions (Ruhl et al., 2016).

Spatio-Temporal Slip Evolution
To develop the kinematic fault model, we first constructed the fault geometry derived from a nonlinear fault inversion of InSAR wrapped phase observations, solving for uniform distribution on rectangular faults (Jiang & González, 2020). A geodetic inversion directly using the interferometric wrapped phase avoids any potential phase unwrapping error ( Figure S6 in Supporting Information S1). The data variance-covariances describing the noise level are calculated based on the covariograms ( Figure S7 in Supporting Information S1) and are used to weight the wrapped phase residuals in the likelihood function as illustrated by Jiang and González (2020). Modeling of a selection of interferograms covering the successive phases confirmed that ground motion could be caused by fault geometry with two distinct planes. During the pre-M4.6 stage, the observed ground motion in the RADARSAT-2 interferogram (22 March Figure 7a). Only one single fault is applied in the modeling above, and the phase caused by the northern subfault is modeled well due to its dominance during the co-and post-M4.6 stages. The residual is relatively larger in the south because of the ignorance of the southern subfault, as shown by the residual phases in Figure S8 in Supporting Information S1. Based on modeled fault geometry in Figure 7a, together with ground motion discontinuities digitized from the interferograms, we constructed a smooth fault plane with uniformly discretized triangular meshes in Figure 7d. These were generated by FaultResampler (Barnhart & Lohman, 2010) and mesh2d (Engwirda, 2014), with a near-uniform side length of around 125 m. Then, a fault slip distribution model with associated uncertainties was estimated. We applied the fault slip inversion method based on a prescribed regularization derived from an experimentally validated physics-based crack model (Jiang et al., 2022). To further investigate the temporal evolution of fault slips with a higher temporal resolution, we invert the fault slip time series using all available interferograms with clear deformation signals. Figure 8 presents the temporal evolution of cumulative slip and slip rate during the 2011 Hawthorne seismic swarm, and Figure S9 in Supporting Information S1 shows the modeled phase and phase residuals. The findings from the inversion results are listed as follows: 1. There were three areas with different spatio-temporal slipping behaviors: a narrow (5 km 2 ) slip area on the southern fault with a high rate (with a lower bound: 1.5 cm/day, or 1.7 × 10 −7 m/s) occurring during the pre-M4.6 stage, a wider (15 km 2 ) slip area with lower average slip (10 cm) on the northern fault that ruptured during the co-M4.6 stage, and a shallow slip area (depth = 1 km) just above the second area during the post-M4.6 stage with a slower average slip rate (with a lower bound: 0.2 cm/day, or 2.3 × 10 −8 m/s). 2. Our results show the aseismic slip mainly occurred on the southern subfault during the pre-M4.6 stage, while the most significant seismic slip hit the northern subfault during the co-and post-M4.6 stages. The results are more consistent with a cascade model of discrete slip patches, rather than a slow-slip model considered as a growing elliptical crack. 3. During the early pre-M4.6 stage (26 February-22 March), the cumulative geodetic moment is 1.7 × 10 16 Nm (equivalent to an M w 4.7 event), 45 times as large as the cumulative seismic moment (0.04 × 10 16 Nm). The cumulative geodetic/seismic moment ratio reduces over time, but remains larger than 3 during the co-and post-M4.6 stages.

On the Spatial Complexity of Fault Slip Distributions
Fault slip most likely has nonuniform spatial distribution due to spatial heterogeneities of rock strength and stress state on the fault, with well-known dependence on depth and the less understood along-strike variations. Seismic and geodetic inversions can reveal how fault slip is distributed on the discretized fault plane. However, to explore all possible models consistent with observations, the parameter space scales up rapidly to a large number of unknowns, increasing the problem's null-space, which means there are many vectors in the model space that are unconstrained by the data. Therefore, it is reasonable to consider our understanding of the complexity of slip distribution in natural earthquakes. The reasonable approach can allow for fault-slip heterogeneity while keeping the problem null-space as small as possible. Mai and Beroza (2002) compiled published finite-source rupture models and proposed the fractal pattern in slip distributions. It is true for large earthquakes, and multiple fault segments with several rupturing centers are revealed by geodetic and seismological observations, for example, the 2008 M w 7.9 Wenchuan earthquake (Shen et al., 2009), and the 2016 M w 7.8 Kaikoura earthquake (Hamling et al., 2017). However, solving a huge number of parameters has a high computation cost. Computation complexities in their algorithms depend greatly on the number of discretized fault patches. For example, when studying a 40 km-long and 20 km-wide fault with slipBERI, there are 200 patches if the patch size is 2 km and the parameter's dimension is 400. The latter would rapidly increase to 1,600 if the patch size is 1 km. This is possibly the reason why the number of imported fault patches has upper bounds in practice, particularly if a Bayesian sampling strategy is employed. Though techniques like parallel computing have been introduced to improve computation efficiency, sampling such high-dimensional problems is still computationally challenging and does not solve the size of the null-space.
In this research effort, we favored a method that dramatically reduces the number of free parameters to solve; the drawback is that it results in compact fault slip distributions. However, our inverted slip distribution patterns are supported by the observations. This is a reasonable approach because many inversion results support fault-slip distributions that are spatially compact, especially for small-magnitude earthquakes (Ainscoe et al., 2017;Barnhart et al., 2014;Champenois et al., 2017;Taymaz et al., 2007;Xu et al., 2016). Many studies have successfully modeled the majority of surface displacement signals using only one single fault with uniform distribution (Biggs et al., 2006;Nissen et al., 2007;Walters et al., 2009). For SSEs across the global subduction zones, distribution patterns usually follow an elliptical shape with one slipping center (Fukuda, 2018;Villegas-Lanza et al., 2016;Wallace et al., 2012), and the fractal pattern is not required.
Benefiting from the online database of finite fault rupture models, SRCMOD (Mai & Thingbaijam, 2014), we were able to quantitatively evaluate how well a single elliptical model fits the available slip distributions across various tectonic settings and magnitudes. We retrieved 300 slip distributions on a single fault from SRCMOD and intended to model the slip distributions with the one-ellipse model. Our experiments showed that for 85% of M w ≤ 7.5 events, the RMSE of the slip residual is less than 20% of the peak slip ( Figure S10 in Supporting Information S1). In addition, a simple circular crack is also the widely accepted assumed model in stress drop estimation based on seismic spectra (Kaneko & Shearer, 2014;Madariaga, 1976). Though only small degrees of freedom are allowed in the one-ellipse model, complexity could be added by incorporating multiple ruptures. As we showed in Section 2.2, a half-moon pattern was retrieved by two containing or overlapping elliptical crack models. Similarly, it is possible to overlap multiple ruptures to simulate multiple peak slips or more complex patterns.
The compact slip distribution in this new elliptical model is also favorable to evaluate the statistics of small earthquakes. Earthquake source parameters' characterization of small earthquakes is important for understanding the physics of source processes and might be useful for earthquake forecasting (Uchide et al., 2014). A wide-used source model to analyze the source parameters of small earthquakes is a circular crack rupture (Brune, 1970;Madariaga, 1976) with stress singularity at the crack tip, and we hope our new elliptical slip model, which avoids this stress singularly, can be an alternative source model in the future (Shearer et al., 2006). Furthermore, by taking advantage of the improved method for estimating slip rates during temporally overlapping InSAR time frames, one can image the fault behavior over a long period in a relatively high temporal resolution. This new method is expected to be applied to investigate the temporal evolution of slow fault slip, for example, transient slow slip (Khoshmanesh et al., 2015;Klein et al., 2018;Kyriakopoulos et al., 2013), afterslip (Thomas et al., 2014) and SSEs in subduction zones (Bletery & Nocquet, 2020;Ozawa et al., 2019;Rousset et al., 2019).

Time-Dependent Fault Kinematics During Continental Seismic Swarms and Other Slow Earthquakes
During the initial stage of the 2011 Hawthorne seismic swarm, a substantial amount of aseismic slip ruptured on the southern subfault without strong seismicity (e.g., the first two periods in Figure 8b), with peak slip rates of 1.1 ∼ 5.4 cm/day, average slip rate 0.4 ∼ 1.9 cm/day, and migration velocity 0.05 km/day. Note that these values are lower bounds, as the time between two neighboring epochs (Δs n ) of SAR image acquisition time may be longer than the duration of SSEs, preventing capture of short events with higher velocities. The limitation due to the temporal sampling of InSAR could be improved by combining all of the InSAR data sets, or incorporating other high-temporal resolution observations, for example, GNSS or strainmeter observations. We anticipate that the current InSAR temporal sampling limitation will be reduced over the second half of this decade (2020s). Our approach will be well suited to fully utilize the multi-constellation of InSAR capable satellites (Sentinel-1, CosmoSky-Med, PAZ, TerraSAR-X, ALOS-2, and NISAR,etc.). The phenomena potentially driven by aseismic slip are widely explored, for example, ETS, Rapid Tremor Reversals (RTRs), SSEs, fault creep, and fluid injection. To better compare this precursory aseismic slip with other identified phenomena in the slow slip family, we compile the slip rates and migration velocities found in the literature (list below and in Table S1 in Supporting Information S1).
1. The peak slip rate. SSEs show a wide range of peak slip rates among subduction zones, for example, 0.27 cm/ day for the Cascadia subduction zone (Bletery & Nocquet, 2020), 0.3 cm/day for South Central Alaska Megathrust (Rousset et al., 2019), and 0.6 ∼ 2.8 cm/day for Japan trench (Hirose & Obara, 2010;Ozawa et al., 2019). During the early stage of the 2011 Peloponnese seismic swarm (Greece) (Kyriakopoulos et al., 2013), the fault behavior was dominated by aseismic slip inferred from the geodetic and seismic moment, and the peak slip rate was 0.26 cm/day. The maximum slip rate in fault creep events is very low, for example, 0.5 cm/year on the Hayward fault (Schmidt et al., 2005), 0.5 cm/year on the Haiyuan Fault (Jolivet et al., 2012;Song et al., 2019), 0.8 cm/year on the North Anatolia Fault (Hussain et al., 2016), and 3 cm/year on the San Andreas Fault (Johanson & Bürgmann, 2005;Khoshmanesh et al., 2015;Scott et al., 2020). However, in the fluid injection experiment the slow aseismic slip during the early stage was much higher, 4 × 10 −3 mm/s (35 cm/ day) (Guglielmi et al., 2015), potentially because the measurement in the fluid injection is real time, and the duration uncertainty is much lower than SSEs observations. 2. The average rate of slip increment. Research on the 2010-2014 seismic swarm in southern Italy (Cheloni et al., 2017) is consistent with our findings. This research revealed that the average slip rate started to increase 2 months before the largest shock (M w 5.1) and reached the highest value, ∼0.1 cm/day, a few days before the largest shock. It then decreased to zero in the following months. This highest average slip rate was at the same level with ∼0.4-1.9 cm/day in our research. The aseismic slip rate inferred by SSEs is lower, ∼0.03-0.14 cm/day (Radiguet et al., 2011), and this value is much lower inferred by RE, ∼0.3-3 cm/year (Mesimeri & Karakostas, 2018;Nadeau & McEvilly, 1999;Turner et al., 2013). 3. Migration velocity. These velocities of ETS and SSEs vary with subduction zones (Yamashita et al., 2015), but the generally reported migration velocity along the strike of the plate geometry is ∼10 km/day (Wech et al., 2009), while RTRs propagate "backwards" 20-40 times faster than ETS advances (Houston et al., 2011). The large-scale features of ETS propagation with RTRs are reproduced and supported by numerical experiments (Liu et al., 2020;Luo & Liu, 2019). Similarly, migration velocity in TES varies over a wide range, from 0.5 to 14 km/day (De Barros et al., 2020;Passarelli et al., 2018).

Spatially Variable Mechanical Response of the Hawthorne Swarm Faults
As shown in Figure 8b, the southern segment is active during the pre-M4.6 stage, and the fault behavior is mostly dominated by aseismic slip, inferred from a very high geodetic/seismic moment ratio, 45 (Figure 8c), while the general cumulative geodetic/seismic moment ratio remains larger than three for the whole seismic swarm. This significant portion of aseismic slip identified here has been reported to explain the discrepancy between the geodetic moment and the seismic moment in a handful of continental seismic swarms (Cheloni et al., 2017;Gualandi et al., 2017;Kyriakopoulos et al., 2013;Lohman & McGuire, 2007;Wicks et al., 2011). In 2005, a tectonic swarm of over a thousand earthquakes occurred in the Salton Trough, California (USA) and Lohman and McGuire (2007) revealed the geodetic moment of the modeled fault system was about seven times the cumulative seismic moment of the swarm. Wicks et al. (2011) studied a swarm in southeastern Washington (USA) and also found the geodetic/seismic moment ratio was about seven. During the 2011 Peloponnese Peninsula seismic swarm (Greece), Kyriakopoulos et al. (2013) revealed a big discrepancy in moment release, where the geodetic moment was ∼5 times the cumulative seismic moment for the interval July 3-October 1. For the 2013-2014 Northern Apennines seismic swarm (Italy), the moment associated with aseismic deformation/the seismic moment ratio is between 70% ± 29% and 200% ± 70% (Gualandi et al., 2017). For the 2010-2014 Pollino seismic swarm (Italy), Cheloni et al. (2017) found that 70% of the moment was released aseismically. Above all, previous studies require aseismic slip to explain the discrepancy between the geodetic moment and seismic moment for seismic swarms, with the estimated ratio of ∼1.7-7. Furthermore, the compact fault slip identified during the pre-M4.6 stage is favored by our improved methodology as demonstrated in Section 2. The previous finding of fractal distribution of fault slip is based on M5.9+ earthquakes (Mai & Beroza, 2002), while small-to-moderate-magnitude ruptures would have a more compact slip distribution with low complexity as observed in the rupture models SRCMOD (Mai & Thingbaijam, 2014). Therefore, we hope that our improved method can be used to improve the detection of similar small-to-moderate-magnitude aseismic transients in future seismic swarms.
The large disagreement between the geodetic moment and the seismic moment indicates that aseismic slip dominates the fault behavior during the early stage of the 2011 Hawthorne seismic swarm, and seismic slip cannot solely explain the observed surface deformation successfully. Thus, the nucleation of the M4.6 event does not follow the cascading model, which only depends on the stress transfer caused by neighboring foreshocks and aseismic slip is not involved. Here, we test whether the nucleation of the M4.6 event follows another earthquake nucleation hypothesis, the preslip model, where the stress transfer caused by aseismic slip is responsible for the largest shock's occurrence. We utilize the cumulative slip distribution from our inversion model and compute the Coulomb stress change on the fault geometry as shown in Figure 8. The cumulative fault slip caused a Coulomb stress increase over the seismic rupture region of the M4.6 event and the maximum value is 4.1 MPa >0.01 MPa, which is enough to trigger an earthquake (King et al., 1994). In addition, we compute the Coulomb stress change caused by a seismic slip of foreshocks on the fault geometry, and we find that the maximum Coulomb stress increase over the seismic rupture region of the M4.6 event is 1.5 MPa >0.01 MPa, so the cascading model may also play a role in the nucleation process. Note that the stress change analysis based on foreshocks' locations can be affected by many factors, for example, the precision of earthquake hypocenter and the stress drop calculation method. For example, an M w 4.3 foreshock occurred 2 hr before the 1992 M w 6.1 Joshua Tree earthquake, and there are opposite conclusions on whether the mainshock is triggered by the foreshock, by using different spatial resolutions in the foreshock-location-based analysis performed by Dodge et al. (1996) and Mori (1996). Therefore, because the Coulomb stress increase caused by aseismic slip is larger than that caused by seismic slip, 4.1 >1.5 MPa, we interpret that the largest M4.6 event could have been triggered by earthquake nucleation initiated by aseismic slip, but the nearby preceding foreshocks likely also contributed to the nucleation process.
The aseismic slip mainly occurred on the southern subfault during the pre-M4.6 stage, while the most significant seismic slip hit the northern subfault during the co-and post-M4.6 stages. Here, we discuss the possible underlying mechanisms of contrasting behaviors on the two subfaults. One potential cause of the precursory aseismic slip on the southern segment is various dilatancy properties along the strike. Many authors have studied the shear-induced dilatancy, which could increase the effective normal stress and thus favor fault stability (Ciardo & Lecampion, 2019;Segall et al., 2010;Segall & Rice, 1995). For example, to explain abundant microseismicity and aseismic transients in barrier zones on the Gofar transform fault, Liu et al. (2020) proposed a numerical model where strong dilatancy strengthening effectively stabilizes along-strike seismic rupture propagation and results in rupture barriers where aseismic transients arise. If this is also true for the 2011 Hawthorne seismic swarm, the shear-induced dilatancy would explain the aseismic transients on the southern fault and the seismic rupture on the northern subfault. What's more, the requirement of enhanced fluid-filled porosity for the dilatancy strengthening might be filled for the 2011 Hawthorne sequence. The 2011 Hawthorne sequence is close to the Aurora-Bodie volcano (Lange & Carmichael, 1996), and geothermal fluids have been found in this area (Hinz et al., 2010), so it is possible that excess fluids can be persistently supplied and lead to large fluid-filled porosity and high pore pressure. Therefore, the dilatancy strengthening might be one of the underlying mechanics that govern the partitioning between aseismic and seismic slip during the 2011 Hawthorne earthquake swarm.
In addition, the fault geometrical complexity could favor the lateral variation of slip and aseismic slip. First, Romanet et al. (2018) proposed that two overlapping faults can naturally result in a complex seismic cycle without introducing complex frictional heterogeneities on the fault. They found that for two mildly rate-weakening faults with a small distance between the faults, a complex behavior with a mixture of slow and rapid slip can be observed. This finding is consistent with the mixture of slow and fast slip close to the connecting region of two subfaults during the 2011 Hawthorne swarm (triangular subfault in Figure 8). Second, Cattania and Segall (2021) highlight the effect of long-wavelength fault roughness on a range of fault behaviors, foreshocks, and precursory slow slip, during the preparation stage of an energetic event. Their numerical simulation suggested the preparation stage is characterized by feedback between creep and foreshocks: episodic seismic ruptures break neighboring asperity groups and favor the creep acceleration, which loads other asperities leading to further foreshocks consecutively. The coexistence of foreshocks and precursory slow slip, as well as their migration toward the hypocenter of the energetic event in Cattania and Segall (2021), also matched our observation during the pre-4.6 stage ( Figure 8). Therefore, we think fault geometrical complexity might contribute to the precursory slow slip during the 2011 Hawthorne earthquake swarm.

Conclusion
This study developed a new methodology for estimating time-dependent fault slip distributions by incorporating a physics-based crack model as a regularization term. We first introduce two propagation patterns of fault ruptures and then propose a method to solve the complex slip distribution with multiple physics-based crack models. Finally, the performance of the proposed methodology is analyzed with simulated experiments and geodetic observations during a real seismic swarm case. The advantages of the proposed method are as follows: 1. The estimated fault slip solutions describe a compact slip distribution, due to the use of a laboratory-derived crack model. This choice significantly reduces the number of parameters to solve, independently of the subse-quent level of fault discretization. Though the slip complexity is less than in the previous methods, the additional complexity in the slip pattern can be incorporated by incorporating multiple partially or totally overlapping elliptical cracks. 2. The robustness of our method has been analyzed by (a) its capability to reproduce synthetic simulated cases with various slip patterns and by (b) the ability of elliptical slip patterns to reproduce published slip distribution from the SRCMOD database. (c) Our proposed method is applied to estimate a detailed time-dependent fault slip distribution model for the 2011 Hawthorne seismic swarm (Nevada, USA). Our results indicate that the seismic swarm was caused by activity on a two subfault network with different orientations. The results also show that aseismic slip on a southern subfault dominates the fault behavior during a pre-M4.6 stage; after the aseismic pulse (during the most energetic stage), the largest event occurred on a northern subfault. Our results are consistent with an overlapping fault slip migration during the pre-M4.6 stage along the southern fault, followed by larger triggered coseismic ruptures of fault patches along the northern fault. Our model is consistent with small-scale spatially compact fault slip distribution and allows us to estimate lower bounds for the peak and average value of fault slip rates. These lower-bound estimates are consistent with reported values for SSEs and other continental swarms.
The new inversion method presented here is complementary to the existing methodologies to estimate fault-slip distributions using geodetic data. We hope that this approach will be particularly useful with current and near-future multi-constellation InSAR satellite radar interferometry missions. In this near-future context, this tool could improve the identification of similar precursory (aseismic) slow slip during other long-lasting earthquake sequences (swarms) and help understand the driving mechanisms of earthquakes.