Research ArticleSeismology

Buried shallow fault slip from the South Napa earthquake revealed by near-field geodesy

See allHide authors and affiliations

Science Advances  28 Jul 2017:
Vol. 3, no. 7, e1700525
DOI: 10.1126/sciadv.1700525


Earthquake-related fault slip in the upper hundreds of meters of Earth’s surface has remained largely unstudied because of challenges measuring deformation in the near field of a fault rupture. We analyze centimeter-scale accuracy mobile laser scanning (MLS) data of deformed vine rows within ±300 m of the principal surface expression of the M (magnitude) 6.0 2014 South Napa earthquake. Rather than assuming surface displacement equivalence to fault slip, we invert the near-field data with a model that allows for, but does not require, the fault to be buried below the surface. The inversion maps the position on a preexisting fault plane of a slip front that terminates ~3 to 25 m below the surface coseismically and within a few hours postseismically. The lack of surface-breaching fault slip is verified by two trenches. We estimate near-surface slip ranging from ~0.5 to 1.25 m. Surface displacement can underestimate fault slip by as much as 30%. This implies that similar biases could be present in short-term geologic slip rates used in seismic hazard analyses. Along strike and downdip, we find deficits in slip: The along-strike deficit is erased after ~1 month by afterslip. We find no evidence of off-fault deformation and conclude that the downdip shallow slip deficit for this event is likely an artifact. As near-field geodetic data rapidly proliferate and will become commonplace, we suggest that analyses of near-surface fault rupture should also use more sophisticated mechanical models and subsurface geomechanical tests.


The nature of earthquake-related fault slip in the shallowest portion of Earth’s crust is essentially unknown. Although myriad surface mapping [for example, the study by Lawson (1)] and paleoseismological [for example, the study by Sieh (2)] studies document a rich spatiotemporal history of surface- or near-surface breaking earthquakes, they typically concern themselves only with measurements made at the surface or, in trench excavations, of a few-meters depth. Inferring shallow fault slip any deeper than that has been limited by sparse surface observations in the near field (within ~1 km) of a fault’s trace. For instance, seismological and Global Positioning System (GPS) geodetic stations are typically spaced too widely (tens to hundreds of kilometers), whereas satellite-based interferometric synthetic aperture radar (InSAR) tends to decorrelate in the near field because of the effects of large ground accelerations. Because the resolving power of slip with depth is proportional to observational distance across the fault [for example, the book by Stein and Wysession (3)], seismological and geodetic fault slip analyses typically do not infer fault slip above 1-km depth.

Without better quantification of shallow fault slip, we are likely limiting not only our understanding of fundamental fault-slip physics but also seismic hazard associated with near-surface rupture. If, for instance, near-surface vertical gradients in fault slip are not negligible and it is assumed that surface observations of fault-parallel displacement are equivalent to slip across a fault at shallow depths, then slip rates used in probabilistic seismic hazard assessments (4) could be underestimated. In addition, source models of recent continental strike-slip earthquakes document an apparent deficit in shallow fault slip compared to slip at depths of a few kilometers (5). It has been unclear, however, whether this reflects an artifact of inadequate observation and modeling of shallow fault slip (6, 7) or whether it represents the fundamental mechanical behavior of Earth’s shallow crust. Possible explanations include a tendency toward velocity strengthening friction of the shallow crust that would inhibit earthquake rupture (8, 9) and the possibility of distributed “off-fault” deformation (10, 11) associated with “compliant” zones of damage around narrow fault zones (1217) or anelastic yielding of the country rock immediately surrounding a seismogenic fault (18, 19).

Recently, a new generation of imaging techniques is filling in the observational gap and providing unprecedented, spatially dense, near-field observations. Space-based optical image correlation (2025), airborne InSAR (26), airborne laser scanning (11, 2729), terrestrial laser scanning (30), and mobile laser scanning (MLS) (31) can document near-field displacements at the decimeter to subcentimeter level with observational density ranging from 1 to 10,000 points/m2 and spatial coverage from meters to kilometers. To date, the earthquake-related studies using these types of data have focused on arithmetic characterizations such as the difference of “far” and “near” field surface displacements (21, 22, 24, 32, 33). The assumption is that these metrics quantify off-fault deformation (OFD) that, in turn, scales with any deficit in shallow slip. To the best of our knowledge, no studies have used near-field data in concert with geodetic inverse methods [for example, the work by Xu et al. (7)] to develop explicit slip models in the shallowest crust (<100 m) that allow for variation in both slip magnitude and depth of faulting.


The South Napa earthquake

The 24 August 2014 earthquake initiated at ~10-km depth and ruptured updip and north along a near-vertical, approximately north-northwest-striking plane (3441) that reactivated a previously mapped strand of the West Napa fault system (Fig. 1A) (42). Dextral surface rupture, expressed mostly as discontinuous en echelon mixed-mode fractures, offset curbs, and buckled pavement, was distributed over a zone ~15 km long and ~5 km wide (34, 43, 44). Slip was concentrated dominantly on the westernmost strand, and our study addresses only this principal rupture (Fig. 1A). Hand measurements made by us and others in the first few hours after the event documented 5 to 50 cm of coseismic displacements in the northern half of the rupture (43, 44). Similar displacement measurements in the southern portion of the rupture were all less than 5 cm. Over the southern half of the rupture (more than ~8 km), vigorous postseismic slip had begun by at most ~3 hours after the main shock (43). By 24 hours after the event, we measured ~20 cm of fault-parallel surface displacement at a location ~ 6 km north of the epicenter. Decimeter-scale postseismic surface displacements accrued in the subsequent months with a spatial concentration in the southern part of the rupture, a pattern complementary to coseismic slip that was focused in the northern part of the rupture (38, 41, 45, 46). The northern, coseismic portion juxtaposes primarily Cretaceous Great Valley sedimentary rocks with Pleistocene alluvium, whereas the southern postseismic portion juxtaposes Plio-Pleistocene and Pleistocene alluvium in a southward-deepening basin reaching depths of ~3 km (47).

