Fault Roughness at Seismogenic Depths and Links to Earthquake Behavior

The nonplanarity of fault surfaces, or roughness, may control earthquake behavior such as nucleation, rupture propagation, and slip distribution (e.g., Lindh and Boore, 1981; Okubo and Aki, 1987; Power et al., 1987; Fang and Dunham, 2013). Fang and Dunham (2013) showed through simulations that fault roughness imposed a primary control on local stress heterogeneities as well as the frictional and slip resistance along a fault. Large ruptures tend to start near restraining bends or higher fault complexity that result in a stress asperity (Lindh and Boore, 1981; Goebel et al., 2012; Allam et al., 2019), whereas similar conditions can cause ruptures to stop in regions where stress conditions preclude continued propagation (Fang and Dunham, 2013).

The influence of fault geometry on earthquake behavior may be reflected in the magnitude–frequency distributions of earthquake sequences. The Gutenberg–Richter relationship describes the magnitude (M) distribution of a set of earthquakes (N) and is often formulated as: $log10N∝−bM$⁠. The parameter b, or b‐value, is the slope of the distribution in log‐linear space (Gutenberg and Richter, 1944) and characterizes the relative frequency of larger to smaller quakes. Although b‐values are close to unity on average, variations of b‐values in space and over the earthquake cycle have been widely reported (Schorlemmer et al., 2005; Ogata and Katsura, 2014; Tormann et al., 2014; van der Elst, 2021). Furthermore, laboratory studies have shown that rougher faults may have higher b‐values, whereas lower b‐values were observed near geometric asperities that host large slip events (Goebel et al., 2012, 2015, 2017). Although laboratory studies show a relationship between the b‐value of acoustic emissions and the heterogeneity of the slip surface, it has been considered impossible to observe this correlation in nature (Goebel et al., 2015).

Fault topography or roughness has been inferred from mapping the traces of surface ruptures (e.g., Okubo and Aki, 1987; Wesnousky, 1988) or exhumed fault outcrops (e.g., Power et al., 1987). Recently, fault scanning techniques (light detection and ranging, laser profilometer, etc.) provided high‐resolution images across a wide range of length scales of fault surfaces (e.g., Sagy et al., 2007; Brodsky et al., 2011; Candela et al., 2011, 2012; Parnell‐Turner et al., 2018). These studies confirm the nonplanar nature of fault trends and surfaces. Observations of fault corrugation and geometric anisotropy are common, with faults being somewhat smoother in the direction of slip (Power et al., 1987; Sagy et al., 2007; Candela et al., 2011).

These studies demonstrate that faults deviate significantly from planar surfaces; however, the observations made on exhumed faults may not provide complete representations of active faults at seismogenic depths. For example, imaged fault surfaces may be degraded by weathering (Power et al., 1987). Furthermore, due to changing material composition (e.g., soft sediments, greater heterogeneity) and lower confining pressures near the surface, fault geometries may be more complex near the surface, and become smoother and simpler at greater depths (Sylvester, 1988). Fault slip surfaces are likely to be 3D with anastomosing branches, stepovers, and other complexities that slip simultaneously during rupture (Wesnousky, 1988; Shipton and Cowie, 2001); imaging of surface traces (1D) or exhumed fault outcrops (2D) do not capture the 3D nature of fault surfaces. In fact, Shipton and Cowie (2001) argue that the evolution of the fault geometry through time (i.e., in 4D) is the key to understanding the dynamics of fault slip. Seismic data may provide the best opportunity to map 3D fault roughness, yet imaging active fault surfaces at sufficient resolution has been beyond the capability of available imaging tools or datasets. Most studies of fault structure using aftershocks have been limited to analyzing the more diffuse clouds of seismicity that are often interpreted as damage zones (e.g., Powers and Jordan, 2010), rather than imaging the fault slip surface itself because of insufficient resolution.

Here, we use a prolific earthquake swarm with well‐resolved earthquake locations to probe 3D roughness, for the first time, across an active fault plane at depth and explore its influence on earthquake behavior.

We study an earthquake swarm near Cahuilla, California, that was notable for its productivity, with 18,250 relocated events (−0.27 < M < 4.4), and long duration, lasting approximately 4 yr from early 2016 to late 2019 (Hauksson et al., 2019; Ross et al., 2020). Here, we use the relocated catalog from Ross et al. (2020) (catalog is available from Caltech [2020]). The swarm occurred approximately midway between two major fault systems in southern California—the Elsinore and San Jacinto faults (Fig. 1)—in a region with numerous, albeit generally smaller swarms (Hauksson et al., 2019; Ross and Cochran, 2021). These swarms are spatially dense, long duration, and occur ubiquitously in the region. Their migration patterns suggest that they are likely driven by natural fluid migration through the crust at depths of 5–12 km (Ross et al., 2020; Ross and Cochran, 2021).

