Bridging earthquakes and mountain building in the Santa Cruz Mountains, CA.

Relative crustal motions along active faults generate earthquakes, and repeated earthquake cycles build mountain ranges over millions of years. However, the long-term summation of elastic, earthquake-related deformation cannot produce the deformation recorded within the rock record. Here, we provide an explanation for this discrepancy by showing that increases in strain facilitated by plastic deformation of Earth's crust during the earthquake cycle, in conjunction with isostatic deflection and erosion, transform relative fault motions that produce individual earthquakes to geologic deformations. We focus our study on the data-rich Santa Cruz Mountains, CA, USA and compare predicted and observed quantities for rock uplift, apatite (U-Th)/He thermochronology, topographic relief, 10Be-based erosion rates, and interseismic surface velocities. This approach reconciles these disparate records of mountain-building processes, allowing us to explicitly bridge decadal measures of deformation with that produced by millions of years of plate motion.

INTRODUCTION

The canonical view of the earthquake cycle posits that tectonic plate motion produces elastic straining of the crust until faults rupture in an earthquake, which unloads the crust as the cycle begins anew (). Once the cycle is complete, all nonrecoverable deformation is resolved as finite slip along faults, except in places where fault discontinuities or changes in geometry allow strain to persistently accumulate (, ). Geodesy has confirmed the essential aspects of this cycle in which interseismic elastic strain within the crust is recovered, in large part, during earthquakes (–). In contrast to this view, the million-year-long record of earthquakes and mountain building preserved in areas around plate-bounding faults often reveals large, spatially complex deformation that differs from that expected from simple elastic straining and recovery (–). Hence, reconciliation of decade-long geodetic measures of strain accumulation and release with the long-term strains recorded in geologic structures [e.g., (–)] has become a grand challenge in the field of tectonics and crustal deformation (–). To address this problem, we combined records of rock uplift, exhumation, topographic relief, erosion rates, and interseismic surface velocities to span the time scales over which mountains are built, and we compared these observed quantities with predictions from tectonic, geomorphic, and earthquake cycle models. We found that deformation of Earth’s crust observed over decades can be reconciled with observations of deformation spanning millions of years by a simple modification of the traditional earthquake cycle, such that elastic loading between earthquakes is complemented by plastic yielding of rocks adjacent to faults, which modulates stresses within the crust over many earthquake cycles. Modeling suggests that whereas uplift develops both during and in between earthquakes in spatial distributions that are distinct from one another, plastic strain, and, hence, mountain building surrounding bends in strike-slip faults, mostly accumulates between, rather than during, earthquakes.

Earthquake cycle deformation in the Santa Cruz Mountains

The Santa Cruz Mountains (SCM), situated within the San Francisco Bay Area in CA, USA, have hosted tens to hundreds of thousands of large earthquakes since the initiation of their formation roughly four million years (Ma) ago (–). Hence, topography (Fig. 1A) (), geologic structure (Fig. 1B) (–), and total rock uplift over the past 4 Ma (Fig. 1C) () in the SCM record the cumulative effects of earthquake cycle–related deformation summed over millions of years. In this region, numerous large, active strike-slip faults produce earthquakes while they accommodate the 38 mm/year of relative right-lateral motion between the Pacific Plate and the Sierra Nevada/Great Valley microplate (). The largest structure, the San Andreas fault (SAF), accounts for ~15 to 20 mm/year of this motion (, ). In the southern San Francisco Bay Area, the trend of the SAF decreases by an average of ~12° over a distance of ~80 km. Near San Jose, the SAF bisects the SCM while the fault’s trend locally decreases by as much as an additional 8° to 10° over the course of 50 km. This sharp left-bend in the SAF (the SCM bend; see Fig. 1A) converts horizontal to vertical motions along regional faults that localize contraction within the bend (Fig. 1B) (, ). The seismic risks posed by earthquakes along these faults to the >7 million residents and Silicon Valley’s globally important economy () have focused detailed geophysical and geologic studies on this area. Surface deformation (, , , ), landscape analyses (, , ), and more than a century of geologic and geophysical characterization (, , –) constrain the kinematics of the SCM over a range of spatial and temporal scales. We use these data to examine how the accumulation of earthquake cycle deformation and concomitant erosion build the topography and geologic structures adjacent to a tectonic plate boundary.

Fig. 1.

Topography, geology, rock uplift, and predicted rock uplift in the SCM.

(A) Relationship between topography of the SCM () (location delineated in inset) and the right-lateral SAF (black line). Distance northwest of the SCM bend southern terminus (red triangle, distance shown along scale bar) shows geographic extent of axes in Figs. 2 and 4 (see Materials and Methods and fig. S4). (B) Regional geology of the Pacific Block (southwest of the SAF) and Bay Block (northeast of the SAF) (–). (C) Total post–4-Ma rock uplift recorded by shallow-marine geologic units southwest of the SAF (), with white hatched line inscribing the range within which rock uplift quantities represent minimum values. Modeled total rock uplift over the last 4 Ma predicted by flexural-isostatic compensation combined with predictions from (D) a linear-elastic dislocation model whose segments slip at geodetically inferred slip rates, (E) a finite-element model with linear elastic blocks separated by a frictional contact surface, and (F) a finite-element model with strain-hardening elastoplastic block rheologies separated by a frictional contact surface. Black dashed lines in (E) and (F) represent resultant model frictional contact surface geometries upon model completion (SAF drawn in solid, bold black line for comparison).

Fig. 2.

Apatite (U-Th)/He ages (± 1𝜎 errors) in the SCM.

(A) Sample locations, with weighted mean A-He ages. Weighted mean ages for samples within the Pacific (B) and Bay (C) Blocks plotted with distance from the SCM bend southern terminus [red triangle; x axes correspond to scale bar distance in (A)]. Sample data and locations are listed in tables S2 and S3.

Fig. 4.

Elastic and elastoplastic model predictions compared to SCM observations across datasets.