Fig. 1 Study area of the 2014 South Napa earthquake and vine row images.

(A) Airborne light detection and ranging (LIDAR) base map of the study area. Buhman Avenue, the northern trench (T2); South Avenue, the southern trench (T1). Inset: Regional location map with major fault strands of the San Andreas Fault system. (B) Photograph of offset vine rows. Arrows indicate right-lateral sense of motion. (C) Map view examples of MLS point clouds from blocks 1 and 7. (D) Example vine row point cloud and elastic dislocation models. Inset: Same model curves as main figure, except with larger across-strike distance. (E) Cartoon illustrating the inversion schematic. S, fault slip; zl, lower dislocation limit; zu, upper dislocation limit; uy(x), fault-parallel surface displacement.

Immediately following the event, we made vehicle-mounted MLS scans over much of the rupture in two epochs (1 to 2 September and 28 to 30 September) that yielded more than 1010 laser reflection points three-dimensionally (3D) referenced to better than a few centimeters (Fig. 1 and Materials and Methods) (31). The rupture crossed multiple Napa Valley vineyards nearly perpendicular to vine rows that had been planted in straight lines (Fig. 1, B to D); thus, any curvature of the rows can be confidently attributed to the South Napa earthquake. A typical MLS “point cloud” from a vine row comprises over 50,000 laser returns with centimeter spacing and images two relatively straight limbs offset parallel to the rupture trace by tens of centimeters and connected by a sigmoidal portion ~10 to 30 m wide. Within this curved zone, there is almost always a narrow zone of surface disruption comprising en echelon sheared extensional fractures and linear “mole tracks” (Fig. 2). Nowhere, including in two excavated trenches (Fig. 2, D and E, and Materials and Methods), do we observe a through-going discrete rupture plane.

Fig. 2 Surface and near-surface expression of faulting.

(A) Base map, optical image collected on 9 September 2014. Red lines, en echelon fractures mapped by hand. Blue line, location of trench in (D). (B and C) Field photos of surface expression of faulting. Yellow lines, en echelon fractures. (D) Trench log from the South Avenue trench. Orange lines, fractures. Black square, location of sample for mechanical tests. (E) Interpreted photomosaic of trench at Buhman Avenue. Brown lines, fractures. See Fig. 1A for location of trenches.

A shallow slip model

Our objective, as in other geodetic or seismological analyses, was to formally invert the near-field data for a model of shallow fault slip. We recognize that there are a multitude of potential mechanically distinct models describing this process (12, 18). Computational limitations require us to find acceptable combinations of model physics and complexity, data set size, and the thoroughness with which we wish to examine model parameter space. For instance, the South Napa near-field fault-parallel displacement data exhibit maximum strains of 1 to 5% (50 cm over 10 to 50 m) (Fig. 1, C and D) that are somewhat higher than the ~0.5% strains for which elasticity is commonly assumed. However, for inverting the exceptionally dense near-field data, an elastic forward model for which there are analytical solutions (48) is computationally tractable, whereas a more mechanically realistic elastoplastic model solved via the finite element method, for example, requires too much computational expense. To justify our choice of an elastic model for our inversion and gain insight into the resolving power of the near-field data, we compare forward models of vertical, strike-slip dislocations (hereafter referred to as “faults”) that follow elastic, laterally varying elastic, and elastoplastic constitutive laws (Materials and Methods).

Generally, surface displacements are sensitive to the depth of the fault’s upper edge: For the same amount of slip, fault edges reaching closer to the surface produce greater horizontal fault-parallel displacement gradients (duy/dy) (Fig. 1D). Moreover, each of the models approaches the same displacement values ~150 m from the fault. Without data from the near field, the different models would be indistinguishable. For the same amount of slip, an elastoplastic model with material properties most appropriate for the South Napa rupture produces more near-field fault-parallel displacement than the elastic model (Materials and Methods and fig. S1). The peak of the displacement profile is closer to the surface projection of the fault, resulting in a narrower fault zone that reflects the zone of high plastic strain in the region immediately above the dislocation edge. Alternatively, a compliant zone model with laterally varying shear modulus of 50% yields surface displacements in the nearest field very similar to the elastic half-space model: There is essentially no difference out to ~10 m from the fault trace, and out to ~50 m the difference is a ~10% decrease (fig. S1) (12).

For a given vine row, we invert the surface displacements uy for fault slip S using the analytical solution for a 2D screw dislocation in a homogeneous elastic half-space between upper (zu) to lower (zl) dislocation depths (Materials and Methods) (3): Embedded ImageWe fix zl conservatively at 5000 m to prevent edge effects and perform a grid search on zu solving for S using a weighted least-squares solution (49):m = [G′ * C−1 * G]−1 * G′ * C−1 * d, where m is a vector of the estimated parameters, d is a vector of fault-parallel surface displacements, C is the matrix of variances, and G is the matrix of Green’s functions relating slip at depth to surface displacement. We assess potential bias in the choice of an elastic model for G by running the inversion on the elastoplastic and compliant zone forward model examples. As to be expected from its narrower zone of surface displacements, for the elastoplastic model, the elastic inversion slightly overpredicts fault slip (by 5%) but underpredicts the depth to the upper fault edge by as much as 60% (fig. S2). In contrast, the elastic model underestimates the compliant zone model slip by ~15% and its depth by ~5% (fig. S2). The depth bias renders our depth of burial estimates conservatively shallow.