The Cahuilla swarm, as identified and precisely relocated (38 and 87 m relative location error in the horizontal and vertical directions, respectively) by Ross et al. (2020), begins at a depth of ∼7 km, and is well fit by a plane striking 343° and dipping −82°, matching the right‐lateral, strike‐slip focal mechanism of the largest event (M 4.4) of the sequence. The Cahuilla swarm is inferred to occur on reactivated joints or foliations with the minimal total slip (Hauksson et al., 2019), so should be considered an immature fault. Earthquakes are distributed across much of the fault surface, although with somewhat higher densities in a 500 m wide zone extending from the presumed injection point up dip to the top edge of the swarm (Fig. 2; Ross et al., 2020). The event densities suggest channeling of events along strike away from the high‐density “pipe” of events extending up dip from the injection point, similar to channeling inferred from depth histograms by Ross et al. (2020). We use the 10th‐percentile rupture time in a 150 m × 150 m grid (i.e., the 10th percentile of the range of earthquake origin times in each grid cell) across the fault plane to show the migration of events (Fig. 2). Over a period of 4 yr, the earthquakes migrate approximately 1 km down dip and 3 km up dip, as well as approximately 2 km bilaterally along strike. Events initially migrate steadily for over an ∼2 yr period along strike and dip. However, once earthquakes reach an apparent permeability barrier located approximately 1.75 km up dip from the injection point, migration slows. Then, following the M 4.4 mainshock (day 957), rapid migration takes place across the ∼1 km wide and 4 km long upper section of the fault.

The M 4.4 earthquake occurred following the breach of the inferred permeability barrier (Fig. 2; Ross et al., 2020). Interestingly, stress drops of events above the barrier were reported to be distinctly lower than for the rest of the sequence, perhaps suggesting a difference in the material properties across the barrier (Ross et al., 2020). Waveforms from the M 4.4 earthquake were too complex to relocate with waveform cross‐correlation techniques, so we approximate the M 4.4 rupture area relative to the relocated catalog using the distribution of earthquakes that occur one day after the M 4.4 (Fig. 1). The estimated rupture area of the M 4.4 matches the size of the fault area expected to slip assuming a circular source with a stress drop of 8 MPa—the median value reported for the sequence (Ross et al., 2020).

We define fault roughness to be the deviation of earthquake locations from a linear or planar fault surface following Malinverno (1990). We first estimate a 2D fault roughness in the along‐strike and along‐dip directions. We assume a general fault orientation identical to that of the M 4.4 earthquake (strike 17°, dip 82°). The rake of the M 4.4 is 179°; therefore, we assume that slip parallel and orthogonal to be in the direction of strike and dip of the fault orientation. We divide the fault into a set of along‐strike and along‐dip profiles using nonoverlapping bins every 150 m. We require at least 250 earthquakes to retain the profile. The 2D roughness is estimated as the mean distance to the best‐fit line given by $1n∑i|axi+byi+c|a2+b2$⁠, for points (⁠$xi$⁠, $yi$⁠) and a line defined by $ax+by+c=0$⁠.

We also estimate roughness in 3D at all earthquake locations with at least 100 neighboring earthquakes within 500 m. For each roughness estimate, we calculate the best‐fit plane using principal component analysis for all points (⁠$xi$⁠, $yi$⁠, $zi$⁠) within 500 m of the earthquake. The 3D roughness is then estimated as the mean out‐of‐plane distance to this best‐fit plane, given by $1n∑i|axi+byi+czi+d|a2+b2+c2$⁠, for a plane defined by $ax+by+cz+d=0$⁠. The method provides a measure of roughness at a single‐length scale (i.e., 1 km). We also estimate roughness at multiple length scales, using the same method but considering seismicity at different radii from an earthquake, to measure the roughness scaling exponent. The results of the simple roughness‐length method of Malinverno (1990) should be interpreted with some caution, because it subsumes several types of fault complexity including point cloud diffusivity, degree of branching or anastomosing, and fault bending. In particular, for seismic data it should only be interpreted as fault roughness if near‐contiguous fault surfaces are apparent in the data; otherwise, the metric would only measure how diffuse are the earthquake locations.

