Resolving Fine-Scale Heterogeneity of Co-seismic Slip and the Relation to Fault Structure

Fault slip distributions provide important insight into the earthquake process. We analyze high-resolution along-strike co-seismic slip profiles of the 1992 Mw = 7.3 Landers and 1999 Mw = 7.1 Hector Mine earthquakes, finding a spatial correlation between fluctuations of the slip distribution and geometrical fault structure. Using a spectral analysis, we demonstrate that the observed variation of co-seismic slip is neither random nor artificial, but self-affine fractal and rougher for Landers. We show that the wavelength and amplitude of slip variability correlates to the spatial distribution of fault geometrical complexity, explaining why Hector Mine has a smoother slip distribution as it occurred on a geometrically simpler fault system. We propose as a physical explanation that fault complexity induces a heterogeneous stress state that in turn controls co-seismic slip. Our observations detail the fundamental relationship between fault structure and earthquake rupture behavior, allowing for modeling of realistic slip profiles for use in seismic hazard assessment and paleoseismology studies.


Validation of displacement measurements using SPOT data
To validate our displacement measurements derived from the 1 m National Aerial Photography Program (NAPP) air photos we used an independent dataset that can also measure the deformation field and slip distributions of both the Landers and Hector Mine earthquakes.
Our independent data consists of pre-and post-event 10 m SPOT imagery that cover both earthquakes in space and time, giving an independent constraint on the surface deformation field and slip distribution, which is produced from the same COSI-Corr procedure as that used for the air photos.
The Landers SPOT correlation ( Supplementary Fig. S1 a,b) is produced from a pre-event, 1991, 10 m SPOT 2 image, and a post-event, 1993, 10 m SPOT 2 image. For Hector Mine we produced the correlation maps ( Supplementary Fig. S1 d,e) from pre-event 1993, 10 m SPOT 4 imagery and post-event imagery from a 2000, 10 m SPOT 4 image. For both earthquakes we used a multi-scale correlation window of initial 64 pixels and final of 32 pixels with a step size of 8 pixels, similar correlation parameter as that used for the air photo data, and applied a nonlocal means filter to the correlation result to suppress noise. To measure the displacement we used stacked profiles of 2.5km width, spaced with the same distance along the surface rupture.
Comparing the slip distribution from the SPOT correlation to the air photos we find an excellent agreement at the first-order scale ( Supplementary Fig. S1). The overall smoothness of the SPOT result is a function of a (i) coarser pixel resolution (10 m of SPOT data compared to 1 m air photos), (ii use of a non-local means filter to reduce noise in the correlation result and (iii) use of wide stacked profiles (2.5 km compared to 138 m in the air photos) with also the same coarse measurement discretization. Therefore the excellent agreement between these two different geodetic datasets at the first-order scale, which is what the SPOT imagery can resolve, provides a robust validation of our displacement measurements.

Synthetic tests
To quantify the measurement precision and any possible bias in estimating displacement we employed synthetic tests. Specifically, we design these tests to capture and quantify how the estimation of displacement measurement is effected by (i) noise and geometrical artifacts within the correlation maps and (ii) the subjective manner when interpreting the deformation signal within the stacked profiles. To do this we create a synthetic rupture through the air photos with pre-determined displacement that is constant along the fault, which is then measured using stacked profiles in order to quantify how well we can recover the spatially uniform values. Thus when measured if any spatial variation of displacement does arise that deviates from the known constant value, this directly constrains the effect of noise, long-wavelength geometrical artifacts (such as scanning, thermo-mechanical warping or radial distortion) from the correlation maps and the subjective interpretation of estimating displacement.
To simulate a synthetic rupture on a 'fault' we simply shift one half of a post-event, aerial image by a pre-determined amount. To generate a displacement for the synthetic fault rupture without it being known to the user when measuring displacement we use MATLAB's random number generator, where the value is recorded and only revealed until after the measurements are complete, therefore allowing us to conduct measurements of displacement without knowing the true value and biasing the results. We correlate the artificially dextrally sheared image with a different non-disturbed pre-event image ( Supplementary Fig. S2) using the same processing procedure in COSI-Corr that we use to produce our real dataset. We measure displacement using stacked profiles of 138 m width and lengths up to a maximum of 3 km. Where multiple fault strands exist within a single profile we measure displacement from each fault independently and simply sum these to give the total displacement accommodated across the system of faults.As the images are aquired over areas that contain real tectonic deformation, which are areas we want to neglect from our tests (as this would contaminate our synthetic known ruptures), we simply used information from the real correlation result (Fig. 1 of the main text) to help avoid such areas. To shift image pixels by a non-integer amount, we applied a cubic interpolation to one side of the deformed image. In total we used five image pairs, one pair from the Landers (Supplementary As the synthetic offset is determined to be constant along the entire fault length, any deviation of the measurements from the synthetic, known values gives a direct quantification of the measurement bias and precision, which reflects measurement subjectivity and variation in image quality and texture. Furthermore repeatedly measuring the displacement along-strike of our synthetic fault across different parts of the correlation image gives multiple independent measurements of displacement. We estimate measurement bias by simply subtracting the average displacement (from a large sample) that we measure from the correlation maps, to the true, known synthetic displacement value. To avoid biasing the measurements, the true, predetermined synthetic fault offset values were chosen using a random number generator and were not revealed to the user until the measurements from the correlation maps were completed.
In total we performed three types of synthetic tests that utilized an array of faulting styles, widths and displacements expected in a real environment from imagery acquired over both the Landers and Hector Mine events, that are outlined in detail below ( Supplementary Fig. S2).

Results
The first test involved imagery associated with the Landers event with the aim of understanding whether COSI-Corr's sub-pixel correlation behaves differently when evaluating fault displacement that exceeds a pixel dimension i.e., 'supra' pixel fault movement. To test whether COSI-Corr gives consistent results and does not artificially alter displacement if the fault movement exceeds a pixel dimension, we used one fault with supra-pixel movement of 2.5 m and a second with sub-pixel movement of 0.4 m displacement (where the image pixel resolution we use is 1 m). We used fault lengths of ~7 km that span the entire image and aquired 45 displacement measruements for each fault. We note this first test is similar to that conducted by ref. 40 found in the supplementary file, and we re-performed this test in order to establish consistentcy between the two studies.
The first synthetic test reveals an overall bias of 0.02 m with a precision of + 0.04 m (1σ) for our displacement measurements ( Supplementary Fig. S3). Specifically, both the sub-and supra-pixel displacement tests yield the same behavior with a 1σ measurement precision of + 0.041 m and + 0.038 m, respectively, in agreement with previous studies 1,2 .
The second and third tests involved imagery used for the correlation of the 1999 Hector Mine event, again, from locations that contain no real tectnoic deformation. Specifically, the second test invovles performing additional displacement measurements, similar to the first test, in order to understand whether the measurement uncertainity is different from the correlation produced by air photos associated with the Hector Mine earthquake in comparison to air photos used for the Landers event (even though the post-Landers and pre-Hector Mine are derived from the same 1994, NAPP flight mission). For test 2 we implemented six faults of equal length and spacing but with different magnitudes of displacement, that were set constant along-strike (Supplementary Fig. S2 d,e,f). In total we collected 434 displacement measurements, finding an overall bias of 0.01 m and precision (1σ) of + 0.06 m ( Supplementary Fig.S 4c). Importantly, we found no discernable diffference of the measurement precision or accuracy derived from the Landers air photos (test 1, with a bias of 0.02 m and 1σ of + 0.04 m) or the Hector Mine imagery (see Table S1 and S2 and Supplementary Fig. S3 and S4).
The third test involved changing the width of deformation to understand whether a wider fault zone causes the measurement estimation to be more amibigous and therefore more difficult to interpret and adequatley extract the true displacement. To simulate 'faults' of various widths we simply varied the number of multiple parallel fault strands within a synthetic 'fault zone'. We simulated six different main 'fault zones' with widths of 1, 25, 50, 100, 125 and 150 m, that are of constant width and displacment along-strike. We simulate widths that are similar to those observed in the real correlation results and that observed from 1,3,4 . From this third test we collected 435 displacement measurements and found an overall measurement bias of 0.02 m and precision of 1σ + 0.06 m ( Supplementary Fig. S5). Importantly, we found the displacement bias and precision is not significantly affected by the width of deformation ( Supplementary Fig. S6).
The similar level of bias and precision found in test 3 to that of test 1 and 2, again indicates the width of deformation has an unnoticeable effect on the estimation of displacement in the stacked profiles.
For all three tests we found the spectrum of these results is close to but does not follow white noise ( Supplementary Fig. S7), in agreement with a previous and similar analysis 2 and importantly, nor do they have the same slope in the power spectrum as the real displacement measurements reported in the main text. The slope of the slip spectra for tests 2 and 3 are 0.31 (corresponding to a fractal dimension [D] = 2.34) and 0.47 (D = 2.26), respectively (Supplementary Fig. S7 b and d). Furthermore, we note that both the slope and the power of the synthetic slip in the spectra are significantly diferrent (and less in power) than that found in the real results presented in the main text, again indicating variation of slip induced by measurement subjectivity and noise within the images, is insufficient in amplitude across a broad range of frequencies to cause the observed amplitude of slip variaiton that is observed in our real correlation maps (Fig. 1).
From all three synthetic tests, we found a consistent value for the measurement bias and precision of displacement when estimated from the stacked profiles. From a total of 959 displacement measurements from the three synthetic tests, we found an overall measurement bias of 0.02 m and precision (1σ) of + 0.06 m with no observable systematic variation as a function of the air photos used (Landers or Hector Mine flight missions), fault width or magnitude of displacement. Using these results we derived an empirical error distribution for the displacement measurements ( Supplementary Fig. S7c), which we propagate through the pre-existing error in COSI-Corr, that forms the basis of the measurement uncertainty found in the main text. The small precision and measurement bias is reflective of the robustness of the correlation procedure, the quality of the aerial images, and the use of profile stacking to suppress noise. These tests also confirm the width of our stacked profiles is appropriate, a width that reduces noise to a low level of uncertainty (1σ = 0.06 m), the same as that observed from ref. 17, that found using the same Landers images used here, this value of uncertainty is the base level of noise in the correlation result between the before and after images. Thus, these results and their agreement with previous work indicate a stack profile of 138 m width is sufficient to reduce noise to its floor level, but not unnecessarily wide to 'over-smooth' displacment along the surface rupture. We note that Figs. For a total of 20,000 slip distributions (10,000 for each earthquake), we found that the estimate of the fractal dimension follows a Gaussian distribution with a 1σ of + 0.02 and + 0.03 for Landers and Hector Mine, respectively ( Supplementary Fig. S9). Using a T-test we compared the statistical significance of the difference between the fractal dimension of the Landers (1.72 + 0.02) and Hector Mine slip distribution (1.62 + 0.03) using the values obtained from the Monte Carlo simulations. We found the T-test rejects the null hypothesis that the fractal dimension of the Landers and Hector Mine slip distributions are drawn from the same distribution, with a pvalue of 0.0014, well below our confidence level of 5%.

Boxcounting of the Landers and Hector Mine fault systems
Fault systems and fractures have been shown to follow fractal distributions and previous studies have characterized the fractal properties of faults using boxcounting methods 41,42 . Here we used a boxcounting method to estimate D in order to quantify which fault system, Landers or Hector Mine is more complex and whether the geometrical complexities that compose these faults systems can be considered scale invariant and fractal structures. The boxcounting method allows one to draw log-log plots of the number of boxes required to cover the object (N r ), in this case the mapped fault traces, against the size of the box (r). The slope of the log-log plot provides an estimation of D, that is between 1 and 2 in this case.
To be considered a fractal object the boxcounting curve needs to follow a straight line in log space, defined by the following relation: log(N r ) = a + D * log(1/r) 5 . To estimate the fractal dimension we used boxes of decreasing size and count the number of them that cover the mapped surface fault traces ( Supplementary Fig. S10). We use boxes of sufficient size that adequately covers a range of spatial scales similar to that of the distance between individual fault strands.
The fault traces were acquired from the USGS (http://earthquake.usgs.gov/hazards/qfaults/google.php), and were produced by careful mapping immediately after both earthquakes 3,8 . Therefore importantly the fault traces analyzed here represent the faulting directly involved in the rupture and are not inactive and therefore irrelevant structures.
To minimize any possible distortions or avoid introducing any potential artifacts, we take a single high resolution image that includes both the fault systems of Hector Mine and Landers, as they are fortuitously only 20 km apart. Therefore both fault traces are projected in the same map projection system (UTM WGS-84, Zone 11N) causing minimal distortion to the geometry of the fault traces and imaged at the same resolution (700 dots per inch) and scale (1:400,000).We then split the image into two in order to appropriately separate the Landers and Hector Mine fault systems and then apply the boxcounting analysis separately to each image ( Supplementary Fig. S10). Thus from this procedure we can be confident any difference in the estimation of the fractal D from the boxcounting method is primarily a reflection of fault geometrical complexity rather than a reflection of any spurious artifacts.

Results
We analyze the boxcounting curve over more than 2 orders of magnitude (10 2 -10 4 pixels). From a linear regression to the boxcounting curve we find a fractal dimension of 1.29 + 0.01 with R 2 = 0.98 for the Landers fault system and 1.15 + 0.01 with R 2 = 0.99 for the Hector Mine rupture, which is seen in Supplementary Fig. S11 as a steeper curve. Thus, the boxcounting method clearly indicates the Landers surface rupture has a more complex geometry than the Hector Mine rupture, which can also be easily validated with a simple visual comparison of the two surface fault traces of the two earthquakes ( Supplementary Fig. S10). We also tested whether the initial box size affects the estimation of the fractal dimension by rotating the faults ( Supplementary Fig. S10). We find an insignificant difference when varying the initial box size  Table S2 for statistics for each fault.  Measurements of displacement for the 1992 M w 7.3 Landers earthquake. Column 1 and 2 contains gives the location of the displacement measurements from the stacked profiles in Latitude and Longitude, column 3 details the displacement and column 4 the 1σ uncertainty of the displacement measurement.