We create a shallow slip model for each epoch of data collection by inverting, independently of one another, 691 vine rows distributed along the trace of the fault (Fig. 3). To the best of our knowledge, such a model is unprecedented. For each row in a vineyard block, we plot the posterior probability distribution function (PDF) for S (Fig. 3, middle panels). We also estimate fault-parallel surface offset uy (Fig. 3, black circles, middle panels; fig. S3; and Materials and Methods), similar to what would be measured by a geologist with a tape measure. Generally, the peak values of the PDF for S range from ~0.5 to 1.25 m, whereas uy varies from ~0.25 to 0.5 m. Because they are estimated independently, the difference between uy and S is an outcome rather than a requirement of our analysis. The lower panel (Fig. 3) displays, for each row, the posterior probability distribution for zu. Peak values of the PDF for the upper depth of the fault, zu, vary from ~3 to 25 m below the surface and tend to be shallower in the dominantly postseismic portion of the fault (Fig. 3, B to H) than for the coseismic portion (Fig. 3A). Where zu is deepest, such as in blocks 1 to 5 (Fig. 3, A to E), uy values lie toward the lower tails of the S PDF. Where zu is shallower, such as in blocks 6 and 7 (Fig. 3, F and G), uy values are closer to the peaks of the S PDF. Generally, the vine row curvature across the fault zone is reproduced exceptionally well by the elastic model (Fig. 3, top panels). We stress, however, that the good fit between data and model does not require the elastic model to be the unique representation of shallow faulting. We choose the elastic model for its computational efficiency and because, with the potential biases described above, it allows us to place constraints on slip and depth-of-faulting estimates.

Fig. 3 Inversion results.

(A to H) Results from the indicated vineyard block. Results are for the second-measurement epoch (28 to 30 September) for each vineyard block (see Fig. 1A for block locations and fig. S4 for similar plots for first-measurement epoch). For each block: Top, point cloud (gray points) of example row with optimal model (red line). Middle, probability distribution for slip S. The colors represent probability p. Black circles, estimate of fault-parallel surface displacement uy. Bottom, probability distribution for upper edge of dislocation zu. White dotted line, location of the example row.

Because the earliest measurements were collected 7 days after the event, our shallow slip model for each epoch reflects cumulative coseismic and postseismic slip. By considering both the field measurements taken by us and others in the hours immediately following the main shock, as well as the subsequently established alinement arrays, we determine that block 1 reflects dominantly coseismic slip (Fig. 3A); this is also where the maximum coseismic offset was measured (34, 43). The remaining seven blocks reflect dominantly postseismic slip (Fig. 3, B to H, and fig. S4).

A buried slip front

Our analysis demonstrates that despite clear surface disruption in the form of mole tracks and discontinuous en echelon fractures, significant fault slip across a dislocation discontinuity did not reach the surface either coseismically or postseismically. This is confirmed in the two trenches: one in each of the dominantly coseismic and postseismic portions of the rupture (Fig. 2, D and E, and Materials and Methods). In neither of the trenches is a discrete shear dislocation plane with decimeters of slip observed despite the obvious surface disruption and decimeter-scale surface displacements immediately adjacent to the trenches. In the southern trench, three discrete shears attributed to the 2014 rupture could be traced to the surface, but no significant slip was resolved on any of them. In the northern trench, no shear dislocations could be traced to the surface or attributed to the 2014 event. Thus, significant fault slip must have terminated deeper than the ~3-m floor of each trench. This is consistent with our estimates of zu ≥ 3 m for each of the blocks nearest the trenches (Fig. 3, A and G).

The fault zone that ruptured in the South Napa earthquake had been previously recognized as a strand of the West Napa Fault zone (42). Accordingly, we interpret the estimates of zu as a map of the position of a slip front—either coseismic or postseismic—that propagated on an extant fault plane. In both the coseismic and postseismic portions of the fault, the slip front reached close to the surface very rapidly. For the region of coseismic slip, this was presumably contemporaneous with the earthquake shaking—geological observations of offset vine rows in block 1 were made within hours of the main shock and are equivalent to our estimates of uy made after 1 week (Fig. 2A) (43). For the postseismic portion, within a few hours of the main shock, there were observations of road displacements associated with vigorous afterslip (43), and so, the front must have propagated to a location close to the surface within, at most, a few hours of the main shock. If, below the postseismic portion, coseismic slip reached to within ~2 km of the surface (3841) and postseismic slip was driven by the corresponding stress gradient imposed at those depths (50), then the postseismic slip front propagated upward at a velocity of order 1 km/hour.

