Distribution of Aseismic Deformation Along the Central San Andreas and Calaveras Faults From Differencing Repeat Airborne Lidar

Fault creep reduces seismic hazard and serves as a window into plate boundary processes; however, creep rates are typically constrained with sparse measurements. We use differential lidar topography (11–13 year time span) to measure a spatially dense surface deformation field along a 150 km section of the Central San Andreas and Calaveras faults. We use an optimized windowed‐iterative‐closest‐point approach to resolve independent creep rates every 400 m at 1–2 km apertures. Rates vary from <10 mm/year along the creeping fault ends to over 30 mm/year along much of the central 100 km of the fault. Creep rates are 3–8 mm/year higher than most rates from alignment arrays and creepmeters, likely due to the larger aperture of the topographic differencing. Creep is often focused along discrete fault traces, but strain is sometimes distributed in areas of complex fault geometry, such as Mustang Ridge. Our observations constrain shallow seismic moment accumulation and the location of the creeping fault trace.


Introduction
Some faults in the San Andreas fault system release a portion of their accumulated seismic moment aseismically via fault slip, lowering the frequency and/or magnitude of damaging earthquakes (Peng & Gomberg, 2010;Steinbrugge et al., 1960). Characterizing fault creep is a vital part of seismic hazard assessment. The San Andreas fault is the main tectonic boundary between the Pacific and North American plates, and the Pacific-Sierra Nevada Great Valley boundary relative velocity is 39 ± 0.2 mm/year (Argus & Gordon, 2001). The 140 km-long Central Creeping Section of the San Andreas fault (CSAF) and adjacent 40 km of the southern Calaveras fault (CF) reach from San Juan Baptista and Hollister, CA to Parkfield, CA, and are bounded by locked fault segments that ruptured in large earthquakes in 1906 (northwest) and 1857 (southeast; Figure 1). The maximum alignment array creep rate along the CSAF is~33 mm/year (Burford & Harsh, 1980). The CSAF releases most seismic moment as aseismic creep, but questions remain regarding the possibility that a rupture could propagate through the CSAF and the frequency of moderatemagnitude earthquakes (Chen & Burgmann, 2017;Field et al., 2013;Harris, 2017;Noda & Lapusta, 2013;Schwartz, 2018).
We measure deformation along the CSAF and CF from differential airborne lidar topography acquired between 2005 and 2018 (Figures 1 and 2). Having preexisting data before a surface rupturing earthquake and documenting fault creep were primary justifications for collecting the 2005 B4 and 2007 EarthScope data sets (Bevis et al., 2005;Prentice et al., 2009). In addition to the synoptic analysis, we compare the strain field to tectonic landforms at Mustang Ridge ( Figure 1). The lidar's 1-2 km aperture is at an intermediate scale to near-field creepmeters and alignment arrays and far-field InSAR and Global Navigational Satellite System (GNSS) data ( Figure 3). Differencing results are sensitive to tectonic and nontectonic landscape changes and data set errors. The net fault offset from creep in this study (~0.3 m) is lower than the 1-2 m offset in the~M7 earthquakes previously studied with differencing but is larger than typical lidar error (~0.05-0.15 m), demonstrating the sensitivity of this approach.