Comparison of modeled (elastic, blue; strain-hardening elastoplastic, red) and observed (black) quantities measuring crustal deformation within the SCM over millions of years to decades. (A) Total rock uplift recorded by a 4-Ma datum exposed southwest of the SCM, averaged over distances within 10 km (black) of the SAF. The dashed line and arrows demarcate minimum values provided by the geologic record. (B) Observed weighted mean A-He ages ±1𝜎 uncertainties (open circles), predicted A-He ages [using model thermal paths and QTQt (); diamonds], assuming a variable geothermal gradient that initiates at 35°C/km and thermal model isothermal contours (10°C increments) spanning the typical partial retention zone for He in apatite (40° to 80°C, dashed lines) that show the evolution of model temperature through time, where 0 Ma on the y axis corresponds to present day and 4 Ma corresponds to the start of the model. We list observed and predicted A-He ages for these samples in table S4. Whereas the strain-hardening elastoplastic predictions (red diamonds) resemble observed A-He ages, the elastic model (blue diamonds) fails to produce sufficient exhumation to predict apatite (U-Th)/He ages younger than 4 Ma at the surface (blue arrows indicate that QTQt model predictions are minimum values). (C) SCM relief (within 10 km of the SAF) and (D) 10Be-based erosion rates [1𝜎 error bars; ()] measured southwest of the SAF. (E) Predicted and observed interseismic velocities from elastic (blue) and elastoplastic (red) finite element model calculations (68% confidence ellipses). All distances are measured from the southern terminus of the SCM bend (0 km on scale bar in Fig. 1A, see Fig. 3). Model comparisons northeast of the SAF are exhibited in fig. S9.

Within the SCM bend, off-fault deformation distributed across a complex plate boundary zone (i.e., along SAF subsidiary regional faults) contributes to mountain range uplift. Here, instead of explicitly modeling all faults in the SCM, we simplify the plate boundary to a single, vertical structure and allow forward-modeled crustal deformation to dictate the magnitude, timing, and distribution of crustal thickening (see Materials and Methods and Supplementary Text). In this conception, analyses of these existing data demonstrate that neither the currently observed surface deformation nor the accumulation of deformation produced by repeating elastic earthquake cycles along the SAF can explain the topography and geologic configuration of the SCM. Specifically, geologic measures of rock uplift over the past 4 Ma inferred from deformed, shallow marine Pliocene deposits preserved southwest of the SAF (Fig. 1C) (, , ) far exceed those predicted by the accumulation of deformation by repeated earthquakes at geodetically calibrated fault slip rates (Fig. 1D). Similarly, the expected long-term accumulation of strain within an elastic crust that hosts a slipping SAF produces little rock uplift and contorts the surface trace of the SAF into geometries that are not currently observed (Fig. 1E). Furthermore, neither of these scenarios, which are meant to capture the key aspects of the widely deployed elastic earthquake-cycle model, produce patterns of rock uplift and exhumation represented by the observed topography (Fig. 1A). Thus, whereas linear elasticity may adequately model deformations over an individual earthquake cycle, it is clear that the first-order characteristics of this mountain belt cannot be explained by simply summing the elastic strains expected from repeated earthquake cycles along the SAF.

Bridging disparate records of deformation within the SCM, CA

Although short-term measures of deformation [geodesy over decades (, )] are clearly incompatible with deformation inferred from observations spanning longer time scales [erosion rates over millennia () or with topographic development over tens to hundreds of millennia (, , , )], the responsible processes are obscured by ill-defined kinematic relationships pertinent to SCM formation since 4 Ma. We bridge these disparate records in the SCM using apatite (U-Th)/He thermochronology (A-He) (Fig. 2). We use these data, which approximately record the time since rocks cooled below 65° to 70°C (–), and, by inference, the patterns and rates of exhumation to constrain kinematics over time scales spanning hundreds of thousands to millions of years (). A-He ages suggest that localized deformation in the vicinity of the SCM bend has exhumed rocks from depths exceeding 2 km since 4 Ma (Fig. 2). Exhumational patterns inferred from A-He ages are consistent with the apparent lateral translation of uplifted crust along the southwestern side of the SAF as inferred from rock uplift (Fig. 1C) (), topography (Fig. 1A) (, ), and 10Be-based erosion rates (). Neither localized deformation nor the along-fault advection of irrecoverably deformed crust are features predicted by continuous, purely elastic crustal deformation in the vicinity of a fault bend (). Furthermore, the spatial asymmetry in topography (Fig. 1A) and cooling ages (Fig. 2) suggests that more complicated rheologies, whose material properties vary across the fault (texts S1 to S2), might provide a dynamical explanation that reconciles deformation observed across these different spatial and temporal scales.

To explore this idea, we developed a model that extends the earthquake cycle to include the effects of strain-hardened plastic yielding (, ), isostatic and flexural compensation (–), and surface erosion (–). We used the Abaqus finite element modeling suite (), a commercial engineering software package designed to simulate the mechanical behavior of materials under prescribed loads and model boundary conditions, to calculate long-term crustal thickening histories and current thickening rates induced by the SCM bend. We built a simple model geometry in which a frictional SAF that separates two deformable crustal blocks hosts a bend along its length (fig. S1). This geometric framework assumes that crustal deformation and thickening in the model encapsulates both coseismic slip along regional faults and distributed off-fault deformation in the SCM (text S1). We first use a purely linear elastic rheology, which provides long-term predictions of deformation due to repeated elastic earthquake cycles along the SAF. We then consider an elastoplastic crustal rheology, which allows for the development of irrecoverable off-fault deformation while coherent block motion accrues across the SAF. In strike-slip fault systems, crust surrounding geometric asperities, such as fault bends, may strain harden despite demonstrable evidence of deformation localization due to the persistent production of faults and fractures across a variety of spatial scales [text S2 and (, )]. In prescribing these strain-hardening elastoplastic rheologies, we assert that the fold-dominated Pacific Block (see Fig. 1B), which largely consists of thick, multilayered sedimentary rocks (–), moderately hardens while it advects into and through the SCM bend. In contrast, we suggest that the fault-dominated Bay Block—which has remained fairly immobile relative to the SCM bend (, ) and consists of sedimentary, metasedimentary, and ophiolitic rocks (, )—hardens rapidly due to the fact that existing regional thrust faults rotate out of favorable planes for slip as new faults develop with accumulating deformation (). We prescribe these rheological differences across the SAF to capture the role that contrasting lithologies, geologic histories, and crustal structure on either side of the SAF appear to play in the development of deformation in the SCM.