The slip front reached closer to the surface (~3 m) in the southern part of the postseismic portion (Fig. 3, F and G) than it did in the coseismic portion, where it terminates ~5 to 7 m below the surface (Fig. 3A). However, there is not a monotonic decrease in slip-front depth from north to south. Toward the middle of the rupture, blocks 2 and 5 exhibit the deepest zu estimates of ~25 m (Fig. 3, B and E). Within each block, the slip front undulates along strike, with depth gradients (dzu/dy) varying by meters over distances of meters. The undulating nature of the depth gradients is the same regardless of the surficial lithologic units through which the fault cuts (fig. S5). The vine rows from blocks 2 to 4 and the northern portion of block 5 are in a portion of the fault that puts Cretaceous sedimentary units against Quaternary alluvium, whereas in all of the other vineyard blocks, the fault contact is entirely in Quaternary alluvium. We know of no other studies with which we can compare the along-strike depth gradients. A full 3D approach (out of the scope of this paper because of excessive computational time requirements), however, is probably more appropriate with depth gradients as strong as meters over distances of meters; we leave further examination of slip-front depth variation for future studies. We infer that along-strike gradients in S are typically less than tens of centimeters over ~100 to 200 m. This is somewhat less than the results for recent surface-rupturing strike-slip earthquakes documenting along-strike gradients of meters over hundreds of meters (21, 22, 24, 29, 32).

Slip deficits and OFD

Our results demonstrate the detailed complementary spatial relationship between coseismic and postseismic slip processes for the South Napa event (Fig. 4) (38, 41, 45). The hand measurements made at the day of the event (34, 44) document a peak in coseismic surface offset in the northern portion of the rupture and a deficit in the south (Fig. 4A). In the south, our analysis portrays progressive postseismic slip filling in the shallow coseismic slip deficit (Fig. 4). By the first week after the event, in the northern portion of the postseismic zone, median slip for each block exceeded ~25 cm and peak values approached ~70 cm. By 1 month after the event, median postseismic slip had extended further to the south, increasing by at least 10 cm everywhere and reaching peak postseismic slip values (~80 cm) close to the maximum observed coseismic surface offset values. We calculate the cumulative slip for each MLS measurement epoch and find that by 1 week, the cumulative slip of the postseismically slipping portion had reached at least 47% of that of the coseismic portion. By 1 month, near-surface postseismic slip accounted for at least 87% of near-surface coseismic cumulative slip. This is in agreement with GPS (38) and alinement array time series (45) showing that by 30 days after the event, shallow postseismic slip had reached close to its maximum amplitude. We note, however, that because of the effect of uy, described above, being less than S for a buried fault, the alinement array measurements tend to underestimate significantly the amount of shallow slip (Fig. 4A).

Fig. 4 Comparison of shallow (this study) and deeper (3941) slip estimates.

(A) Colored columns, stacked PDFs of S from each vineyard block for the second-measurement epoch (Fig. 3, middle panels). The width of the column represents the along-strike length-scale of the vineyard block. Black circles, hand measurements of fault-parallel surface displacement made immediately following the earthquake. Gray circles, alinement array measurements 60 days after the event (45). (B) Colored columns, stacked PDFs of zu from each vineyard block for the second-measurement epoch (Fig. 3, bottom panels). (C) Coseismic, kinematic slip model from Wei et al. (41). (D) Normalized cumulative slip from three coseismic slip models. W, the study by Wei et al. (41); J, the study by Ji et al. (36); M, the study by Melgar et al. (40). Colored rectangles in the shallowest portions are cumulative slip S from all the vineyard blocks in this study, normalized according to the maximum value of each of the kinematic models.

This rapid, near-surface postseismic slip equilibration erased a slip deficit that was present along the strike of the rupture. Downdip, however, all of the coseismic kinematic models (that do not include postseismic measurements) (3941) display an apparent shallow slip deficit of 40 to 80% relative to the maximum slip values at depth (Fig. 4, C and D). The shallowest slip estimates from these models (at depths of ~1 km) are roughly in agreement with our near-surface slip estimates even after ~1 month of postseismic slip (Fig. 4D); therefore, rapid postseismic slip was not sufficient to erase the apparent shallow slip deficit. Formally including near-field geodetic data in slip inversions, however, significantly reduced apparent shallow slip deficits from some recent continental strike-slip earthquakes (7). Thus, either (i) the South Napa apparent shallow slip deficit is an artifact of the inversion process (7) or, if the shallow slip deficit from the kinematic models is real and our results show that postseismic slip is likely not enough to erase it, then (ii) significant “off-fault” deformation associated with the fault slip must have occurred (18).

The term “off-fault deformation” has been used to describe either plastic deformation in a zone surrounding the fault in dynamic rupture models (8, 18) or a zone of distributed deformation surrounding the fault including some combination of plastic deformation and brittle deformation from slip on subparallel faults or block rotations (10, 21, 24, 33, 51). Recent attempts to quantify OFD have used some combination of the difference of far-field (measured outside the zone of high shear strain near the fault trace) and near-field (measured at the fault trace) fault-parallel offset (21, 33). However, following these measures, the reduction in surface displacement due to a buried fault would yield apparent OFD without needing to appeal to plastic or secondary faulting (Materials and Methods). Thus, only in the end-member case when all slip reaches the surface can surface measurements alone quantify OFD. In the more general case when some degree of buried faulting and gradient in slip with depth are expected, additional information is needed to quantify OFD.

For the South Napa rupture then, without undertaking comprehensive mechanical analyses of the near-surface materials via widespread coring and laboratory testing, we cannot assess the degree of plasticity and the potential for OFD in the shallow section. We can, however, make a qualitative argument to suggest that there is negligible OFD. We know that there is relatively high clay content in the Napa Valley sediments (5254) and that this is consistent with relatively high cohesion values and relatively narrow distributions of plastic strain in elastoplastic formalisms of faulting (Materials and Methods and figs. S1 and S2) (18). In addition, we performed mechanical tests on three samples from one of the trenches (Materials and Methods) and found 50% clay content with a coefficient of friction of ~0.40. These values along with the excellent agreement between our elastic model and the data are consistent with relatively low amounts of plastic deformation (18). These factors, then, in conjunction with the lack of observations of significant subparallel secondary faulting or rotated blocks, lead us to suggest that OFD associated with the South Napa rupture is most likely insignificant.