We estimate b‐values using the “b‐positive” method of van der Elst (2021), which uses the positive magnitude differences between successive earthquakes and is robust to transient changes in the completeness. We estimate b‐values at the same points in which 3D roughness is calculated, but for this calculation we use the nearest 150 M ≥ 0.6 neighboring earthquakes with positive magnitude differences. Thus, b‐values are determined over radii from 400 to 800 m around each point with a median value of 540 m, depending on the density of earthquakes.

The along‐strike and along‐dip profiles delineate a clear set of fault surfaces (Fig. 3a). The along‐dip profiles show bends, anastomosing branches, and stepovers, especially toward the southeast portion of the rupture. In the along‐dip profiles, we observe extensive branching and a bend in the fault surface near the location of the inferred permeability barrier (Fig. 3a). The bend is the most prominent along the northwest side of the fault and has an amplitude of ∼500 m out of plane from the remainder of the profile. Near the rupture area of the M 4.4 earthquake on the southeast side of the fault, we observe discontinuities in the fault surface due to both branches and stepovers. We measure variable roughness across the fault profiles with mean out‐of‐profile distances between 29 and 127 m along the ∼4 km profiles (Fig. 3a).

The along‐strike profiles are not only more similar to each other in appearance and have a smaller range of mean out‐of‐plane distances (27–90 m), but also show evidence of complexity. The typical amplitudes of these features in the both along‐strike and along‐dip profiles are generally less than 200 m, with wavelengths (i.e., along‐profile lengths) of about 500 m to 2 km. Overall, the along‐dip profiles are 50% rougher than the along‐strike profiles, with an average out‐of‐plane distance of 62 m for the along‐dip profiles compared with 43 m for along‐strike profiles (Fig. 3b). Thus, this strike‐slip fault is smoother in the direction parallel to slip.

We next extend the analysis to measure complexity across the 3D fault surface. Figure 4a shows the mean out‐of‐plane distance, or 3D roughness, that is, the root mean square (rms) residuals to a plane fit to all events within 500 m of an individual earthquake location. We observe that 3D roughness varies by a factor of 8 across the fault surface (0.01–0.08 km). The highest roughness values (∼80 m) at this length scale are found within the estimated rupture area of the M 4.4 earthquake, where the 2D profiles show significant fault complexity. Similar to the 2D profiles, the 3D roughness also shows fault corrugation or broadly varying patterns of fault roughness. The corrugation is oriented subparallel to oblique to the strike (or rake) direction such that it is somewhat inclined relative to strike from southeast to northwest. The range of roughness values does not systematically change with depth or along strike.

We similarly explore the distribution of b‐values across the fault (Fig. 4b). We require 150 positive magnitude differences and impose a magnitude of completeness cutoff (⁠$Mc$ 0.6) to estimate b‐values; therefore, b‐value estimates sample a range of radii around the earthquake of interest, typically extending 400–800 m around each point (the median radius is 540 m). The b‐values vary by a factor of ∼2 (0.8–1.8) across the fault. The data suggest similar, albeit somewhat weaker, corrugation in b‐values as were observed in the 3D roughness. Qualitatively comparing Figure 4a and 4b, we find that areas of the fault with higher roughness tend to have higher b‐values. The exception is in the M 4.4 rupture area, where we observe the highest roughness values with significant fault branching and bending, but the corresponding b‐values are around 1.

We quantitatively compare 3D roughness and b‐values for the whole sequence by plotting roughness and b‐value estimates for each earthquake. We compare roughness and b‐values before and after the M 4.4 in Figure 4c,d. Roughness values pre‐M 4.4 are between 10 and 50 m (Fig. 4c), whereas post‐M 4.4 roughness values are larger (15–110 m) (Fig. 4d). The differences in roughness that we report for events before versus after the M 4.4 mainshock primarily reflect changes in the spatial sampling of roughness (i.e., events occurring in new locations), rather than temporal changes in roughness. Pre‐ and post‐M 4.4 b‐values are similar, but post‐M 4.4 b‐values have slightly higher maximum b‐values (from ∼1.5 to ∼1.65). Although a weak correlation (correlation coefficient of 0.34) exists between b‐value and roughness for the whole sequence (Fig. S1, available in the supplemental material to this article), post‐M 4.4 we observe an area of high roughness values with corresponding lower b‐values (∼1.0–1.2); these anomalous values are primarily within the estimated rupture area of the M 4.4 earthquake.