Thickening of crust in the vicinity of the SCM bend (fig. S2) is driven by SAF rate shear velocities (16.8 mm/year) applied to the edges of the model, which are communicated to the modeled frictional SAF via a shear-traction–free base located at 20 km in depth (vertical velocities along this surface are required to be zero) (fig. S1) [e.g., ()]. The thickening history is then used with a simple isostatic flexural model (, ) to calculate both the predicted 4-Ma-long history and instantaneous rates of rock uplift [flexurally modulated, uneroded rock uplift grids are shown in Fig. 1 (E and F)]. Because crustal thickening, rock uplift, and surface elevations depend on the history of erosion of Earth’s surface (, ), we used the power-law incision rule (–) to model erosion rates and surface elevations over time and space (Fig. 3). Last, we coupled the history of erosion with a simple thermal model in which the active geothermal gradient evolves from an initial 35°C/km state (, ) in response to exhumation-driven one-dimensional (1D) heat advection and conduction to model the cooling history of rocks as they are exhumed to the surface. We then input these cooling histories into the QTQt kinetics model (text S3) (), which models the temperature-dependent diffusion of helium in apatite through time, to predict A-He ages at individual points that have been exhumed from depth and are currently exposed at Earth’s surface. Together, the combination of these models predicts the history of rock uplift, rock uplift rates, erosion rates, surface uplift, topography, and A-He cooling ages throughout the mountain range. We compared these predictions with the multifold data from the SCM to determine whether simple modifications of the rheology of the crust used in earthquake cycle models might reconcile deformation observed across disparate time scales spanned by different geologic and geophysical measurements (Fig. 4).

Fig. 3.

Elastoplastic tectono-geomorphic model evolution through time.

The temporal evolution of equally spaced river catchments i to v (location and down-stream orientation delineated by white arrows) spanning the length of the SCM bend southwest of the SAF highlight changes in rock uplift, relief, area-weighted mean catchment elevation, and exhumation through the course of SCM formation since 4 Ma for the strain-hardening elastoplastic model. (A) The areas of analysis that encompass the rock uplift, topographic, A-He and 10Be datasets we consider extend 10 km from the SAF in map view (20 km in stream profile distance; see Materials and Methods), and are delineated by the rounded polygons outlined in solid black (Pacific Block; Fig. 4) and dashed gray (Bay Block; fig. S9). Bend-relative catchment distances within the Pacific Block in (B) are measured relative to the southern terminus of the SCM bend (red triangle) and are measured in a reference frame perpendicular to the SAF northwest of the SCM bend (orientation depicted by blue arrow). Each catchment model illustrates the evolution of predicted stream profile topography from 4 Ma (when our model begins) to 0 Ma (present day). Our model subsumes the Black Mountain bend into the larger SCM bend, so the effects of this structural feature are not captured in our model results (see Discussion).

Investigators have used elastoplastic rheologies to model crustal deformation over a variety of temporal and spatial scales (–), although the manner by which plastic strain accumulates during individual, predominantly elastic earthquake cycles remains unclear. We examined temporal variations in the accrual of uplift and plastic strain during an individual earthquake cycle along the SAF in the deformed elastoplastic model. In this procedure, we allowed the length of the modeled SAF to rupture, analogous to the moment magnitude (M) 7.9 San Francisco earthquake in 1906 (). We then extracted horizontal interseismic velocities late in the earthquake cycle along the SAF and superimposed them with those produced by slip along the four other major fault systems within the area. We reference these velocities to the LUTZ GPS station to provide a direct comparison with geodetic observations (Fig. 4E) (). Last, we analyzed distributions of uplift (Fig. 5) and plastic strain (fig. S8) during the coseismic and interseismic phases of the earthquake cycle to assess how predominantly elastic deformations, incurred during earthquakes and interearthquake periods, combine to build mountain ranges over geologic time scales.

Fig. 5.

Seismic cycle rock uplift and rock uplift rates in the SCM.

Modeled quantities for (A) coseismic rock uplift, (B) interseismic rock uplift rates, and (C) integrated seismic cycle (coseismic + interseismic) rock uplift rates show the distribution of rock uplift and subsidence in the SCM during an earthquake cycle along the SAF for the elastoplastic model. Comparison between the integrated seismic cycle rock uplift rate field and the (D) present-day rock uplift rate field active after 4 Ma of deformation in the long-term elastoplastic model reveals that these rates are approximately equal.

RESULTS

Apatite (U-Th)/He thermochronology

We collected and analyzed 15 sedimentary bedrock samples: six from the southwest of the SAF and nine from the northeast (Fig. 2 and table S2). A-He ages observed in these samples reveal three fundamental kinematic properties of the range that informed our tectonic thickening model design. First, southwest of the fault, a monotonic northwestward increase in A-He ages with distance from the southern tip of the SCM bend suggests that this left bend has caused rapid uplift and exhumation in its vicinity while rocks were translated relatively northwestward through the structural knot in the fault (, , , ). Second, the lack of correlation between cooling ages and along-fault distance from the SCM bend northeast of the SAF suggests that deformation in this zone has persisted within the bend without relative translation of this material to the southeast. This “trapping” of crust within the northeastern block bounded by the SCM bend appears to produce asymmetric topography and exhumation across the fault in which doming northeast of the SAF contrasts with northwest-elongated, roughly cylindrical folding of crust to the southwest (, , ). Third, post–4-Ma deformation to the southwest of the SAF suggests that deformation has localized in a 5- to 10-km-wide zone adjacent to the SAF over the course of SCM formation. These ages, their trends, and the kinematic relationships they reveal are largely consistent with previous conceptions of SCM formation (, ).

Reconciliation of disparate datasets