The ability of the MLS platform to acquire extremely-high-resolution near-field geodetic data over a large portion of the South Napa rupture permits us to estimate a model of coseismic and postseismic shallow slip of exceptional resolving power and detail. Combined with evidence from trenches, our analysis leads us to conclude that despite its classification as a “surface-rupturing” earthquake and despite surface disruption and measurable surface displacements of offset features, slip on the principal strand did not breach the surface. Rather, slip terminated at varying depths along strike from ~3 to as much as ~25 m. For larger events, M 7 and above, that undoubtedly produce surface breaks (55), our study implies the possibility that surface fault slip measurements should not necessarily be considered representative of fault slip at depths of even a few meters below the surface. This should not necessarily be surprising; however, in the most thorough study known to us of the expression of faulting in trenches, 67 to 70% of 72 strike-slip fault strands associated with individual faulting events could not be traced to the surface, but rather, they died out centimeters to meters from the ground surface at the time of faulting (56).

An important hazards implication of our study is that surface-faulting analyses equating surface measures of displacement, uy, with actual fault slip, S, can be susceptible to underestimation error (fig. S6 and Materials and Methods). This potential error can introduce an underestimation bias not only to previous characterizations of coseismic surface rupture (57) but also to empirical slip scaling relations (58) and fault slip rate calculations, especially shorter-term ones, used in seismic hazard assessments [for example, the report by Field et al. (4)]. We demonstrate the magnitude of the potential error by nondimensionalizing Eq. 1 by zu and plotting the ratio of uy and S for a range of zu from 0 to 100 m (Fig. 5). If measurements are made within a distance of x/zu < 5, then the error of any individual measurement will exceed ~15%. Without high-resolution near-field geodetic data, recognizing the subtle curvature of the surface displacement field could be difficult in the field.

Fig. 5 Plot of the ratio of fault-parallel surface offset uy and slip S.

Derived from the 1D screw dislocation solution (3) and plotted versus distance x normalized by upper edge depth zu.

Why would a slip front propagate through kilometers of Earth’s crust only to arrest meters from the surface? One possibility is that shear waves, which refract into nearly vertical ray paths in the near surface, cause plastic yielding before the dynamically propagating slip front arrives at that depth (59). Alternatively, as the front approaches the surface, it may preferentially promote mode III fracturing. The result would be that the primary rupture would remain buried and that the surface shear displacement would be accommodated through a combination of frictional failure and secondary fractures, such as observed in the South Napa rupture (Fig. 2) and in a myriad of other surface disturbing events (57). Furthermore, sediments and poorly consolidated near-surface materials, when combined with decreasing normal stress, lead to the likelihood of more ductile material in the shallow fault. Thus, as the rupture front approaches Earth’s surface, there is a tendency to shift from a discrete fault plane to distributed deformation (44). For the postseismic slipping portion of the rupture, this mechanism could also be valid, albeit with a slip front propagating more slowly than in the coseismic case. A deeper understanding of the physics controlling slip-front arrest would have direct implications for earthquake engineering efforts. For instance, ground motion prediction equations for the western United States are explicitly dependent on the depth to the top of rupture (60). Earthquakes that do not rupture the surface systematically produce greater short- and intermediate-period ground motions than their surface-breaking counterparts (61, 62). However, the data constraining these depths for historic earthquakes are especially sparse, and, as described above, the physical mechanisms relating the arrest of a rupture to ground motions are only beginning to be explored.

Although our study presents the most detailed estimate of shallow slip for a near-surface rupturing earthquake that we know of, we are still not able to say much quantitatively about an apparent shallow slip deficit. We find that by ~1 month after the earthquake, shallow postseismic “afterslip” had equilibrated along strike of the portion of the fault that ruptured coseismically. This is similar to the result using GPS and InSAR data (38) and spatially sparse alinement array data (45). Although comparison with the coseismic kinematic models indicates an apparent vertical shallow slip deficit, our analysis shows that postseismic slip is not enough to make up the deficit and that there does not appear to be the significant OFD needed to create a real shallow slip deficit. Our reasoning, then, leaves us with the explanation that some combination of elastic Green’s functions (18) and lack of near-field data in the inversions (7) probably exacerbates the apparent shallow slip deficit in the kinematic models. More thorough exploration via joint inversion of the near-field and farther-field data (GPS and InSAR), which is out of the scope of this paper, is the subject of ongoing work.

With the proliferation of near-field geodetic data, we are entering a new phase of exploration of earthquake-related fault slip and physics in the shallowest portion of Earth’s crust, paradoxically a place we know little about. We anticipate that each new surface-disturbing earthquake will be accompanied by different types of near-field geodetic data. A high priority will be to merge these data with intermediate-field and far-field measurements (for example, strong ground motion, InSAR, GPS) in joint inversions that permit more complete analyses of the vertical distribution of slip from kilometers to meters depth. Our work on the South Napa earthquake, however, also illustrates that the new geodetic data and analyses should be accompanied by supporting mechanical data. Although certain parameters such as depth of faulting and slip magnitude can be constrained with dense near-field data and elastic models, a deeper understanding of fundamental processes such as slip-front arrest and the distribution of on-fault deformation and OFD requires the application of more sophisticated mechanical and dynamic models. In practice, this should then allow for significantly improved ground motion prediction equations for the near-field (60, 62) and postseismic slip forecasts in areas traversed by critical lifelines (45, 63).