We calculate errors for our roughness and b‐value estimates with bootstrapping (resampling the points with replacement, and in the case of the roughness calculation, refitting the plane). About 95% confidence bounds on our estimates are shown (for selected points, given that they are highly spatially correlated and overlap) in Figure S1. This plot shows that b‐value and roughness estimates in different areas of the fault are distinct, given errors due to sample size.

Next, we examine how roughness evolves in space as events migrate across the fault surface by plotting the roughness values for the whole sequence at the time the event occurred (Fig. 4e). We observe a broadening in the range of mean out‐of‐plane distances through time as more of the fault surface is sampled. Early in the sequence (<600 days) roughness values are 20–35 m, expanding to 10–60 m prior to the M 4.4 earthquake as the fluid migration continues, and we can measure roughness values across a larger fault section. A clear increase in the maximum observed 3D roughness values is observed after the M 4.4 earthquake as events expand across new sections of the fault. This increase in mean out‐of‐plane distance values is due to the high roughness measurements within the inferred rupture area of the M 4.4 reflecting the complex fault structure that includes a number of bends and stepovers (see Fig. 3a). We are unable to determine temporal changes in roughness across the fault because sections of the fault are active during different subsets of the study period due to spatial migration of the earthquakes.

The distribution of event densities across the fault may suggest fluid channeling (Fig. 2a) that could be caused, at least in part, by variations in fault roughness. The orientation and scale of the channeling inferred from the event densities (Fig. 2) are similar to the longer wavelength spatial variability in 3D roughness (Fig. 4). However, correlation values for event density and roughness (or b‐value) are low. This is likely because the total number of events along a section of fault is primarily controlled by the pore pressure changes due to the natural injection, with the highest event densities in an approximately 1 km wide vertical section above the inferred natural injection point.

We repeat the 3D roughness calculation for a range of radii (from 125 m to 1.5 km) to sample the roughness across different length scales. Figure 5 shows the distribution of roughness measurements at different length scales (i.e., twice the radii) as well as their mean and median values. We find that the roughness measurements across length scales display power law behavior consistent with a self‐affine surface; that is, the surface is self‐similar with different scaling factors in the along‐strike and along‐dip directions. By fitting a line to the roughness and length scale estimates, we determine how the roughness of the fault surface changes with surface length, with the slope of the line referred to as the Hurst exponent (⁠$ζ$⁠) (e.g., Fardin et al., 2001; Beeler, 2021). We find a slope (⁠$ζ$⁠) of 0.52 across length scales from 250 m to 3 km consistent with a Brownian surface.

We explored fault roughness and earthquake behavior across a fault illuminated by a long‐duration earthquake swarm near Cahuilla, California. The swarm occurred on a 4 km × 4 km immature fault at 4–8 km depth in low‐permeability plutonic rocks (Hauksson et al., 2019). The dense, low‐magnitude seismicity and high‐resolution locations allow us to map stepovers, fault branching, and corrugation. Corrugation of the fault zone is subparallel to the strike direction with a wavelength of approximately 1–2 km and amplitudes typically of a few hundred meters but as large as 500 m. Our observations are similar to those of John (1987) who examined an exposed normal fault system and found corrugations with wavelengths between 0.2 and 10 km and amplitudes from 30 to 400 m. Examining a smaller fault surface, Sagy and Brodsky (2009) found a broadly undulating fault surface with small (10–40 m), quasi‐elliptical bumps protruding ∼1 m out of the surface. Our data are not of sufficiently high resolution (relative location errors in the tens of meter) to show meter‐scale resolution of the fault surface, but the amplitude to wavelength ratios are similar between our results and the previous studies.

We also estimate the Hurst exponent (⁠$ζ$⁠) to quantify how roughness changes with scale (Candela et al., 2011; Beeler, 2021). The previous studies found fault profiles or surfaces have Hurst exponents ranging between Brownian and self‐similar (⁠$0.5≤ζ≤1$⁠) (Candela et al., 2012; Beeler, 2021). Here, we find $ζ=0.52$ consistent with a Brownian surface. This can be interpreted as the surface being somewhat more correlated over short distances than long distances. The Hurst exponent found here is at the low end of what is typically seen in natural faults, such that the out‐of‐plane amplitudes of closely located points on a fault are somewhat less correlated than is generally found along natural faults and lower than for a self‐similar surface. The relatively low $ζ$ determined here may reflect the immature fault surface with lower total slip (assuming slip acts to smooth the fault surface), or that we are able to better measure the 3D nature of the fault surface that includes fault branches and stepovers. Overall, we find that fault roughness at seismogenic depths is consistent with measurements made using different techniques on mapped surface ruptures and exhumed faults.