Using our model framework, we track the coevolution of rock uplift, exhumation, topographic relief, and erosion in response to the crustal thickening derived from the mechanical model. To do this, we model roughly fault-perpendicular river catchments extending 10 km away from the SAF (20 km in profile length; see Materials and Methods) and quantify the evolution of these river catchments in response to tectonic thickening through model time (Fig. 3). Modeled catchments spanning the length of the SCM bend (e.g., profiles i to v) record changes in isostatically compensated rock uplift, exhumation, topographic relief, and catchment-averaged erosion rates since 4 Ma, which we use to make predictions for observed quantities in the SCM. Comparisons of modeled (blue, linear elastic; red, strain-hardening elastoplastic rheologies) and observed (black) total post–4-Ma rock uplift (Fig. 4A), A-He cooling ages (Fig. 4B), topographic relief (Fig. 4C), and modern, 10Be-based erosion rates (Fig. 4D) show that observations are generally predicted by the strain-hardening elastoplastic models, but not the linear elastic rheologies. In the latter case, modeled rock uplift values underpredict observations by >1 km, and modeled exhumation is insufficient to produce cooling ages younger than 4 Ma. Similarly, topographic relief is underpredicted by the elastic model, and modeled erosion rates are <½ of those currently observed. In contrast, strain-hardening elastoplastic models predict rock uplift magnitudes consistent with the minimum bounds imposed by observations. These models also predict the magnitudes and trends in observed A-He ages and modern erosion rates, although modeled relief is larger than that observed in the bend and is less than that observed northeast of the bend (see Discussion). The ~7× increase in crustal strain facilitated by plastic deformation (relative to elastic model strains; fig. S3) since 4 Ma drives the observed differences between model predictions. Last, the asymmetry in the topography and long-term deformation within the bend could only be reproduced using material properties that varied across the modeled SAF for both elastic and strain-hardening elastoplastic models (figs. S4, S5, S16, and S17, and texts S2 and S5). Deformation that develops over the course of the model causes a progressive migration of the SAF location and a decrease in the bend angle with model time in both cases; the linear-elastic model predicts a contorted SAF geometry with a bend angle that is gentler than that currently observed in the SCM (Fig. 1E). In contrast, differing cross-fault strain-hardening properties (fig. S4) in the elastoplastic model localize deformation within the bend, which ultimately produces a bend geometry similar to that observed (Fig. 1F and animation S1).

Long-term elastic earthquake cycle and elastoplastic model predictions diverge, although it is unclear whether these differences are expressed within geodetic observations measured over decades. We examine the impact of elastic and strain-hardening elastoplastic rheologies on interseismic velocities by comparing extracted velocities for each case with regional observations (Fig. 4E) (). Observed horizontal surface velocities are reproduced by both the linear-elastic and strain-hardening elastoplastic scenarios, which produce similar goodness-of-fit measures to the observed velocity field (table S1 and text S6). Whereas geodetically measured horizontal surface velocities cannot discriminate between these two first-order rheologies, peak vertical interseismic surface velocities predicted by the elastoplastic model locally exceed those predicted by the elastic model by >300% (fig. S6). This rate discrepancy mimics the long-term record, in that the rock uplift predicted by the purely elastic, SAF earthquake-cycle model accounts for only ~30% of that produced in the elastoplastic model (fig. S6).

Rock uplift and plastic strain development over an earthquake cycle

The long-term elastoplastic model reveals how deformation develops over disparate scales of time and space, although within this model framework, the role of individual earthquakes and interearthquake periods in dictating long-term deformation patterns remains enigmatic. To predict the distribution and magnitude of deformation that develops during and in between earthquakes, we model a 250-year-long interseismic period around the frictionally locked SAF, followed by an individual earthquake characterized by slip associated with a reduction in friction along the modeled fault to a depth of 10 km after the final time step of the long-term model (see Materials and Methods and fig. S7). Earthquakes in the SCM (Fig. 5A) generate rock uplift concentrated near the southern tip of the SCM bend southwest of the SAF and adjacent to the northern tip of the SCM bend northeast of the SAF. The tallest peaks on either side of the SAF subside during modeled earthquakes. Interseismically (Fig. 5B), broad, diffuse rock uplift persists about the ends of the bend, while plate motion deforms Earth’s crust around a frictionally locked fault. We add the components of rock uplift associated with both the coseismic and interseismic phases of the seismic cycle to image the cumulative rock uplift-rate field over an earthquake cycle (Fig. 5C). We find that rock uplift is concentrated near the southern end of the bend within the crust southwest of the SAF, whereas rock uplift about the northern end, as well as east of the SCM bend northeast of the SAF, is more diffuse. Single earthquake cycle rock uplift rates are consistent with rates derived from the long-term model (Fig. 5D), and the observed distribution and magnitude of rock uplift persist over multiple earthquake cycles (animation S2). Last, whereas coseismic and interseismic rock uplift contributes to rock and surface uplift within the SCM over time, most of the plastic strain accumulates between earthquakes (fig. S8).

DISCUSSION

Numerous studies have used complex nonlinear rheologies to model deformation in Earth’s crust over disparate time scales (, , ), but none have attempted to reconcile records of deformation that capture geodetic, geomorphic, thermochronologic, and geologic responses to deformation that bridge time scales spanning decades to millions of years, as we present in this contribution. Model results indicate that plastic strains mostly develop interseismically, suggesting that it is largely the steady motions between earthquakes along the SAF that build the SCM rather than deformation that develops during or immediately following large earthquakes themselves (–). Postseismic relaxation from viscous flow in the lower crust and upper mantle may also produce minor, long-wavelength deformation in the SCM [e.g., )], although we do not model these effects. In addition, our model does not incorporate dynamic, time-dependent earthquake rupture processes or aseismic afterslip, which could serve to localize irrecoverable deformation adjacent to the SAF (). Nevertheless, the patterns of rock uplift produced over the course of this idealized earthquake cycle mimic the deformation accrued over many millions of years in which asymmetric long-term rock uplift develops around the SCM bend. This behavior arises because both are manifestations of the underlying crustal rheology—strain-hardening dictates the spatial distribution of plastic strain accumulation and rock uplift during the earthquake cycle. For example, the concentration of persistent rock uplift and plastic strain development near the southern SCM bend tip (Fig. 5, C and D) results from the juxtaposition of strong, strain-hardened crust northeast of the SAF with relatively weaker crust to the southwest (fig. S22). This zone of pronounced rock uplift is coincident with sharp increases in rock uplift, young A-He ages, the onset of elevated topographic relief, and high erosion rates observed in the SCM (Fig. 4). Parcels of crust that entered the SCM bend deformed, uplifted, and strain-hardened as they advected along the SAF (animation S1). Yet today, they no longer reside within the zones of concentrated rock uplift and plastic strain accumulation. These model characteristics appear required to explain the long-term kinematics proposed for SCM formation (, , ) in which crust deforms and uplifts as it enters the SCM bend and advects laterally while plate motion continues. Thus, strain-hardening induced by strain accumulation between individual earthquakes dictates the long-term kinematic properties of the SCM over geologic time scales.