This section contains detailed descriptions of the methodology used to support the main text of the paper. The following sections are in order of appearance in the main text: Fracture mapping, Mobile laser scanning, Trench excavation, Comparison with elastoplastic model, Estimation of fault-parallel surface offset, Shallow slip model, OFD calculations, and Laboratory testing.

Fracture mapping

Ground breaks were mapped on a 0.3-m resolution 2011 NAIP (National Agriculture Imagery Program) imagery orthophoto base map at 1:300 scale (1 inch = 25 feet or 1 cm = 3 m), which is the practical limit of the imagery resolution in this area (Fig. 2A). Rupture mapping was digitized using ArcGIS. Absolute accuracy location of the features is estimated to be 1 m in the built environment (school parking lot or backyards) and 2 to 3 m in vineyards, although relative accuracy of the mapped features is much higher.

Mobile laser scanning

To provide the highest-quality geodetic point cloud, we time-tagged all MLS and navigation data and recorded them to hard drives for post-mission processing. The core of an accurate MLS point cloud is the combined global navigation satellite system (GNSS)/inertial navigation system (INS) used to precisely tag laser observations and provide high-frequency estimates of laser scanning platform position and orientation. For the Napa missions, we used an iXblue LandINS with an internally integrated Trimble BD960 dual-frequency GPS/GLONASS (Globalnaya Navigazionnaya Sputnikovaya Sistema) receiver card. The LandINS recorded raw accelerations and angular rotation rates, along with dual-frequency carrier phase GPS observations. The LandINS also provided a PPS (pulse per second) signal aligned with GPS time to time tag the raw laser observations from the Riegl VZ-400 laser scanner. Post-mission, the INS and GNSS observations were processed using carrier phase processing using a tightly coupled Kalman filtering approach (64) in the software package Inertial Explorer (65), resulting in a 100-Hz estimate of platform position and attitude. The processed trajectory was then combined with the time-tagged laser observations and boresight calibration parameters to produce a georeferenced point cloud (66). The calibration parameters were iteratively refined using a least-squares planar matching optimization technique (67). Finally, the georeferenced point cloud is then segmented into ground and above-ground (that is, vine rows) features using the commercial software package TerraScan. TerraScan uses a ground filtering algorithm similar to that described in the study by Axelsson (68).

Trench excavation

South Avenue. A 20-m-long trench was excavated at the South Avenue site (Figs. 1 and 2) to examine the paleoseismic record and near-surface structure of the principal strand of the South Napa earthquake rupture. The site was selected based primarily on accessibility; the property owner was replanting the vineyard that occupies the property and was amenable to us conducting a paleoseismic investigation during November 2014. The site is located on an alluvial surface mapped as Pleistocene alluvium (Qoa) (69), and the hope was that a record of past surface-rupturing earthquakes would be recorded in these sediments. A ~2.5-m-deep trench was excavated, and the south wall of the trench was cleaned and gridded with a 1-m-wide × 0.5-m-high string grid. The exposure was logged at a scale of 1:200 on grid paper, and photos of the exposure were also taken (Fig. 2D).

The trench exposed a sequence of relatively flat-lying fluvial clays, silts, sands, and gravels faulted by a ~2-m-wide zone of well-developed shears (Fig. 2D). Only three discrete shears attributed to the 2014 rupture could be traced to the surface, rupturing through a ~30-cm-thick unit disrupted by agricultural plowing. Other faults, indicating previous episodes of movement, appeared to be truncated at the base of the plow zone or lower. Below the plow zone was a 30- to 60-cm-thick argillic horizon with a 30- to 40-cm down-to-the-west component of slip. This argillic horizon is attributed to soil development processes and, given the thickness and degree of development, indicates that the geomorphic surface the trench was placed across has likely been stable for tens of thousands of years. However, no datable material (carbon) was found in the trench to provide radiometric dating of the deposits.

Although the well-developed shear zone is indicative of repeated episodes of faulting, no evidence of discrete faulting events (for example, fissure fills, colluvial wedges) were identified in the exposure. It is likely that the lack of deposition at the site during the time it took to develop the argillic horizon prevented the preservation of individual earthquakes in the stratigraphy, making this site less than ideal for recording past earthquakes. A 20-cm-thick gravel deposit is located on the west side of the fault zone, about 30 cm above the bottom of the trench, and is truncated at the westernmost edge of the fault zone. This gravel is not present on the east side of the fault, and this relationship is interpreted as evidence of lateral slip that occurred before 2014, because it would take a significant amount of lateral slip, likely multiple meters, to laterally fault this gravel out of the plane of the trench exposure on the east side of the fault zone.

Buhman Avenue. The Buhman Avenue trench site is located across a short, preexisting ~2-m-high scarp formed on colluvial deposits. The rupture is weakly expressed in the shallow subsurface and difficult to identify; however, well-developed shear surfaces with local clay development and juxtaposition of units are better expressed at depth, which we interpret as the result of an unresolved number of previous earthquakes. Radiocarbon dates confirm that the faulted deposits are Early Holocene in age.

Deposits in the trench are poorly stratified, consisting primarily of massive clay with localized sand and gravel interbeds. The section appears to reflect colluvial deposition near the base of the local slope.

Comparison with elastoplastic and compliant zone models