Lidar Data
We use four airborne laser scanning data sets: (1) 2005 B4 (wavelengh = 1,064 nm; Bevis et al., 2005Bevis et al., , 2017 data with~2 ground pt per m 2 , (2) 2007 EarthScope (1,047 nm; EarthScope, 2008) with~2.5 ground pt per Colors indicate the right-lateral creep rate derived from lidar topographic differencing in this study. Black lines are from the USGS/CGS Quaternary Faults and Folds database (https://www.usgs.gov/natural-hazards/earthquake-hazards/ faults) (U.S. Geological Survey, n.d.). The 1857 and 1906 labels show rupture locations of large SAF earthquakes in those years. Blue dots show M5.7-6.5 earthquake locations from 1882-1885 (Toppozada et al., 2002). m 2 , (3) 2018 Parkfield with~9 ground pt per m 2 (532-1,550 nm; Vadman, 2019), and (4) 2018 Salinas with 4 ground pt per m 2 (1,064 nm; U.S. Geological Survey, 2018), as shown in Figure S7 in the supporting information. The flight lines are all parallel or subparallel to the CSAF. The National Center for Airborne Laser Mapping collected the B4, EarthScope, and Parkfield data and performed ground classification using TerraScan software. The U.S. Geological Survey contracted the FEMA-supported Salinas data collection, which were processed and classified by the vendor (Woolpert). We transformed the point clouds to UTM Zone 10N NAD83 coordinates with ellipsoid heights.

Surface Deformation
We use a windowed implementation of the point-to-plane Iterative Closest Point (ICP) algorithm to calculate the 3-D displacement field over the 11-13 year lidar timespan (Besl & McKay, 1992;Chen & Medioni, 1992). Prior to the ICP analysis, we remove points where elevation changes between the two data sets exceed 40 mm/year, determined by subtracting 1 m resolution gridded elevation data (lower panel in Figures S8-S21). These areas are indicative of landslides (prominent in the NE rate fields in Figures S14-S15), earthflows, fluvial system, or anthropogenic landscape changes that complicate detecting tectonic deformation. For the displacement rate and slip calculations, we rotate our coordinate system 45°c lockwise to be aligned with the average CSAF trend. The ICP algorithm iteratively aligns windows of both data sets using a rigid deformation (translation and rotation; scale is constant). We use the LIBICP (Library for ICP) algorithm (Geiger et al., 2012) to perform 3-D differencing on 40 × 40 m 2 windows from the 2005 and 2007 topography using a 20 m sliding window ( Figure S1A). We add a 5 m buffer to the 2018 data sets so that the rigidly deformed (e.g., translated) 2005 and 2007 data sets align with 2018 data. The displacement rate fields are shown in Figures 2 and 4 and are included in the Supporting Information S1.
There are three primary error sources in topographic differencing. (1) Point cloud errors from the laser measurement, inertial measurement unit (IMU), and GNSS position errors (Glennie, 2007;Goulden & Hopkinson, 2010;Toth et al., 2007). The 2005 B5 and2007 EarthScope data sets were acquired with older technology and have larger errors compared to the 2018 data. Flight line overlap alignment errors create centimeter-to decimeter-scale linear offsets in the NE and vertical displacement rate fields orthogonal to the sensor flight direction, while the flight-parallel NW results have lower error. GNSS errors in the B4 data set from atmospheric delays are spatially correlated and most pronounced over high topographic relief (Shan et al., 2007). Boresight error can be corrected if full mission instrument data are available (e.g., Glennie et al., 2014), but was not performed here. (2) In low relief landscapes, there may not be a unique result, because the point clouds can be well-aligned with multiple rigid deformations. (3) Contrasting point classification, data density, and anthropogenic change may not be removed with the 0.5 m vertical difference criteria filter.
We quantify error due to sources (2) and (3) using the displacement correlation error, following Scott et al. (2018). This approach assumes that errors are spatially uncorrelated and captures error associated with low relief landscapes and local issues, but not longer wavelength issues such as flight line alignment error. We calculate the displacement rate correlation error from the eight adjacent displacements that lie within a 30 m radius ( Figure S1B) and remove a planar ramp from each displacement component. The remaining standard deviation is the displacement correlation error. Over the 11 year-timespan of the EarthScope-FEMA data sets, the median horizontal and vertical correlation errors are 5.1 and 1.1 mm/year, respectively ( Figure S2). Error varies along strike ( Figures S3-S4), with larger errors in regions with lower topographic relief and some correlation to lower creep rates.

Creep Rate
To calculate the creep rate, we manually map the active trace of the CSAF from sharp discontinuities in the northwestern displacement rate field. We calculate creep rate from the displacement discontinuity method (Scott et al., 2018). Every 100 m along strike, we select 3-D displacements over 400 m along-strike and across the full 1,200-2,400 m data set width. We mask rates with correlation errors exceeding 50 mm/year, Figure 3. Right-lateral creep rate from topographic differencing in this study ("Topo" or black circles) and other geodetic data sets with a range of apertures. (a) Creep rate measurements from cultural features (Brown & Wallace, 1968;Burford & Harsh, 1980), creepmeter (CM; Burford, 1988), alignment arrays (AA; Burford & Harsh, 1980;Titus et al., 2006), short-range electronic distance measurements (EDM: Lisowski & Prescott, 1981), InSAR (Tong et al., 2013), and geodolite measurements (Lisowski & Prescott, 1981). Numbers in the legend indicate the timing of the first and last measurement of that data set; Most rates do not span the full length of time. (b) Differencing rates (black circles) compared to shallow right-lateral slip from InSAR and GPS slip inversions from Jolivet et al. (2015) and Khoshmanesh and Shirzaei (2018a). (c) Topographic differencing versus other rates. Explanation values are the other geodetic data sets minus the topographic differencing right-lateral offset rates (mm/year). A negative value indicates greater topographic differencing rates. horizontal displacement rates exceeding 100 mm/year, or that cover yet unmasked landslides (manually mapped by us). We calculate the inverse-error-weighted mean of the 3-D displacement rate components ([dNE, dNW, dZ]) along both sides of the fault, where w i is the normalized displacement rate correlation error (by component, ne, nw, z) from M rate measurements. The displacement rate error (Δd) is, The fault normal, right-lateral, and vertical creep rates (CR) are calculated by a vector subtraction of the mean rates (Equation 1) on southwestern (dsw) and northeastern (dne) sides of fault, We calculate creep rate error (ΔCR) from the displacement rate error (Equation 2) on both sides of the fault (Δdsw and Δdne), We assume that the fault strikes 315°, though the strike varies locally by ±15°( Figure S5). When calculating slip rates, we do not account for local strike variation because the overall plate boundary motion parallels the overall fault trend, and large flight line errors in the northeastern rate field would contaminate the estimated creep rate. The right-lateral creep rate error is 0.4 mm/year. This is lower than the correlation displacement error, because we assume that errors are spatially uncorrelated and use~450 rates on either side of the fault with an inverse error weighting.

Strain Field
We calculate the 2-D horizontal strain field to examine the distribution of deformation. We present results at Mustang Ridge in Figure 4 and results elsewhere in the Supporting Information S1. We calculate the 2-D horizontal displacement gradient rate tensor ( _ e ij ) following an approach similar to Allmendinger et al. (2012) and Scott et al. (2018). At each displacement rate, we solve for _ e ij from the displacement rates within a 45 m radius. We apply no distance or error-based weighting and remove rates with error exceeding 10 mm/year.

Creep Rate
The ICP-derived northwest rate field shows that the CSAF has accommodated right-lateral offset between Parkfield and San Juan Baptista over the 11-13 year time span, as shown in Figure 2. Along the length of the CSAF, the shape of the creep profile is asymmetric with higher creep rates in the southeast and lower rates in the northwest. At Parkfield, the fault creeps at a rate of 8.0 ± 1.6 mm/year, showing a sharp increase in creep rate to the northwest (all locations are in Figure 1). The maximum rate of 41.6 ± 4.9 mm/year occurs 5 km southeast of Slack Canyon. At Mee Ranch, the average rate is 33.3 ± 1.9 mm/year. The 20.4 ± 1.1 mm/ year rate at Dry Lake Valley is relatively low. Fault creep reaches a local maximum of 29.6 ± 3.2 mm/year at 90-105 km northwest of Parkfield near Pinnacles National Park. The CSAF creep rate decreases when the CSAF and CF overlap. The rates decrease on both faults from 105-125 km northwest of Parkfield. The CSAF rate increases at 125 km as the fault trace bends westward. We do not interpret the vertical and northeastern rates as the relatively low rates are indistinguishable from noise and are influenced by flight line error.

Comparison With Other Geodetic Creep Rates
We compare the 2005/2007-2018 lidar topographic differencing creep rates to other data sets and fault slip inversions (Figure 3). Aperture is likely the primary reason for rate differences, although the varying epochs capture slip rate transients such as those that occurred in the study area following the 2004 Parkfield Earthquake (Titus et al., 2006) and elsewhere due to pore fluid pressure changes (Khoshmanesh & Shirzaei, 2018b), earthquakes, or other causes (McFarland et al., 2019). All topographic differencing rates are larger than the 30 m cross-fault aperture creepmeter rates, with an average difference of 8.8 mm/year. The differencing-based rates are an average of 3.2 mm/year faster than alignment array rates (150 m median aperture). The electric distance measurements (EDM) in multiyear surveys from 1974-1980 with apertures of 2-3 km are an average of 1.7 mm/year lower than the differencing-based rate. At Mustang Ridge, the topographic differencing rates are 4-7 mm/year higher than the alignment array and EDM rates, likely because the alignment array does not span the full~1 km wide fault zone.
The L-band InSAR creep rates measured by Tong et al. (2013) from 2006-2010 over a 20 km cross-fault aperture are an average of 3.4 mm/year lower than the topographic differencing rates. Tong et al. (2013) used two independent line-of-sight look directions, assumed no vertical surface deformation, and estimated a 1.5 mm/ year rate error. Error in the InSAR and lidar topographic differencing results may explain the difference, but not the bias. Both data sets show a maximum rate near Slack Canyon and an asymmetric rate decay. Geodolite rates with line lengths often in excess of 15 km are an average of 6.2 mm/year faster than the topographic differencing rates and support the inference of increasing inferred creep rate with instrument aperture. However, particularly toward the ends of the creeping section, where the geodolite rates often exceed the other data set, the discrepancy may reflect elastic strain accumulation as the fault locking increases.
We compare topographic differencing rates to the shallowest fault patches in two slip inversions constrained with GNSS and InSAR data (Figure 3b). Jolivet et al.'s (2015) shallowest fault patches extend from 0-4 km depth. We use rates from Khoshmanesh and Shirzaei (2018a) fault patches at 0-2.5 km depth from 2005-2010, after the 2004 Parkfield earthquake. Neither study predicts a maximum slip that coincides with the topographic differencing peak creep rate nor the rapid rate decay approaching Parkfield.

Fault Creep in Area of Complex Faulting: Mustang Ridge
Topographic differencing reveals plate boundary deformation in structurally complex areas. Fault zone complexity is sometimes indicated by topographic features of apparent tectonic origin. For example,

Geophysical Research Letters
at Mustang Ridge, the fault crosses relatively high topography, as shown in Figure 4. Here, the fault steps rightward~1 km through a zone of distributed faulting indicated by en echelon scarps and sag ponds (Figure 4). Rymer et al. (1984) interpreted this area as a zone of distributed faulting responsible for the narrow-aperture creep rate being lower than those measured nearby along the CSAF. Using airborne lidar, DeLong et al. (2010) characterized the en echelon faulting network and the CSAF rotation along subsidiary shallowly dipping faults and proposed that the distributed faulting is responsible for the landform patterns. Topographic differencing supports these interpretations: The subsidiary faults that rotate short sections of the CSAF's main trace accommodate fault creep and transfer slip to a northeastern, parallel, fault trace (Figure 4b). Most or all fault creep steps~1 km to the right over~4 km along strike despite en echelon scarps and other evidence for faulting that span~10 km along strike (Figure 4a). The differencing results do not allow for precise partitioning of creep along each subsidiary fault but indicate that perhaps 3-5 of the subsidiary faults spaced less than 1 km apart accommodate most creep in the stepover region. The main trace of the CSAF is therefore~1 km northeast of the fault as mapped in the USGS/CGS Quaternary Faults and Folds database and by Rymer et al. (1984), along the northwestern portion of Mustang Ridge. This more active fault trace is less distinct in the lidar than the now-quiescent fault strand to the southwest previously thought to be the active CSAF.

Off-Fault Deformation
In earthquakes and along creeping faults, shallow deformation is often not localized along the main fault but is instead distributed in the shallow crust (e.g., Dolan & Haravitch, 2014). Burford and Harsh (1980), Lisowski and Prescott (1981), and Titus et al. (2006) interpret distributed off-fault deformation along the CSAF from higher deformation rates at larger apertures. With topographic differencing, we directly observe the localization of deformation along the fault trace ( Figure S6) as well as the structures that accommodate distributed deformation, such as the step-over at Mustang Ridge. Dolan and Haravitch (2014) attribute variations in the distribution of deformation to fault structural maturity and show that mature faults (>85 km of offset) accommodate 85%-95% of deformation along a narrow trace. The CSAF has over 315 km of offset (Mathews, 1976). Based on creepmeter and differencing rate comparisons, the CSAF accommodates 20%-50% of deformation off fault beyond a 30 m aperture. This exceeds the expectation from Dolan and Haravitch (2014) for the mature fault, possibly due to a slip rate dependence on the distribution of deformation. Off-fault deformation varies with geometric fault complexity: Deformation is distributed at Mustang Ridge but is localized at Dry Lake Valley where the fault is straight, consistent with observations from Scott et al. (2020).
There are several applications of this analysis: (1) Ongoing efforts to mitigate hazard from fault displacement at critical infrastructure (Baize et al., 2020;Petersen et al., 2011) require information about fault location, width, and slip from recent earthquakes. Applying these analyses to the CSAF would facilitate developing similar relations along the creeping fault, assessing differences between seismic and creeping deformation, and understanding the physical processes that guide off-fault deformation.
(2) These observations can inform research on how the mechanical properties of the shallow crust modulate deformation. Nevitt et al. (2020) used near-field geodetic data and finite element models to explore rheologic and other controls on shallow deformation for the 2014 M6 South Napa, California, Earthquake. Applying similar methods to Dry Lake Valley and Mustang Ridge could illustrate rheologic differences for the varying distribution of deformation. (3) By comparing inferred fault location from topographic differencing and geomorphology, we can understand which visible fault traces are currently active and where creep is being accommodated in absence of clear geomorphic indicators. Such analyses can demonstrate fault zone evolution and facilitate better approaches for understanding which faults are active based on the geomorphic record.

Along-Strike Fault Creep Asymmetry
The along-strike creep distribution is asymmetric with the maximum rate occurring 20 km northwest of Parkfield. In earthquakes, the slip distribution is also often asymmetric due to regional stress field gradients, varying damage or rock type along strike, geometric bends, and the partition of slip along adjacent faults (e.g., Burgmann et al., 1994). We consider how such features may control creep asymmetry along the CSAF.

Geophysical Research Letters
The CSAF separates Salinian Block crystalline basement to the southwest from basement rocks of the Franciscan Complex and Great Valley Sequence to the northeast (Irwin, 1990), and at shallower depths is bordered by a wide variety of Paleogene and Neogene marine and terrestrial sedimentary units that contrast across the fault. However, the regional geological fabric is slightly oblique to the CSAF. The rock types surrounding the CSAF may modulate the creep profile's shape by changing what geological materials are along the fault surface, particularly as the Great Valley Sequence loses width near Parkfield and Cholame.

Strain Accumulation and Seismic Potential
Given the 39 ± 0.2 mm/year relative Pacific-Sierra Nevada Great Valley plate boundary motion (Argus & Gordon, 2001), our results show variable and generally low moment accumulation potential along the CSAF, with at most several mm/year moment accumulation near Slack Canyon. Creep rates are the lowest from Bitterwater to Dry Lake Valley. This area hosted three M5.7-6.5 earthquakes from 1882-1885 (Toppozada et al., 2002) as shown in Figure 1 suggesting that accumulated strain in the absence of larger creep rates may be released in moderate-sized earthquakes or that creep rates vary in time. Creep rates decrease to the southeast at Parkfield where historical M6 events have released seismic moment (Lienkaemper & Prescott, 1989;Toke & Arrowsmith, 2006).
Fault creep rates indicate reduced seismic moment accumulation (e.g., Field et al., 2013). However, a major challenge to constraining the reduction due to creep along partially locked faults is uncertainty in how creep rates change with depth (Weldon et al., 2013). Existing slip models for the creeping SAF show differences in slip with depth and along strike (Jolivet et al., 2015;Khoshmanesh & Shirzaei, 2018a;Maurer & Johnson, 2014;Ryder & Bürgmann, 2008). For the M7 2016 Kumamoto, Japan earthquake, Scott et al. (2019) show that near-field topographic differencing, optical correlation, and InSAR displacements add resolution to shallow fault slip. CSAF topographic differencing results are likely to contribute toward better resolving shallow slip and constraining seismic moment accumulation in the shallow crust.

Conclusions
We quantify the tectonic creep rate along the CSAF by differencing lidar topography acquired between 2005-2018. We measure creep rates and strain with unprecedented spatial resolution. We show that creep is asymmetrically distributed along the length of the CSAF and that strain localization varies along strike. Our analyses highlight new approaches for measuring fault slip using topographic differencing when the slip magnitude is just above the data set error, extending airborne lidar differencing to creeping faults and earthquakes with decimeter-scale offset. Our lidar topographic differencing results add new constraints to deformation patterns along the CSAF and support future work to estimate seismic moment accumulation in the shallow crust, decipher which geomorphic fault traces are actively creeping, and probe the mechanical properties of the shallow crust.