Our results suggest that the geometry of the SAF evolved since the SCM bend formation and was strongly influenced by the patterns of deformation in strain-hardening crust that developed in the SCM (fig. S22 and animation S1). We assert that the SCM bend likely formed as the result of plate motion reorganization from 4 to 6 Ma (, ), which caused a shift in the location of the predominant plate-bounding fault and subsequent abandonment of fault segments poorly oriented with the new direction of plate motion [i.e., the Pilarcitos fault; (, )]. At ~4 Ma (when our model begins), deformation began to accumulate surrounding the newly formed and initially steep-angled SCM bend (fig. S1 and animation S1). Over time, the bend angle decreased due to the steady migration of the SAF in response to deformation within both the Pacific and Bay Blocks. However, the Bay Block hardened particularly rapidly due to (i) concentrated deformation focused around the steep northern terminus of the SCM bend (θ2; fig. S1) and (ii) enhanced hardening, which may result from persistent distributed thrust faulting and block rotation. The eventual juxtaposition of the strong, rapidly strain-hardened Bay Block and relatively weaker Pacific Block crust near the southern terminus of the SCM bend limited additional migration of the SAF, fixing the SCM bend angle into the currently observed geometry (animation S1).

In the SCM, the accumulation of off-SAF crustal deformation in both the Pacific and Bay Blocks has been accommodated, at least in part, along numerous mapped regional faults and folds [Fig. 1B; (, , )]. Because the strain-hardening elastoplastic continuum crustal model we use does not allow for localized deformation along structural discontinuities (text S1 and fig. S1), we cannot explicitly predict and compare the respective roles of irrecoverable deformation that is broadly distributed and that which is accommodated along discrete faults. However, the steeply rising topography west of Santa Clara Valley (Fig. 3A), the spatial discontinuity of uplifted and preserved Plio-Pleistocene marine units across faulted blocks in the Pacific Block (Fig. 1B), and the distinct increase in rock uplift adjacent to the SAF (Fig. 1C) all suggest that localized deformation along regional faults and within fault-bounded blocks strongly influenced the evolution of crustal deformation and rock uplift in the SCM. The prescribed rheological properties of the modeled Pacific and Bay Blocks attempt to capture the aggregate behavior of complex, localized features, such as faults and folds. (text S2). This approach appears to have been successful, as the first-order attributes of the observed deformation field and the associated geomorphic responses appear to be generally reproduced by the elastoplastic model predictions.

Some of these findings may be applicable to other large restraining bends, such as the “Big Bend” in the SAF, near Los Angeles, CA and the Denali bend, AK. Restraining bends evolve in their mechanical efficiency through time while deformation accrues (), but strain hardening of crust acts to stabilize restraining bends because evolving bend geometries anchor within hardened crust that is resistant to additional deformation, which allows for long-term kinematic behaviors to develop and persist through time. Subtle differences in the crust’s material properties and strained state before restraining bend deformation may cause the rheology of crust on opposite sides of the strike-slip fault to evolve differently, ultimately producing asymmetries in the long-term and earthquake cycle deformations in these areas. These differences may be expressed as asymmetric topography and structural history around restraining bends (), linking the steady accrual of deformation over millions of years to the topographic expression and earthquake cycle dynamics through the evolution of constitutive behavior while crustal strain accrues.

Our model generally predicts observations of rock uplift, exhumation, landscape evolution, and erosion in the SCM, although there are some discrepancies between predictions and observations. First, modeled relief and erosion rates underpredict observations for bend-relative distances ranging from 40 to 50 km (Fig. 4), which are spatially coincident with a small, unmodeled restraining bend structure that resides ~45 km northwest of the SCM bend (Fig. 3; ()). Second, our model erroneously predicts that the A-He sample closest to the southern tip of the SCM bend was never exhumed from sufficient depths to expose fully reset samples at the surface (Fig. 4B). This underprediction of cooling may result from our use of a 1D thermal model, which does not allow lateral heat conduction that might elevate geothermal gradients in the upstream direction of the bend. Alternatively, the discrepancy may reflect the model idealization of crust as a purely plastic continuum in which strain does not explicitly localize along crustal discontinuities, such as faults. There are many such mapped faults in the SCM (Fig. 1B) (, , , ), which could serve to localize deformation and enhance exhumation experienced by the southernmost A-He sample, reducing the predicted age.

Our results demonstrate that disparate observations of crustal deformation spanning five orders of magnitude in time (decades to millions of years) and seven orders of magnitude in space (millimeters to tens of kilometers) can be reconciled by extending the earthquake cycle model to include plastic yielding with strain hardening [e.g., ()]. This work addresses a long-standing grand challenge in solid-earth science (–) and provides a dynamical framework that may be generally deployed to understand crustal deformation over vastly different time scales. We find that observed contemporary geodetic horizontal surface velocities are consistent with a wide range of crustal rheologies, and so, these measurements provide limited information about the long-term behavior of Earth’s crust. Predicted interseismic vertical surface velocity fields differ between the two rheological models (fig. S6), but current geodetic measurements of rock and surface uplift over decades include the effects of numerous nontectonic processes and postseismic relaxation from recent earthquakes, which obfuscate the spatially continuous interseismic rock uplift signal (, ). Instead, the integration of geologic, thermochronologic, geomorphic, and geodetic data in an internally consistent dynamic framework reveals the rheologies that are capable of bridging earthquakes and mountain building in the SCM.

MATERIALS AND METHODS

Apatite (U-Th)/He thermochronology—measurement procedures and data analysis