As discussed, for instance, by Kaneko and Fialko (18), and reproduced by us with a finite element model, elastoplastic models of shallow fault-related deformation using a Drucker-Prager constitutive relation yield characteristic surface displacement profiles that reflect relative contribution of three fundamental parameters: (i) rock cohesion, c; (ii) internal angle of friction, ϕ; and (iii) fault friction, μ. Quasi-static 3D finite element models were constructed using Abaqus version 6.11. The vertical, planar model fault measures 400 m in length and 300 m in width, with the upper edge buried 5 m below the surface. The fault sits in a homogeneous model domain that measures 1200 m × 1200 m × 800 m. We use 198,149 second-order tetrahedral elements with biased node seeding to allow for finer meshing adjacent to the fault. Element size ranges from 1 to 100 m. Loading is introduced in two steps: (i) gravitational load to equilibrate lithostatic pressure gradient, and (ii) prescribed uniform slip of 1 m on the fault (uy = ±0.5 m on each fault surface). Assuming linear elastic constitutive properties, this model setup was benchmarked using the Okada (70) analytical solution for a rectangular dislocation, which showed that the numerical solution was accurate to within 1.3% of the analytical solution at ±50 m from the fault.

Elastoplastic models. For elastoplastic modeling, we implement the Drucker-Prager yield criterion with associated plastic flow. On the basis of mechanical testing (see below) and analysis of wells in the region (5254), we assume that the following mechanical properties are representative of the setting for the South Napa rupture: c = 50 kPa, ϕ = 25°, and μ = 0.4. These values were converted from Mohr-Coulomb to Drucker-Prager space following the Abaqus documentationEmbedded ImageEmbedded Image

Near-surface material that permits plastic deformation exhibits, in cross-section, a triangular zone of distributed plastic strain immediately above the fault tip. For the base-case South Napa model, we find that the zone of plasticity above the fault tip is relatively narrow (εplmax < 0.1 at |x| > ±10 from the fault trace at the surface) and approaches the general strain profile of a linear elastic material. This suggests that the contribution of plasticity to OFD may have been negligible and supports our use of an elastic forward model in our inversion.

In addition to the base-case model, we ran a suite of models varying cohesion from 1 Pa to 1 GPa to investigate the effect of this parameter on surface displacements (fig. S1). We find that under these loading and boundary conditions, increasing cohesion leads to a decrease in surface displacement and a broader zone of surface deformation. At c = 1 GPa, the surface displacement profile for the elastoplastic model is equal to that of a linear elastic model (that is, the yield criterion is not met, so deformation is purely elastic).

Compliant zone models. As above for the elastoplastic model, we use Abaqus to test the effect of a compliant zone on surface displacement profiles. The 3D model dimensions were 1200 m × 1200 m × 800 m. The fault dimensions were 400 m × 300 m, with the upper fault edge buried 5 m from the surface. The model used 121,090, 10-node tetrahedral elements with biased node seeding to allow for denser meshing along the fault. Node seeding ranged from 1 to 100 m. The model volume was defined by linear elasticity and partitioned such that a compliant zone within 25 m orthogonal to the fault plane could be defined with a shear modulus, G, distinct from that of the surrounding material, G0. The input material properties were as follows: Young’s modulus, E = 100 MPa; Poisson’s ratio, ν = 0.3; and density, ρ = 2000 kg/m3, which can be used to calculate shear modulus through Hooke’s law for isotropic linear elasticity G = E/[2 * (1 + n)]. Loading, boundary conditions, and contact properties were implemented in the same manner as was done for the elastoplastic models. We ran the models for the two ratios of G/G0: 0.5 and 1 (where G/G0 = 1 is the homogeneous elastic case). In fig. S1, we show only the G/G0 = 0.5 case.

Elastic inversion. In fig. S2, we demonstrate the biases in the elastic inversion by running the elastic inversion on the elastoplastic and compliant zone models with varying cohesion from above (fig. S2). For the elastoplastic models, slip S tends to be overestimated by a maximum of ~5% (fig. S2A), and upper depth zu can be underestimated by as much as 60% (recall from above that true upper depth is 5 m). For the compliant zone model, slip S is actually underestimated by ~14%, and depth is underestimated by ~7% (fig. S2).

Estimation of fault-parallel surface offset

We compare two different methodologies: (i) estimation of fault-parallel surface offset by regression of linear features and (ii) estimation of a shallow slip model. Independent of our estimation of a shallow slip model, for a given vine row, we estimate uy, fault-parallel offset, by first rotating the data set into a coordinate system where the fault trace is parallel to the vertical coordinate axis. We remove outliers using a coarse spatial mask and fit a linear model to the two portions of the vine row determined to be lying more than 10 m from the center of the fault zone. We constrain the slope of each line segment to be the same and solve for it using MatLab’s robust linear regression routine, robustfit, that uses iteratively reweighted least squares with a bisquare weighting function. uy is then estimated as the mean value of the residuals between data and models (fig. S3).

Shallow slip model and first epoch results

To minimize the number of free parameters in the shallow slip model inversion and because the analytical forward model is vertical (screw dislocation), as above, we define the X = 0 position of each vine row as the point inside the fault zone where the curvature changes from concave down (x < 0) to concave up (x > 0). In the main text, we show the inversion results for the second-measurement epoch (Fig. 3). We also show the inversion results for the first-measurement epoch (1 to 2 September) for each vineyard block (fig. S4).

Difference between slip and surface displacement estimates