The distribution of b‐values across the fault follows similar spatial patterns as fault roughness. We observe no clear depth correlation with b‐value over the depth range of the sequence (4–8 km). We do find a weak, positive correlation between b‐value and roughness with a cross correlation value of 0.34. In laboratory studies, Goebel et al. (2017) suggested that b‐values are typically higher on rougher fault surfaces. Simulations have shown that the minimum and the maximum magnitudes depend on small‐scale fault roughness, with shorter wavelength roughness associated with smaller and more numerous earthquakes (Heimisson, 2020). The fault roughness may control the size distribution of earthquakes such that higher roughness is associated with higher b‐values, that is, a greater proportion of smaller events. Given the weak correlation observed here, it may be useful to examine b‐value and roughness using additional earthquake datasets and simulations.

We find that the largest earthquake (M 4.4) occurs in a region with relatively low b‐values (1.0–1.2) and corresponding anomalously high roughness (50–80 m) estimates. The correlation between roughness and b‐value may break down at higher stressing rates or near asperities where larger events are more likely to occur (Goebel et al., 2012, 2015, 2017). Laboratory studies suggest low b‐values correspond to regions where large slip events occur (Goebel et al., 2012, 2015), in agreement with our findings. Larger events observed in the field, laboratory, and simulations are found to preferentially initiate and terminate near fault bends or heterogeneities (Lindh and Boore, 1981; Goebel et al., 2012; Allam et al., 2019). Furthermore, studies have shown that roughness controls the background stress on faults (Fang and Dunham, 2013), nucleation patch size (Okubo and Dieterich, 1984), and the maximum magnitude (McLaskey and Lockner, 2014). In the Cahuilla swarm, the high roughness measurement near the rupture area of the M 4.4 reflects multiple branches and stepovers, and such geometrical heterogeneities may be associated with larger stress (Scholz, 1968; Fang and Dunham, 2013).

We show for the first time that roughness can be measured using a high‐resolution catalog of dense earthquake locations along an active fault. Earthquake simulations and laboratory studies show fault roughness controls aspects of earthquake sequences and rupture processes. However, earthquake simulations currently use generic, randomly generated roughness distributions to develop various rupture scenarios. Using high‐resolution catalogs of small‐magnitude events to estimate the Hurst exponent and to image fault features including stepovers, branches, and bends, we could provide bespoke geometries and roughness scaling for faults as input into simulations. This could potentially lead to improved forecasting of where earthquakes might start or stop, their slip distributions, and other information important to understanding fault‐specific hazard. However, the simple metric (i.e., mean rms residuals to a profile or plane) that we used here to estimate roughness should be applied only to datasets that contain sufficiently dense and well‐located seismicity to image fault slip surfaces.

We measure roughness of an active, strike‐slip fault at depth using earthquake locations from a prolific, multiyear earthquake swarm. We find that the scaling of fault roughness is self‐affine, and we estimate a scaling exponent (0.52) that is consistent with a Brownian surface. Furthermore, our at‐depth roughness measurements are consistent with other measurements of fault roughness estimated from exhumed faults. We find that the strike‐slip fault is approximately 50% rougher in the along‐dip direction than in the along‐strike direction (i.e., slip direction), which is consistent with past observations and the intuition that faults are smoothed in the direction of repeated slip over time. Finally, we find some evidence for a weak correlation between fault roughness and b‐value across much of the fault; however, analyses of data from other seismically active faults would help to confirm this observation.

The supplemental material to this article includes a figure of b‐value and roughness for the entire sequence with 95% confidence ranges for selected points. The Cahuilla swarm catalog is publicly available from the Southern California Earthquake Data Center (Caltech, 2020). All waveform, parametric data, and the conventional catalog are available from the Caltech and U.S. Geological Survey Southern California Seismic Network (Caltech and U.S. Geological Survey, 2023b) and at the Southern California Earthquake Data Center (Caltech and U.S. Geological Survey, 2023a).

The authors declare no competing interests.

This work benefited from discussions with Annemarie Baltay, Nick Beeler, Peter Bird, Emily Brodsky, Tom Hanks, Steve Hickman, and Sarah Minson. The authors thank Ole Kaven, Greg McLaskey, an anonymous reviewer, and a TSR associate editor for their suggestions that greatly improved this article. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.