We collected 15 sedimentary and metasedimentary bedrock samples throughout the SCM and extracted apatite crystals at Stanford University. Bay Block samples were analyzed at the Stanford University Noble Gas Laboratory in 2009 following apatite (U-Th)/He measurement procedures described in (). Pacific Block samples were analyzed in the Noble Gas Thermochronometry Laboratory of Berkeley Geochronology Center from 2015 to 2018 following procedures described in (). We applied an α-ejection correction (FT) after () to all measured single-crystal ages assuming no U-Th zonation. After data acquisition, we used the robust z-score method () to identify individual FT-corrected age outliers and excluded crystals with z scores > 5 from sample populations (Supplementary Methods). We then used A-He age uncertainties to calculate weighted mean A-He ages for each sample. We show the significance of crystal rejection and weighted mean sample calculations in fig. S10 and explore the impact of crystal size and parent nuclide (U, Th, and Sm) concentrations on observed age distributions in fig. S11 and text S3. We provide laboratory measurements and sample locations in tables S2 and S3.

Dislocation model of rock uplift in the SCM

We modeled dislocation-derived crustal displacements by using geodetically derived SAF slip rates () to simulate the long-term rock uplift expected about the SCM bend using triangular dislocation elements (). We first calculated the multiseismic cycle rock uplift rates as a result of SAF slip. We incorporated the Marin, San Francisco, Peninsula, SCM1, SCM2, San Juan Bautista, and Creep1, Creep2, and Creep3 sections, as previously defined (), and applied specified slip values for each segment across the entire fault surface (we prescribed neither locking depths nor depth variations in slip/creep). We integrated these rock uplift rates over 4 Ma to produce the dislocation-derived, time-integrated crustal thickening in the SCM (fig. S2).

Spanning disparate records of deformation with a coupled tectonic-geomorphic model

From 3D models to along-fault predictions—casting results in the SCM bend coordinate system

SAF-parallel trends in rock uplift magnitudes, cooling ages, relief, and erosion rates suggest that these quantities evolve as a function of distance (and time, given the strike-slip nature of the SAF) from the SCM bend southern terminus. We compare measurements for these four quantities to modeled values using a coupled tectonic-geomorphic model framework. In this framework, we first modeled 3D tectonic thickening of Earth’s crust surrounding the SCM bend using the Abaqus software suite (described below) (). We extracted along-fault thickening rate profiles (i.e., surface displacement profiles in a fixed-base model of the upper crust; fig. S12) along the modeled SAF through time and used these to calculate isostatically compensated, along-fault rock-uplift rate profiles. Next, we calculated topographic evolution in response to input rock uplift and the associated erosional response along the SCM bend with a series of 1D river catchment models oriented perpendicular to the SAF. Last, we extracted predicted cumulative rock uplift, exhumational histories and associated temperature histories, mountain range relief, and erosion rates from the erosional models to compare against SCM measured quantities (fig. S13). In this comparison, we collapsed both measured and model-predicted quantities onto a single, fault-parallel profile that encapsulates the 10-km region adjacent to the SAF. We use this profile to cast the relevant SCM measurements and predicted quantities in a bend-relative coordinate system oriented parallel to the SAF north of the SCM bend (azimuth ~323.8°). We fixed the origin of this coordinate system to the current southern terminus of the bend, with bend-relative distance increasing to the northwest (Fig. 3). We use this coordinate system to compare modeled and measured quantities for both the Pacific (Fig. 4) and Bay Blocks (fig. S9).

Modeling crustal deformation with finite element analysis in Abaqus

We modeled SCM bend-induced crustal thickening for elastic and elastoplastic crustal rheologies using the Abaqus finite element modeling suite (). We juxtaposed two solid 3D blocks (analogous to the Pacific and Bay Blocks) across a vertical frictional contact surface [coefficient of friction = 0.2, analogous to weak portions of the SAF (–)] that consists of two parallel but offset segments connected by multiple obliquely trending segments (fig. S1 and Supplementary Methods). The contact surface geometry evolved during the model run while deformation accrued. The prescribed initial bend geometry in the elastic model mimics the current SCM bend; we used the San Francisco, Peninsula, SCM1, and part of the SCM2 fault segment locations [as previously defined () and used in the dislocation model] in conjunction with a stepover width of 16.7 km to designate the SAF surface trace and fault geometry. For the elastoplastic model, we connected these parallel, offset fault segments (offset = 16.7 km, fig. S1) with two obliquely oriented fault segments. The southern segment (fig. S1, segment F1) deviates from the parallel fault orientation by 12.5° (θ1), whereas the northern segment (F2) deviates from the parallel fault orientation by 25° (θ2). Segment F2, which was intended to mimic the steep SCM bend, initiated with a steepened bend angle of 25°, because deformation within the adjoining crustal blocks reduced the obliquity of the connecting segment over model time. All models extend to a depth of 20 km to capture the deformation resolved in the SCM upper crust (). We explored a range of material/geometric combinations for both elastic and elastoplastic models, which are documented in tables S5 to S7 and fig. S5. We extended the lateral model boundaries to >120 km from the bend, >2 times the extent of the deformed region, to minimize the influence of edge effects on our results (geometry shown in fig. S1). Motion within the model was excited by displacing the model boundaries at a total rate of 16.8 mm/year in the direction parallel to relative plate motion between the two simulated crustal blocks for the 4-Ma model duration. This condition is consistent with inferred fault slip rates across the San Andreas in the SCM () and with the ~4 Ma age of rock and surface uplift initiation in the SCM inferred from depositional ages of uplifted Plio-Pleistocene submarine units (, ) and thermochronologic cooling age distributions (this study). We applied fault-parallel displacements along the lateral boundaries of the model and held the boundaries fixed in the vertical and fault-perpendicular directions. At the model base, we allowed for unconstrained, traction-free lateral motion and prescribed a vertical zero-displacement boundary condition. We include the Abaqus input files for all models in the data repository.