In fig. S6, we demonstrate the underestimation bias E that arises from assuming that uy and S are equivalent. In the plot, each column represents the probability distribution of error, E (defined as S-uy), in 20-cm bins for dislocation upper depth, zu, grouped into 1-m bins. For both measurement epochs, the bias is present for all depth increments (Fig. 4, A and C). There is an overall bias for all of the South Napa measurements of 17 ± 30 (2σ cm for the first epoch) and 16 ± 38 (2σ cm for the second epoch) (Fig. 4, B and D). With respect to the estimated slip values, this can represent as much as ~40 to 50% underestimation bias.

OFD calculations

Two recent attempts to quantify OFD have used some combination of the difference of far-field (measured outside the zone of high shear strain near the fault trace, uff) and near-field (measured at the fault trace, unf) fault-parallel offset. Milliner et al. (21) define a percentage, OFD% = (uffunf)/uff, whereas Rockwell et al. (33) define total OFDT = uffunf. Allowing zl to go to infinity in the equation for a screw dislocation Embedded ImageEmbedded Image(1)and substituting uff = S and unf = 2uy yieldsEmbedded Image(2)Embedded Image(3)

Similarly, it can be shown thatEmbedded Image(4)

Thus, each of these OFD metrics is a function only of the distance from the fault at which near-field offset is measured (x) and the depth to the top of the fault (zu). It is independent of the slip of the fault S.

Laboratory testing

We collected three samples from the base of the south-facing wall of the South Avenue trench (see location in Figs. 1 and 2D). X-ray diffraction analysis showed that samples had similar compositions of quartz + plagioclase + montmorillonite + minor amounts of chlorite and calcite. Comparison of frictional strength to mixing law data (71) suggests that the samples are roughly 50% clay. Samples were crushed until they pass through a 100-mesh sieve (149 μm) and mixed with deionized water to form a paste. Following the procedure described in Tembe et al. (71), samples were tested in a triaxial loading apparatus at constant normal stress of 41.0 MPa and constant pore pressure of 1.0 MPa for an effective normal stress of 40.0 MPa. Forcing block pairs of Westerly granite and porous Berea sandstone were prepared from 38.1-mm-diameter cylinders that were cut at 30° to the cylinder axes. Sawcuts were surface-ground and roughened with 120-grit SiC. A 1-mm-thick layer of the gouge paste was sandwiched between the granite/sandstone forcing blocks and slipped into a 4-mm-thick polyurethane tube that isolated the sample from the silicone oil–confining fluid. The lower forcing block was low-porosity granite to minimize fluid pressure transients that might accompany variations in loading stresses during testing. Porous sandstone was chosen for the upper forcing block to assure good pore fluid communication between the gouge layer and the external pore pressure control system. A constant pore pressure of 1.0 MPa was maintained by computer control throughout each experiment. The sample assembly was attached to steel-end pieces and placed in the pressure vessel. Confining pressure was adjusted automatically to maintain constant normal stress on the fault surface as a piston was advanced against the sample column. A greased Teflon shim between the piston and lower end cap minimized resistance to lateral slip that accompanies shearing on the inclined gouge layer. Axial shortening rate was varied between 0.1, 1, and 10 μm s−1 to determine velocity dependence of shear strength of the gouge layer. Coefficient of friction (shear stress/effective normal stress) at 10-mm displacement was 0.40, 0.41, and 0.45 for the three test samples. All samples had similar steady-state velocity dependence of friction: (ab) = dμss/dlnv = +0.0026 ± 0.0004.


Supplementary material for this article is available at

fig. S1. Surface displacement profiles of elastoplastic models for varying values of cohesion C.

fig. S2. Results of running the elastic inversion on the elastoplastic models in fig. S1 expressed as a percentage of the true value.

fig. S3. Example of estimation of shallow slip model (S and zu) and fault-parallel surface offset (uy) for vine rows.

fig. S4. Inversion results for the first measurement epoch (1 to 2 September) for each vineyard block.

fig. S5. Fault zone location for each vine row measurement.

fig. S6. Comparison of estimates of near-surface slip S and fault-parallel surface displacement uy for all vineyard blocks.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank two anonymous referees for their helpful reviews as well as R. Gold and J. Hardebeck for helpful reviews of an earlier version of this manuscript. We thank D. Graves for logistical assistance. We also thank R. Forloine, K. Forloine, M. Beaulac, A. Wagner, D. Zygielbaum, C. Gurney, T. Perez, R. DeLeuze, Z. Berkowitz, E. Luse, T. Hallkovich, A. Ceja, J. Jacobs, and M. Jacobs for providing access to vineyards and private property. We thank C. Ji and D. Melgar for their kinematic inversion results. Funding: We acknowledge FEMA (Federal Emergency Management Agency) funding provided to the U.S. Geological Survey (USGS) through interagency agreement number HSFE09-15-X-0805. Paleoseismic trenching was funded in part by a grant from the USGS’s National Earthquake Hazards Reduction Program (award number G14AP00035). Author contributions: B.A.B. wrote the manuscript, collected all field data, and helped in designing and performing all analyses. S.E.M. helped in designing the elastic inversion and participated in all analyses. C.L.G. processed the LIDAR data. J.M.N. performed the elastoplastic forward models. T.D., R.R., and M.M. excavated and logged the trenches. T.L.E. and K.H. collected field data. D.L. performed the mechanical testing. V.L. contributed regional geological interpretation. A.L. performed the fracture mapping. J.M. and D.S. participated in overall interpretation. D.Z. collected field data and contributed vineyard process expertise. All authors edited the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article