Next, we prescribed linearly elastic and strain-hardening elastoplastic materials to the solid 3D blocks. Linear elastic model materials followed Hooke’s law in which strain is proportional to stress. We assumed a Young’s modulus of 30 GPa and a Poisson’s ratio of 0.25, which were isotropic and homogeneous throughout the model (tables S6 and S7). For the elastoplastic material, we additionally assumed that plastic material behavior follows a linear Drucker-Prager yield criterion (Supplementary Methods) (, ). For both blocks, the material follows a nonassociated flow rule and is nondilatant, and yield stresses in tension and compression are equal. We assigned a Drucker-Prager friction angle of β = 45° and a Drucker-Prager cohesion of d = 50 MPa for the elastoplastic material, which roughly correspond to Mohr-Coulomb equivalent values of 35° and 35 MPa, respectively (Supplementary Methods) (). These material parameter prescriptions are consistent with measured quantities for rock types representative of SCM geology (tables S8 and S9). Materials assigned to the modeled Pacific and Bay Blocks are initially mechanically identical (tables S6 and S7). However, the strengths of these materials, manipulated as a function of Drucker-Prager cohesion, may vary as strain accrues. We used observed differences in discrete deformation styles between the Pacific and Bay Blocks (i.e., pervasive cylindrical folding vs. extensive thrust faulting) in conjunction with inferred kinematic characteristics of the SCM to inform our bulk crustal material property prescriptions. We explored the impact of different user-specified relationships between cohesion and accumulated plastic strain, which we allowed to vary between blocks, on the divergent deformation histories predicted for each side of the simulated SCM bend. We detail these strain hardening behaviors in fig. S5 and discuss prescription justifications in text S2.

We assumed that deformation-induced vertical nodal displacements generated in the mechanical model capture the timing and extent of crustal thickening for the prescribed model scenario (text S4). Following model completion, we extracted bend-relative nodal displacement histories, which we used to approximate tectonic thickening rates, on the model surface (depth = 0 km) within 10 km of the contact surface for both solid blocks. This area best represents the zone of localized deformation adjacent to the SAF captured in the rock uplift field (Fig. 1C), and the zone from which our cooling ages (Fig. 2) and erosion rates (Fig. 4D) () were collected. We used a 5-km-wide, fault parallel moving window to extract averaged nodal displacements as a function of distance along the fault relative to a specified southern bend terminus (unique to each model) at regularly spaced, 250-thousand year intervals. We used vertical displacements to capture crustal thickening and fault-parallel displacements to resolve the rates at which crust advects through the modeled SCM bend (fig. S12). We integrated the space-time history of advection rates to track the time-space history of inferred crustal thickening throughout the near-fault model domain. We used these histories to model the topographic and erosional history of each point along the crest of the range that is currently lying at the surface (discussed below).

Calculating isostatic deflection and rock uplift for a thickened crust in the SCM

We approximated the isostatic response to deformations calculated within the Abaqus models using a simple isostatic compensation model to estimate the rock- and surface-uplift response to thickening (see previous section) and erosion (see following section) of the modeled SCM. To calculate isostatic deflection and resulting rock uplift due to crustal thickening, we first extracted vertical surface displacements (rock uplift in the absence of isostatic deflection and erosion, which we equate to crustal thickening) predicted by each model. This crustal thickening rate was then adjusted for the mass removed by erosion of the surface (see following section). Using this net thickening rate, we then calculated the elastic crustal deflection by using the simple approach outlined in () and exploited in () in which the isostatic flexural response varies with the wavelength of the mean topography. We used the parameter values we prescribed in Abaqus (Young’s modulus = 30 GPa, Poisson’s ratio = 0.25, and elastic crustal thickness = 20 km) and assigned a topographic wavelength of 60 km or twice the 30 km average width of the SCM to calculate the ratio (f) of a combined flexural-isostatic deflection to a pure Airy isostatic crustal response (see Supplementary Methods). These parameter prescriptions yield f ~ 0.19. However, the effective value of f doubles in the presence of a “broken” elastic plate (i.e., along plate boundaries, as is the case in the SCM) (), yielding an effective f ~ 0.38. This f corresponds to the lower bound of elastic crustal properties (E = 10 GPa; v = 0.25 and h = 5 km for an intact elastic plate) previously assumed for the SCM (text S4) (). With f in hand, we calculated rock uplift by subtracting the resultant deflection from crustal thickening inputs.

Rock uplift grids shown in Fig. 1 (D to F) were calculated directly from the final state of tectonic thickening resolved in the time-integrated dislocation model (Fig. 1D) or in the Abaqus models (Fig. 1, E and F) and do not include the effects of erosion. These calculations use a Fourier transform technique (see Supplementary Methods) to calculate isostatic deflection and resultant rock uplift in response to tectonic thickening for an intact elastic plate ().

Modeling the geomorphic response to dynamic, isostatically compensated rock uplift

The isostatic-flexural response and the exhumation history recorded by the A-He thermochronology within the SCM critically depend on the removal of mass from its surface. We designed a coupled tectonic-geomorphic model framework to simulate the evolution of the SCM since 4 Ma using the power-law incision rule (–). We used area-length relationships that we calculated for the SCM and specified an erodibility (K) of 10–4.6 m0.2/year [comparable to volcaniclastic rock erodibilities; ()]. We assigned a model catchment flow length of 20 km to simulate the evolution of the SCM within 10 km of the SAF, because Hack scaling suggests that catchment flow lengths are greater than corresponding map view catchment lengths by roughly a factor of 2 (). We summarize these geomorphic parameter prescriptions and their calculations in table S10 and in Supplementary Methods.

We began by initiating all catchment model elevations at zero and fixing the catchment position furthest from the SAF (x = 20 km) at base level (zero elevation) for the duration of the geomorphic model. Then, we used the interpolated advection and thickening rate surface extracted from Abaqus to calculate bend-relative displacement histories and corresponding thickening histories for each catchment since 4 Ma. We used catchment crustal thickening histories to drive rock and surface uplift in the geomorphic models in which we account for the role of crustal deflection in compensating crustal thickening to produce isostatically compensated rock uplift (above). We first calculated the area-averaged erosion rate for the catchment in its current state. We then calculate the rock uplift rate, Urock [L/t], in response to both the reduction of the crustal thickening component (𝜏) through erosion and the influence of isostatic deflection aswhere ⍴m [M/L3] and ⍴c [M/L3] are the mantle and crustal densities, respectively, f [] is the ratio of a combined flexural-isostatic deflection to a pure Airy isostatic crustal response (), and [L/t] is the area-averaged erosion rate. We assume a mantle density of 3200 kg/m3, a crustal density of 2600 kg/m3, and that f = 0.38. We used the resultant rock uplift rates to calculate the changes in topography in response to rock uplift and erosion following (), which yielded topography as a function of time along the modeled river catchments within the SCM bend (Fig. 3 and Supplementary Methods). We show the impact of variable prescriptions of K, f, and crustal thickening histories on modeled results in figs. S14 to S17.

Comparing modeled and measured quantities in the SCM

With this model framework, we calculated the total rock uplift in each catchment; the evolution of surface uplift, topography, and relief as a function of time; and catchment-averaged erosion rates throughout the model. We extracted total rock uplift by integrating the calculated rates over the duration of the geomorphic model. We compared these quantities to averaged rock uplift values within 10 km of the SAF, which were computed from rock uplift magnitudes inferred from deformed Plio-Pleistocene marine units preserved in the SCM (, , ) using a circular, focal mean metric. We compared topographic relief by subtracting the highest point in the catchment model (the crest elevation) from the lowest point in the model (which we designate as zero) and comparing these quantities to topographic relief in the SCM calculated using a focal-range metric within 10 km of the SAF. We calculated current catchment-averaged erosion rates derived from the final catchment profile and compared these to erosion rates reported in the SCM (). Last, we calculated exhumational histories by subtracting the area-averaged mean elevation of each catchment from the corresponding cumulative rock uplift at each model time step and converted these to time-temperature paths and apatite (U-Th)/He ages (described below).

We converted modeled exhumational histories to time-temperature histories to calculate the suite of apatite (U-Th)/He ages predicted by each model configuration. We began by assuming an initial, spatially constant geothermal gradient of 35°C/km at the beginning of the model at 4 Ma. For each bend-relative position, this gradient evolves in time as a function of 1D vertical heat advection in response to erosion and exhumation calculated by the coupled tectonic-geomorphic model (fig. S18). The SCM average surface temperature since 4 Ma was set to 15°C, and we assumed a constant thermal diffusivity of 40 m2/year [L2/t], as is consistent with sedimentary rock types found in the SCM (). The temperature history of a given catchment model at its mean elevation can thus be expressed aswhere i represents model time and ranges from 0 to 4 Ma, T [T] is the temperature at a given model time, Tsurface [T] is the surface temperature, Ztotal [L] is the total exhumation experienced, Z [L] is the exhumation experienced at the given model time, and ∇geothermal( [T/L] is the geothermal gradient at the given model time. Using this procedure, we extracted thermal histories that correspond to the A-He dataset sample locations and used these to forward-model predicted cooling ages with the QTQt kinetic model (). These models used the input time-temperature history in conjunction with measured crystal dimensions to calculate He diffusivity, He retention, and predicted cooling ages of individual samples currently exposed at Earth’s surface. Because the time-temperature histories only capture the thermal evolution of the modeled crystals since 4 Ma, they provide an incomplete thermal record for the samples we consider. However, although long-term, complex thermal histories can substantially affect resultant cooling ages (), we chose to model only the timeframe that is specifically predicted in our coupled tectonic-geomorphic model (text S3). We used QTQt to calculate the mean corrected age for each sample for comparison with the measured cooling ages (Fig. 4B and table S4).

Seismic cycle modeling in the SCM

Modeling coseismic and interseismic phases of the seismic cycle

We modeled coseismic and interseismic phases of the seismic cycle using our Abaqus model in its final, deformed state. To do this, we first applied 2.1 m of additional right-lateral displacement to each of the lateral model boundaries following the completion of the 4 Ma model duration for a total slip of 4.2 m, as would be consistent with the geodetically inferred slip rate of 16.8 mm/year applied over 250 years. Next, we extracted the nodal velocities for the model base, modeled fault surfaces, and lateral model boundaries during this 250-year period, which we treated as a proxy for long-term crustal velocities integrated over multiple seismic cycles at the end of 4-Ma model time. With this information, we modeled the phases of the seismic cycle as follows. For interseismic periods, we first increased the frictional coefficient of the upper ~10 km of the SAF from 0.2 (the long-term frictional coefficient we use) to 0.3 [reported SCM locking depths range from 9 () to 15 km ()]. We then applied the long-term nodal velocities to the lateral and basal model boundaries and to the lower ~10 km of the SAF to simulate the continuation of plate motion in the presence of a frictionally locked upper fault. We applied these velocities for 250 years to simulate a 250-year-long earthquake cycle, which is broadly consistent with the inferred recurrence interval for large earthquakes along the SAF (). During this time, the upper portion of the SAF remained frictionally locked due to the increased frictional coefficient. To model the coseismic phase, we returned the frictional coefficient along the upper ~10 km to its inferred long-term value by decreasing the frictional coefficient from 0.3 back to 0.2. We fixed nodal positions of the lateral and basal model boundaries and the lower ~10 km of the SAF during the coseismic phase. The reduction in friction along the upper portion of the SAF in conjunction with interseismically accumulated displacements and associated stresses drove the coseismic fault rupture along the upper portion of the SAF. Last, we used the vertical displacement fields associated with each of these phases to assess interseismic and coseismic rock uplift distributions. We graphically summarize this process in fig. S7.

Calculating interseismic velocity fields from deformed Abaqus models

We extracted horizontal surface velocities during the interseismic phase of the seismic cycle in our Abaqus model to produce surface velocity fields that can be directly compared to interseismic velocities in the SCM. We used the Abaqus-derived interseismic velocity field to interpolate velocities at each geodetic station in the BAVU network (). We calculated all station velocities relative to a fixed LUTZ geodetic station (in San Jose, CA) and subsequently added velocity components resulting from movement along other major strike-slip faults in the San Francisco Bay Area to produce interseismic velocities comparable to those reported for the SCM (). In this work, we exclusively considered the set of stations bounded by the geographic extent displayed in Fig. 4E, where we display predicted velocities against measured quantities. We show predicted velocity residuals in fig. S19.