Kinematics and paleogeology of the western United States and northern Mexico computed from geologic and paleomagnetic data: 0 to 48 Ma

Since Wilson (1965) defined plate tectonics, the approximation of oceanic lithosphere as a set of rigid caps whose relative motions are described by Euler poles has been very successful. Many researchers have attempted to apply similar rigid-microplate methods to continental tectonics. Certain first-order faults are selected as microplate boundaries, and the microplates between them are separately restored to their former positions through geologic time. This may be done manually by cutting a geologic map (e.g., Dickinson, 1996; Ingersoll and Rumelhart, 1999; Burnham, 2009) and working in the map-projection plane, or with computerized application of Eulerian finite rotations on a sphere (e.g., McQuarrie and Wernicke, 2005; Wilson et al., 2005).

One problem with this approach is that the choice of boundaries of microplates is difficult where fault traces do not connect. In western North America, most normal faults in the Basin and Range province do not connect to transform faults at their ends, although Stewart (1998) documents some exceptions. Also, the dextral splays of the Calaveras fault in northern California do not connect back to the main Pacific–North America (PA-NA) plate boundary at the coast (Jennings and Bryant, 2010). These facts indicate that continental microplates in our study area are not rigid. Another indication of deformation mode is that restored microplates in most models do not fit together, and authors must artificially separate them or allow them to overlap. Finally, a combined geologic and geophysical model of neotectonics in the western United States by Bird (2009) found that fully a third of PA-NA relative motion in California is taken up by permanent deformation within crust segmented by the named and modeled faults.

A second criticism of the rigid-microplate approach is that during each time step, the motion of each microplate has exactly three degrees of freedom (i.e., latitude and longitude of the Euler pole and its angle of rotation, or angular rate). Typically this means that certain favored constraints are enforced exactly, while other relevant data are not used. The selection of favored constraints is subjective, so the modeling process is not fully reproducible.

For these reasons, one of us attempted to devise an alternative modeling formalism (Bird, 1998), in which the crust is a deforming continuum, and faults are just sites of concentrated weakness (and sometimes measurable strain) within this continuum. This approach greatly increases the number of degrees of freedom in the model, often beyond the number of data. Therefore, stabilizing and regularizing assumptions are added: the internal deformation of microplates is minimized, and it is required to be consistent in orientation with interpolated stress directions. An advantage is that unlimited data, where available, may contribute to the solution with natural weighting by the inverse squares of their standard errors (as discussed below).

As we explore the potential importance of distributed (continuum) permanent deformation in southwestern North America, it would be logically inconsistent to pretend that it does not occur in other plates. Therefore, we regard the results of global plate circuits (e.g., from stable North America through Africa and Antarctica to the Pacific plate) as rough estimates that are likely to include systematic bias. We do not use them to provide boundary conditions on our solutions, nor do we plot them for comparison on our paleogeologic maps. The comparison and reconciliation of our models with global plate circuits is a large topic for future research, but that effort must consider the likelihood of systematic errors in both kinds of reconstruction.

Our new approach (and algorithm) for modeling have the ability to test published kinematic concepts about paleotectonics. They also make available the residual misfit associated with each datum; analyzing these by data class may reveal systematic conflicts between types of data. For example, this study addresses the frequent disagreement between latitude changes inferred from paleomagnetic data and latitude changes obtained by restoring fault slip; this conflict has been the basis of the extended “Baja British Columbia” controversy of the past 40 years (e.g., Cowan, 1994; Butler et al., 2001). Finally, such modeling has discovery potential, as unanticipated kinematic patterns may emerge. Here we apply this approach to compute the paleokinematics of southwestern North America back to 48 Ma, generating paleogeologic maps, and maps of velocity, strain rate, and fault heave rate at intervals of 6 m.y., based on calculations with timesteps of 0.2 m.y.. We discuss the implications of our findings for proposed correlations across the San Andreas plate boundary, for large rotations of the western Transverse Ranges of California, and for the initial response of continental North America to its first contact with the Pacific plate ca. 28 Ma. Further extensions of the model back to 90 Ma will be presented in a future paper.

We have compiled most of the kinematic information on the tectonics of southwestern North America since the mid Cretaceous from scientific journals, special publications, government sources, theses, and dissertations. Information was first recorded as notes with full citations in bibliographic software. Then values for displacements, uncertainties, and timing constraints were copied into summary tables with short-form citations.

Our survey extends from the 20°N parallel (Trans-Mexico volcanic belt) to the 49° N parallel (Canadian border). On the east, it extends to the 96°W meridian in the Great Plains, which was part of stable North America during this time. On the west, it extends to the outer continental edge, where water depths reach ~2 km (e.g., Patton Escarpment); during the time interval between 90 and 28 Ma, the oceanic lithosphere farther west was largely decoupled from North America by rapid slip on the subduction megathrust at the Franciscan trench (e.g., Hamilton, 1969; Atwater, 1970; Engebretson, 1983; Engebretson et al., 1985). Our summary of events since ~90 Ma includes the late Sevier orogeny, the Laramide-Hidalgoan orogeny, Basin and Range extension, and Cenozoic strike-slip tectonics (e.g., Dickinson, 2006).

One goal in presenting these data tables is to clarify the origin of our kinematic targets; another is to assist other researchers in more rapidly attaining an overview of the literature, so that they can expand these data sets and/or analyze the data using other methods. The description below covers our data sets of fault traces, fault offsets, restored sections, rotations and translations based on paleomagnetic data and stress directions, including timing constraints.

Every field geologist understands that a fault trace on a geologic map is a simplification of the complex quasi-fractile fracture sets seen in outcrops, and that the amount of simplification depends on the scale of the map. In a similar way, our data set of fault traces is a simplification of published traces, appropriate to map scales such as 1:1,000,000. Swarms of parallel or en echelon fractures or splays were simplified as one central trace. Intermittent cover by Quaternary alluvium or young volcanic rocks was ignored by including dotted traces which suggest probable fault continuity at seismogenic depths. Strike-slip faults of large offset were extrapolated under covering units to connect to adjacent faults. Early-generation extensional faults exposed in ranges in Nevada, but not in the intervening basins, were extended half-way across each basin to prevent systematic underestimation of the aggregate amount of early phases of Basin and Range extension. Detachment faults that are not exposed (or poorly exposed), but which have been inferred from metamorphic core complexes, were drawn with strike orthogonal to the primary lineation in the core complex, and with a length determined by the dimensions of the complex.

These editorial decisions were recorded as colored traces drawn by hand onto photocopies of published fault maps. Then these colored traces were assigned unique 4-digit serial numbers, and digitized. At very large scales, longitude and latitude coordinates on the published figure were treated as if they were orthogonal Cartesian coordinates while digitizing. At smaller scales, the first-author's utility program Projector (available at http://peterbird.name/oldFTP/Digitise/) was used to convert between coordinate systems. All traces are recorded as sequences of longitude and latitude with °E and °N considered positive, respectively. The file format is the .DIG format documented at http://peterbird.name/guide/dig_format.htm. Our convention for the dip of nonvertical faults is that they dip to the left when progressing along the trace in the direction that we digitized it.

To ensure uniform coverage of the Basin and Range province in the U.S., we included all normal faults mapped in the review papers by Stewart (1978, 1998) and dePolo (1998). To ensure uniform coverage of faults in northern Mexico, we included all faults from the Tectonic Map of North America (1:5,000,000; Muehlberger, 1992). Initially, we had no names for many of these traces, so we assigned generic names, such as “normal fault #S116, 41N, 115W, NV.”

Some of these traces have been previously documented (e.g., Bird, 1998, 2007, 2009). In 2013, the Community Fault Model project of the Southern California Earthquake Center released Fault Model 3.1 for a seismic-hazard study (Plesch et al., 2007; Field et al., 2013); these more accurate traces were used in place of our previous estimates. However, our division of long traces such as the San Andreas fault into “fault trains” (Bird, 2007) of approximately uniform offset history differs from the segmentation in FM 3.1.

One important omission was discovered after calculations were completed: Our data set does not include a major listric normal fault that roots east beneath the Baja California microplate, named the Santa Margarita–San Lazaro fault (Fletcher et al., 2003, 2007; Ferrari et al., 2013). This fault has the same strike length as the Tosco-Abreojos fault, which also has east-directed heave. This shows that the Magdalena platform has undergone transtensional strain and not just dextral strike slip.

More generally, it is inevitable that our compilation is missing other important faults. We hope that readers will find our data compilations useful as a starting point for expansions and extensions. Since fault networks are quasi-fractile and potentially unbounded in length, there is much to do. Yet seismology offers a hopeful argument: A large fraction of the total seismic moment release of the Earth is already captured in our catalogs of great and moderate earthquakes, so the addition of the many small ones not yet detected is unlikely to change global tectonic models, although they are locally important. Similarly, it seems likely that the faults already incorporated in our model were responsible for the great majority of the regional strain.

This fault-trace data set (Fig. 1) is file fUSM_20200901PB.dig and is included in the Supplemental Material¹. It contains 2385 traces with total length of 158,436 km. For the convenience of readers, a KML version of this file is also provided in the Supplemental Material; this can be used to display fault traces in Google Earth. Because fault names are not standardized, and also because many fault names were not known to us, searching for a certain name is not the best way to find a desired fault in this file. Instead, plotting software such as Google Earth or our program RetroMap4 (Supplemental Material) should be used to plot local views of the data set with identifying serial numbers on each trace. The fault-offset spreadsheet fUSM_PRIME, described below, will then give the fault name (or location), together with citations of literature containing further information.

A fault offset is one component of its slip, as defined with conventions and codes of Bird (2007): R and L are right-lateral and left-lateral heave components parallel to the fault trace; D (extensional) and P (compressive) are heave components perpendicular to the fault trace, and N (normal; extensional) and T (thrust; compressional) are fault throws. By use of these six codes, we avoid the inconvenience and possible confusion associated with negative offsets.

In the Colorado Plateau of the USA, and in the Sierra Madre Oriental of Mexico, there are many examples of monoclines that do not expose fault traces. Our approach is to interpret these as forced folds over reverse faults in the metamorphic basement below the sedimentary strata, and to measure a throw (T) as the structural height of the monocline (which is usually greater than its topographic height).

Because our computed models (described below) are models of horizontal velocities and displacements on the spherical surface of a planet, fault throws (N or T) must be converted to strike-perpendicular components of fault heave (D or P) prior to their use in modeling. This conversion is made using the actual fault dip where that is available as a notation (e.g., “dip_degrees 75”) in the fault-traces .DIG file. In other cases, it is made using standard default dips for normal and thrust faults, which are 65° and 25°, respectively. Note that this default dip would be a very poor description of most low-angle detachment faults, and that offsets of type D are strongly preferred to characterize them. Offsets of class D are also the best way to summarize extensional faults that are listric in geometry, or which involve domino-style extension with rotation of faults (e.g., Wernicke and Axen, 1988). Seafloor-spreading at ridge segments in the Gulf of California is also described using offsets of type D, even though dike intrusion may be the dominant process, and there may not be any clear master fault.

The fault-offset-file input to our modeling program can provide either throw or trace-perpendicular heave for a particular dip-slip fault trace, but not both. The data editors decide which to use according to information in the literature. For example, most high-angle normal faults in the Basin and Range province are characterized by their throw (offset type N), but most low-angle (and/or rotated) detachment faults are characterized by offset type D. Any fault with oblique slip will be characterized by one strike-slip offset (R or L) and also one dip-slip offset (N or T or D or P). These two components may have identical, overlapping, or distinct time windows. Any fault that is known to have had separate episodes of slip in its geologic history, or moved at a variable rate, may be characterized by any number of sequential offsets of the same offset type, as long as their time windows do not overlap.

The fault-offset-table input to our modeling software has the following columns: (1) serial number of the fault trace, immediately followed by a one-character offset type indicator; (2) fault name and location; (3) offset target in km; (4) estimated standard error of this offset in km; (5) beginning time (greater age) in Ma; (6) ending time (lesser age) in Ma. Where available, a neotectonic offset rate and its uncertainty may also be entered in columns (7) and (8); however, this information was used only in the “Neotectonic Calibration” phase of our study, described below.

Our effort to summarize offset distances (columns 3 and 4) and time windows (columns 5 and 6) was long and difficult. Many references provide only lower limits, upper limits, or rough estimates of offsets and/or time windows. Also, many published references disagree about these values. For this reason, we collected all published opinions in an expanded version of the fault-offset table, postponing final editorial decisions until the end of our literature search. In the extra rows of this expanded table, the first two columns are filled differently. The first column contains a reliability grade, according to these codes: A = original report of new geologic data from field or laboratory; B = synthesis of geologic data about one fault, basin or range; C = regional synthesis of geologic data in a state or larger region; D = opinion whose basis is unclear; and E = (our) editors' comment. The second column is a short-form citation, possibly followed by clarifying detail such as the original authors' name for the fault, or phase of movement in a complex history. It is this expanded table that we supplementally provide as fUSM_PRIME_20210122PB.xlsx (see ^{footnote 1}). This file has 16,891 rows, including 2841 rows that represent our editors' summaries, plus 14,050 rows with supporting data and short citations. Although it is in a spreadsheet file format, it contains no formulas and performs no calculations. The entries are sorted primarily by state or province, and secondarily by fault name. Within the data block for each fault, supporting citations are listed in order of publication.

Our editorial review of fault offsets was performed jointly by both authors between the literature-search and computational phases of this project. In choosing offset targets, we gave greater weight to A-class references than to B-class, and greater weight to B than to C. Within a class, we gave greater weight to the sources most recently published. If the uncertainty was stated in a recent A- or B-class reference, we tended to accept that. Otherwise, we considered the whole range of published views to characterize the uncertainty. It is important to note that the uncertainties shown are standard errors (σ), and therefore the ±σ range is only asserted to contain the actual fault offset with ~66% confidence. (To attain 95%-confidence theoretically requires a wider range of ±2σ, according to a Gaussian-distribution model.) In merging information on time windows, we gave greater weight to crosscutting relations involving radiometric dates than to those involving fossil-based ages of sedimentary strata.

There are many cases in which our literature search did not yield any time-window information for certain faults. To supply an estimate, we made temporary working maps which displayed the time windows published for nearby faults of the same offset type (e.g., N or T or R or L) and visually interpolated. If neighboring dated faults were missing, few or far away, we relied on class-C regional references such as de Cserna (1989), who asserted that all normal faulting in Sonora is Miocene or younger, and Oldow et al. (1989), who described it as occurring during 20–10 Ma.

Another difficult case that we encountered frequently was that of a fault for which we knew the type of offset, but found no published estimate of the offset distance. Fortunately, our modeling algorithm (described below) can work with very uncertain input data, as long as the uncertainty is quantified. For this purpose, we divided all the faults in our database with known offsets according to their offset type (R, L, N, T, D, and P) and graphed the distribution of these known offsets separately for each offset type. All results are shown in Table 1, and the detailed distribution of offsets of type N is shown in Figure 2. The median of the appropriate distribution was used as the estimate for any fault of unknown offset. However, we did not use the formal standard deviation of each distribution as our uncertainty. Instead, we chose a smaller σ that gave a Gaussian distribution approximating the actual distribution in its main body, ignoring its high-magnitude tail. (Our argument for this method is that the high-offset faults are preferentially named, mapped, and studied to determine their offset histories. Therefore, they are not representative of the “ordinary, anonymous” faults that lack offset information.) Note that for every offset type in Table 1, the assumed standard deviation is at least 50% of the assumed offset; this permits an offset of zero within the Gaussian 95%-confidence bounds.

Another important point about the magnitudes of fault offsets is that almost all offset distances quoted in the geologic literature (and collected into our tables) are maxima along that fault trace. Well-connected throughgoing faults that make up a primary plate boundary may each have uniform offset along their traces, so that there is no important distinction between the maximum and the mean. However, isolated faults are often observed to have offset tapering to zero at each end. If the distribution of offset versus length were elliptical, then the ratio of mean offset to maximum offset would be π/4 = 0.79. It will be useful to recall this overestimation bias when we discuss the fitting (or under-fitting) of fault offset data by our models in a later section.

Several special explanations are needed regarding the system of 11 spreading centers and 12 dextral transform faults in the Gulf of California/Mar de Cortés (Muehlberger, 1992; Martín-Barajas et al., 2013). Unlike most continental faults that we treat as embedded in a deforming continuum, seafloor-spreading systems move rapidly with respect to both surrounding plates, leaving abandoned extensional centers amid new crust, and converting active transform faults to inactive fracture zones as they move away from the ridge. To handle this, we have programmed a different kinematic treatment for the faults that make up symmetric spreading systems; this new algorithm is described in the  Appendix. Another kind of time dependence probably occurred as an initial irregular and disorganized rift in the former Peninsular Ranges magmatic arc evolved into the present stable spreading geometry (Umhoefer et al., 2018). We are unable to model this process because too much remains unknown about the locations of former active faults and their individual slip histories; therefore, we merely project the present ridge/transform geometry back in time without any topological changes. For example, we do not attempt to model the likely jumps in spreading location that occurred near Isla Ángel de la Guarda (Nagy and Stock, 2000). Finally, we do not have enough offset data (pairs of piercing points) to define the slip history of each Gulf fault individually. Instead, we note that most of the Baja California peninsula has been essentially rigid during the Cenozoic, and take advantage of this circumstance by merging together all available constraints on slip history, regardless of latitude (Table 2), and then applying this history to each of the 23 master faults in the Gulf. The result is a crude cartoon of the actual tectonic history of the Gulf, but one that provides the proper kinematic boundary condition on tectonics of the California Borderland and Transverse Ranges.

Figure 3 presents a ribbon-plot of the net fault heave targets in the period 48–0 Ma. Note that the width of each ribbon is plotted as constant along the trace; however, actual fault heaves may have been locally less than shown, and more variable.

In regions of primarily dip-slip faulting and folding, many structural geologists have published cross sections, some of which are accompanied by their interpretation of a restored section with a different original length. We capture this information on the elongation/contraction of sections as an estimated constraint on one component of the relative displacement between the two end points of the section over some window of geologic time. If there was also a component of relative motion between end points that was perpendicular to the section plane (whether due to strike slip, or to rotation), then to first order, this has no bearing on the validity of the datum. Our collection of restored sections (Fig. 4) and their estimated length changes are available as file cUSM_20200701PB.xlsx in the Supplemental Material (see ^{footnote 1}). It has 90 entries, mostly concentrated in the Rocky Mountain foreland of the United States, the Sierra Madre Oriental of Mexico, the Basin and Range province around southern Nevada, and the Transverse Ranges of California.

Data obtained from the orientations of paleomagnetism include poles, inferred from a sample set that averages secular variation, and virtual geomagnetic poles (VGPs), whose averaging of secular variation is uncertain. These can constrain kinematic reconstructions in two ways. Because the apparent polar-wander path for stable North America is known (e.g., Torsvik et al., 2012), declinations to poles and/or VGPs from deformed western regions yield declination anomalies, which are interpreted as net rotations about vertical axes since the time of magnetization acquisition. Inclinations can be compared to a dipole-field model of expected inclinations to yield paleolatitude anomalies, thus providing an estimate of the N-S component of crustal displacements. Our paleomagnetic data set was assembled from three sources: relatively uniform regional coverage from the International Association of Geomagnetism and Aeronomy (IAGA) Global Paleomagnetic Database (GPD, http://www.ngu.no/geodynamics/gpmdb/, accessed 2017.04.10), plus focused coverage of southern California in a UCSB data set contributed by Bruce Luyendyk, and some recent data from Baja California and Sonora.

Several decades of debate about the tectonic significance of different kinds of paleomagnetic data suggest the need for certain quality controls. Data from plutonic rocks are suspect because their lack of reference to the paleohorizontal at the time of remanence acquisition makes tilt-correction difficult (Butler et al., 1989, 1990; Ortega-Rivera et al., 1997; Dickinson and Butler, 1998). Also, many sedimentary rock types are susceptible to inclination shallowing, due to compaction and diagenetic processes and/or remagnetization (Dickinson and Butler, 1998; Tan and Kodama, 1998; Butler et al., 2001; Hillhouse, 2010). Volcanic rocks are less affected. Lacking the expertise to deal with these problems, we limited our data query in the IAGA-GPD to primary magnetizations from volcanic rock sequences with ages less than 100 Ma; this yielded 189 poles and/or VGPs. Four of these had implausible declination anomalies for their locations (outside California) and ages, and we rejected them as probably suffering from inadequate averaging of secular variation; in three of these cases the original authors had stated similar reservations. The UCSB data set adds another 26 poles in southern California; while some of those are from sedimentary strata, all have been screened for tectonic relevance. Finally, we added 44 Baja California and Sonora volcanic poles and/or VGPs from Bennett and Oskin (2014). Our working collection of 255 poles and/or VGPs (Figs. 5 and 6) is available as file pUSM_20200701PB.xlsx in the Supplemental Material (^{footnote 1}).

Our algorithm uses interpolated stress directions to constrain the orientations of strain-rate tensors that are allowed in unfaulted continuum between modeled faults. Geologic indicators of paleo-stress direction include dikes, veins, fault swarms, and stylolites. Although technically each of these data types provides the orientation of paleostrain, it is conventional to call them “stress directions” (e.g., Bird, 2002) and we continue that usage. The convention is also to report stress direction as the (present) azimuth of the most compressive horizontal principal stress that was active when the feature formed. From the literature about the western United States and northern Mexico, Bird (2002) collected 369 such indicators and documented three regional stress-direction changes of high significance since 85 Ma. For this project, we added 66 indicators, primarily from Varga (1993), Bellier and Zoback (1995), Erslev (2001), and Ferrari et al. (2002). This data set (Fig. 7) is provided as file sUSM_20200701PB.xlsx in the Supplemental Material (^{footnote 1}).

Because these paleo-stress-direction data are sparse in many regions for many epochs, we chose to activate the option in our program Restore that treats fault traces as stress-direction indicators in their first phase of movement. Thus, most-compressive horizontal principal stress is assumed parallel to new normal faults, orthogonal to new thrust faults, and at a 45° angle from traces of new strike-slip faults. (The angle of 45° is arbitrary but represents a compromise between angles of 25°~30° typically observed in laboratory rock mechanics, and angles up to ~75° sometimes seen in the field.) Faults that are reactivated for a second phase of slip are considered less reliable indicators because of the preexisting plane of weakness, and not used.

Our algorithm uses Bayesian prior constraints (or “pseudo-data”) stating that the scalar strain rate in all unfaulted continuum between modeled faults is zero. The uncertainties (μ) attached to these goals of zero are a kind of input data because they control how strictly the prior constraint is enforced in each region, and roughly how much continuum strain rate might appear in the solutions. These μ values are in units of scalar strain rate (in/s) and could also be thought of as “continuum compliances.” We use these values to express the rates of deformation that are likely to be missing from our compilation of fault traces and fault offsets, for whatever reason.

The regional or baseline μ in this study is 5 × 10⁻¹⁶/s, for reasons explained below. Other local values are listed in Table 3. A higher value is appropriate in the Franciscan Complex (and associated exotic ultramafic rocks) of California, because these rocks were scraped off subducting oceanic plates and accreted to North America (e.g., Cloos, 1982); the faults and strains associated with accretion make it difficult to detect younger deformation. To quantify this, we note that along all the many northward splays of the Neogene–Quaternary Calaveras dextral fault (which splays lie mostly in the Franciscan Complex), the total slip rate declines from 160 km since 12 Ma at the south end to zero to the north, in a distance of 250 km. Partitioning the implied strike-parallel strain rate on both sides of the system suggests that strain rates have reached 8.5 × 10⁻¹⁶/s without creating structures that compilers (e.g., Jennings and Bryant, 2010) can confidently identify as active in the Quaternary. We also apply a very large μ, one order of magnitude above the baseline, to the terranes now covered by the flood basalts of the Columbia Plateau and Snake River Plain, for times prior to these eruptions, because the older geology is entirely hidden and may contain important faults that were active in pre-eruption times. Similarly, we discovered during the restoration process that important amounts of extension that occurred between the mapped detachment faults in the Colorado River extensional corridor have probably been hidden by Miocene volcanic rocks and younger sedimentary sequences of the Colorado River. The choice of μ for these regions is discussed below.

These μ values enter our calculations when they are assigned to individual elements of our finite-element grid(s). We chose the appropriate elements by using mapping software to overlay these grids on the digital Geologic Map of North America (Garrity and Soller, 2009). Our initial F-E grid, NI01.feg, is in the Supplemental Material (see ^{footnote 1}).

A broad conceptual overview of the algorithm of program Restore includes these steps:

1. A 2-D grid of contiguous but non-overlapping spherical-triangle finite elements is drawn to define the model domain, which must be simply connected in the topological sense.

2. At the start (younger end) of each restoration time step, which has duration 0.2 m.y., the relevant data and pseudo-data are assembled in a list with subscript k. This list includes azimuths of the most-compressive horizontal principal axes of paleostress interpolated to the centers of all finite elements, and μ values for each.

3. All target displacements or rotation angles d_k are divided by the length of their associated time windows to yield rate targets: . In the same way, their associated standard errors Δd_k are converted to standard errors in the rate: .

4. The objective function is written as a sum of quadratic terms involving model predictions p_k compared to target rates r_k, normalized by uncertainties. Each term is a multiple of [(p_k – r_k)/σ_k]². The factor applied to each squared prediction error (p_k – r_k)² could be called the “natural weight” for that datum, because it appears algebraically when the weighted-least-squares objective function is derived from the maximum-likelihood criterion, under the usual assumption of independent data with Gaussian uncertainties. This natural weight means that each datum with a large uncertainty has less importance in the objective function, and that our solutions will be primarily guided by those data that have smaller uncertainties.

5. Quadratic terms involving fault offset-rate targets are multiplied by an additional weighting factor of (L_k/L₀), where L_k is the length of the fault trace in one finite element, and L₀ is a model control parameter with units of length. Also, quadratic terms involving continuum strain-rate targets and their principal-axis azimuths are multiplied by an additional weighting factor of (A_k/A₀), where A_k is the area of one finite element, and A₀ is a model control parameter with units of area. Point data at paleomagnetic study localities and ends of cross sections have multipliers of unity. These extra weighting terms could be called “element-discounting weights” because their purpose is to make the objective function (and the eventual reconstruction solution) independent of the number and size of finite elements used in the grid. In turn, this gives each modeler the opportunity and responsibility to adjust L₀ and A₀ to yield a model that balances the fit to different classes of data (i.e., point data, line data, and area data).

6. Complex algebra is used to replace all model predictions p_k with weighted sums of the horizontal velocity components of adjacent nodes of the F-E grid.

7. Differential calculus is used to find multiple independent constraint equations describing the one point in nodal-velocity space that will yield the minimum of the objective function.

8. This system of simultaneous linear equations is solved by parallel matrix processing to obtain the horizontal components of velocity at each node. Each velocity solution is iterated six times to let two nonlinear terms in the objective function converge. (These terms are related to the loose enforcement of interpolated-stress-direction constraints on the orientations of the principal axes of continuum strain rates, as detailed in section A.7, “Use of Stress Directions,” in the technical appendix of Bird, 1998.)

9. Positions of all nodes are projected back to the earlier end of the time step. All model data are also translated as points embedded in the F-E grid.

10. Using the new node positions, all equations are re-formed and re-solved, yielding slightly different velocities for each node.

11. The position of each node is corrected by: the time step multiplied by one-half of the change in vector nodal velocity from step 8 to step 10. This is the predictor-corrector method of time-integration. The paleo-positions of all model data are also corrected.

12. If the model's target geologic age has been reached, or if any finite element is excessively deformed (so that it would fold over and have a negative area), the program outputs several files and then stops. Otherwise, it loops back to step #2 above.

All equations (especially the algebra of step #6) were published in appendix A: Algorithm in Bird (1998). That reference also describes an option to iterate the entire time-history of a reconstruction many times, in order to adjust the allocation of target rates across time steps; that experimental feature was not used in this project. Some new algorithmic steps regarding the translation of plate-boundary fault traces through time, which are new to version 4 of program Restore, are described in the  Appendix to this paper. The Fortran 90 source code and an executable version of Restore v.4 are provided as Supplemental Material (see ^{footnote 1}), Restore4.f90.txt and Restore4-Win64par.exe.

Values of the continuum compliance parameter μ are listed in Table 3. The data-class weighting parameters L₀ and A₀ also have subtle and systematic effects on the solution. Here we explain how preferred values for μ, L₀ and A₀ were chosen.

During the first time step of any restoration project, starting at 0 Ma, the calculation initially yields neotectonic velocities. These are very comparable to those obtained with related finite-element program NeoKinema, which was used (e.g., Bird, 2009; Field et al., 2013) in systematic studies of the western United States for seismic-hazard models. The primary difference is that NeoKinema tries to fit thousands of highly accurate GPS velocities at benchmarks and does not use any paleomagnetic data. Bird (2009) reported extensive testing of the effects of continuum compliance μ, and found that a value of 5 × 10⁻¹⁶/s is optimal for neotectonics in the western United States, both because it permits fitting all input data sets adequately, and because it is a self-consistent value yielding model output that reports the RMS rate of continuum deformation as almost the same as the input value of μ. Therefore, we adopted this value as the default value for most of the continuum when modeling paleotectonics.

The neotectonic study of Bird (2009) also suggested preferred values of L₀ and A₀. However, before adopting these, we performed a set of 70 neotectonic calculations with Restore, trying all combinations of ten values of L₀ with seven values of A₀. All tested values differed from neighboring values by a factor of two. Various measures of model quality were then plotted as a function of the logarithms of these two control parameters. We investigated errors in fault offset rates, both in relation to their long-term geologic target values and in relation to independent Holocene estimates. We also investigated mean errors in the fitting of principal axes of continuum strain rate to interpolated stress directions. Finally, we investigated various measures of the mismatch between neotectonic velocities from Restore and 6679 interseismic velocities at GPS benchmarks (not used in the Restore calculation). The most useful diagnostic turned out to be the mean ratio of Restore velocity to GPS velocity, averaged across all GPS benchmarks; optimal values of this ratio (close to unity) were found along a particular contour in (log₁₀L₀, log₁₀A₀) space. From these diagrams, we selected the preferred values of L₀ = 2 × 10⁴ m and A₀ = 1 × 10⁹ m² that were used for most of the modeling in this study. The identifier for this reference model, in filenames and in graphics is “NI.”

However, we were concerned that perhaps this reference model does not give adequate weight to paleomagnetic data, because we noted that it did not fit their declination anomalies in detail, but only in general trend. Therefore, we computed two additional restorations to test lower weights on fault and continuum constraints, implying greater relative weight on point data.

Model “HP” was run with L₀ = 3.2 × 10⁵ m and A₀ = 1.6 × 10¹⁰ m². Because each of these parameters is 16× higher than in model NI, the relative weight on all paleomagnetic data (and cross-section elongations) was effectively increased by 16×. However, this model resulted in generally northwestward velocities across coastal North America that were too high (relative to GPS) by ~20%. This defect did not change across the time span of 0–6 Ma, after which this calculation was abandoned. The problem is that merely increasing weight on paleomagnetic data excessively increases the northward “pull” exerted by paleomagnetic paleolatitude anomalies. This “pull” is valuable in our kinematic reconstruction models because it plays the same role as PA-NA plate coupling (Bird, 2002) would in dynamic models based on rheologies and quantitative stress fields; however, the strength of this surrogate effect must be calibrated.

Therefore, we also computed an alternate model “PM” with L₀ = 3.2 × 10⁵ m and A₀ = 2 × 10⁹ m². In this model, the more modest change in A₀ (× 2) had the desired effect of regulating the overall northwestward drift in the model, while the high value of L₀ meant that weight on paleomagnetic data was increased by 16× relative to fault-offset data. This model was run over the timespan of 0–30 Ma. Unfortunately, the unforeseen defect of this model was that continuum stiffness constraints had now been increased by 8× relative to fault-offset data (all relative to reference model NI). This resulted in generally under-fitting all fault offsets, more than before. For example, detachment faults in the Colorado River extensional corridor achieved only ~50% of their target slips, as opposed to ~80% in model NI. Also, the coastal Hosgri fault was under-performing more seriously: Compared to a target of R = 156 ± 10 km in 11–0 Ma, model PM calculated only R = 76 km, compared to R = 114 km in reference model NI. Therefore, we present only maps and graphs of results from reference model NI in this paper.

Restoration proceeds by alternating between machine processing and manual oversight. With our particular data sets and F-E grid, automated processing requires ~14 h of computer time (on a 64-bit, 3.3 GHz computer with 64 GB memory that uses 16 parallel processors when solving linear systems) per 6 m.y. of geologic time. However, the program rarely progressed further than ~0.8 m.y. without an automatic stop for manual correction and renumbering of a distorted F-E grid. Also, each time another 6 m.y. had been modeled, we manually prepare a standard set of maps and graphs to examine output. Auxiliary programs used to prepare input and map output are available in the Supplemental Material (see ^{footnote 1}) and/or at http://peterbird.name. Experience shows that our typical rate of progress is ~16 m.y. of reconstruction per month.

The initial F-E grid was created in program OrbWin as a mesh of equilateral 30-km-wide spherical triangles obtained by the 8^th iterated subdivision of a global icosahedron. Then, program Fault_Corridors was used to modify this F-E grid by inserting 4-km-wide “fault corridors” of thin elements along most fault traces. Wider corridors, up to a maximum width of 50 km, were used along major extensional faults to allow for the anticipated losses of crustal area during restoration. Further manual editing in OrbWin was needed to stitch these F-E grid components together and fix all topological problems.

After any automated stop, the user must survey the deformed F-E grid in OrbWin and edit it to fix topological problems and anticipate further retrodeformation strain. Along major strike-slip faults, deformed elements within the fault corridors are deleted and redefined, combining existing nodes in new triplets. At spreading centers, it may be necessary to delete nodes that have been restored too close together, and then regrid the area with a smaller number of new elements. After each manual regridding, the F-E grid is renumbered by utility program OrbNumber to reduce the bandwidth of the linear systems in Restore; boundary conditions must then be re-stated in terms of these new node numbers. Our boundary conditions are simple: the eastern edge of the model is always fixed with respect to stable North America, and all other model edges are free.

After every 6 m.y. of restoration, maps were produced (as Adobe Illustrator .AI files) by utility program RetroMap4. Available plots of output include fault heave rates (or net heaves) shown with colored ribbons, continuum strain or strain rate shown as swarms of equivalent conjugate faults, principal stress directions, colored contour maps of velocity and rotation-rate, and paleogeologic maps (described below). Inspection of these plots can sometimes catch errors; for example, strike-slip on an isolated fault trace might result from clerical error in recording its time-window of activity. Or, a belt of concentrated continuum strain linking two modeled faults might indicate that some intervening fault has been inadvertently omitted. In addition, the quality of fits to paleomagnetic constraints, fault offsets, and cross-section elongations was assessed with scatter plots of model versus data values.

All codes and data sets needed to reproduce our results are included in the Supplemental Material (see ^{footnote 1}). However, a significant weakness of our modeling method is the need for frequent interruptions of the calculation(s) for manual re-gridding of fault corridors. This causes the time required for most restoration experiments to stretch into weeks or months, and in turn, this discourages more extensive experimentation with different data-weighting or data-editing options. Therefore, such modeling would benefit greatly from the addition of an automated algorithm for intelligent local regridding (anticipating the new strains that the next time steps are likely to bring) that did not require calculations to stop. This would encourage more experiments.

One such experiment that could improve our results would include working with different sets of paleomagnetic data, such as a subset of only poles that include full averaging of secular variation, and/or a superset that made use of declination anomalies from sedimentary rocks while omitting their inclination anomalies.

During computation of reference model NI (0–48 Ma), there were four times and/or places when and/or where the authors needed to intervene by making changes to the input data to achieve better realism. A map of all features named in this discussion is shown as Figure 8.

In the area north of San Gorgonio Pass in southern California, the San Andreas fault system presently makes a left step of 18 km width between the Mission Creek segment and the San Bernardino segment. This has apparently increased resistance to dextral slip since ca. 3 Ma, causing the San Gorgonio Pass–Garnet Hill thrust to appear and the Banning fault to reactivate as dextral, taking up part of the relative plate motion. Our initial assumption was that this left step was caused entirely by ~16 km of left slip on the Pinto Mountain fault since 6 Ma (Langenheim and Powell, 2009), and that restoration of this slip would restore the San Andreas system to a smooth trace at 6 Ma. However, our expectation was not met, for two reasons: (1) model Pinto Mountain fault slip was a bit less, at 14 km overall, and less than that at its western end; and (2) Pinto Mountain slip was too slow to damp out the tightening-S-kink instability that appeared in our San Andreas trace during restoration. We were able to solve the problem by a combination of automated slip-rate-dependent smoothing of strike-slip fault traces (see  Appendix) and also by imposing thrust heave of P = 12 ± 2 km on the enigmatic Dillon fault (Hope, 1969; Rymer, 2000) in the Little San Bernardino Mountains foothills during 6.0–0.3 Ma. Although this problem may be an artifact due to deficiencies in our restoration process, the possibility of thrusting on the Dillon fault during development of the left step seems to warrant further consideration. This thrusting may also explain elevation of the Little San Bernardino Mountains.

Several authors (e.g., Crowell, 1975, 2003; Ingersoll et al., 2014) have advocated that San Andreas slip should be restored so that the San Gabriel Mountains were opposite the Chocolate and Orocopia Mountains ca. 6 Ma. (For more details, see Discussion.) Our first attempts at this restoration fell short by a consistent large deficit, ~40 km. In hindsight, this can be attributed to three causes: (1) Our weighted-least-squares method is likely to underfit the largest slip targets on any fault system, and the problem is especially serious on the southern San Andreas fault because of the left-step problem discussed above. (2) We probably overestimated the contributions of branching dextral faults (San Jacinto, Banning and Eastern California shear zone) when we subtracted them from the offset of piercing points to get the target slip for the San Andreas. (3) Uncertainties such as (2) above probably raised the uncertainties of slip estimates by previous authors as well (Darin and Dorsey, 2013), whereas the uncertainty in locations of piercing points is only a few kilometers. For the final model, we increased slip targets and reduced standard errors for segments of the San Andreas system; for example, the Mojave South target was adjusted from R = 166 ± 10 km to R = 200 ± 4 km, and the Mission Creek target was adjusted from R = 203 ± 10 km to R = 225 ± 4 km; the intervening San Bernardino North and South Branch targets got similar adjustments. After these changes, a reasonable restoration was achieved at 6 Ma (Fig. 9).

The Oceanside detachment fault is an east-dipping extensional structure that was active during 18–6 Ma; its trace in southern California is roughly co-located with the younger Inglewood-Newport–Rose Canyon system, but also extends farther into offshore Baja California. Santa Catalina Island is a well-known exposure of its footwall metamorphic complex. We initially relied on published estimates that its heave was D ≈100 km (Crouch and Suppe, 1993; De Hoogh et al., 2019). However, once we arrived at our 6 Ma reconstruction (Fig. 9), at the end time for Oceanside detachment slip, we reconsidered how much slip would be allowed by geologic constraints. Using the seafloor geologic map of Vedder et al. (1974), we noted that San Nicolas Island (composed of undeformed Eocene turbidites) is the easternmost unambiguous outcrop of non-exhumed geology representing the older stable Outer Borderland terrane. The distance from the eastern tip of San Nicolas Island to the Oceanside trace at 6 Ma was 140 km, and we decided to use that upper limit as the revised heave for the Oceanside detachment: D = 140 ± 10 km, at all latitudes south of its intersection with the Cristianitos fault.

We use the term “Colorado River extensional corridor” to refer to a group of 13 detachment faults (and one tear fault) extending from the Mojave Desert into southern Arizona, which accommodated major SW-NE crustal extension during ca. 25–12 Ma. During data collection, we noted that this fault system is not continuous; there are six gaps in which no detachment faults (or tear faults) have been mapped. As expected, this lack of continuity (and the requirement that our algorithm must minimize deformation in unfaulted continuum) made it difficult to restore this system. Initial results were that extension rates in the 21.0–20.8 Ma timestep were only ~17% of their targets in the mapped detachment fault systems. Our first remedy was to insert high-μ areas (Table 3) in the gaps between mapped metamorphic core complexes, based on the argument that linking structures probably exist, but they have been buried, either by Miocene volcanic and sedimentary rocks, or by younger Colorado River sediments. Also, the fact that younger tectonics (12–0 Ma) involved qualitatively different dextral strike-slip strains may have made first-generation extension harder to recognize. However, the addition of high-μ areas only brought extensional heave rates up to ~35% of their targets. After several experiments, we reduced the standard errors on the target detachment heaves from an average of ±10 km to ±2.5 km. Then, extensional heave rates reached 78% of their targets. This experience highlights one deficiency of our algorithm: because it assumes a Gaussian distribution for the PDF of offset on every fault, it does not provide any convenient way to specify lower limits on their offsets; instead, the operator may have to increase the effective weights on some offsets by artificially shrinking their uncertainties.

In the neotectonic velocity test, reference model NI had mean error of 2.5 mm/a and RMS error of 3.7 mm/a with respect to 6679 interseismic velocity vectors measured at GPS benchmarks. Some of this discrepancy is explained by the fact that the Restore model predicts long-term average velocities rather than the interseismic velocities measured by GPS; the distinction between these two velocities can rise to one-half of fault heave-rate for GPS benchmarks very close to an active fault. The mean azimuth error for principal axes of continuum-strain rate (with respect to interpolated stress directions) in this neotectonic model was 7.5°, comparable to the nominal minimum 10° standard error on the interpolated stress direction. The mean error with respect to target offset rates on all faults (also at 0 Ma) was 0.4 mm/a.

In reconstructing the study area back to 48 Ma, we used 240 time steps of 0.2 m.y.; each time step included two velocity solutions, each of which was iterated six times to allow convergence of nonlinear terms in the objective function (detailed in the appendix to Bird, 1998). Convergence was uniformly good; the mean of fractional velocity changes in the final iteration (RMS change/RMS velocity) was 1.5 × 10⁻³. The mean discrepancy between principal axes of continuum deformation and interpolated values was 6.3° at the end of these iterations.

Considering the 2000 fault offsets that occurred entirely in the age range of 48–0 Ma, and which had positive (nonzero) offset targets, 93.8% were less than fully achieved (“underfit”). In further detail, 54.5% were underfit by more than 20% of their targets, and 17.1% were underfit by more than 50% (Fig. 10). Or, computing prediction errors in nondimensional terms (as numbers of standard errors on each target value, shown in Fig. 11), the mean absolute value of these normalized fault offset errors was 0.899, and their RMS size was 1.631, both of which can be compared to the ideal of unity.

The systematic tendency of our weighted-least-squares method to underfit fault offsets (Fig. 11) is a natural result of the algorithm being simultaneously designed to minimize errors in continuum strain rates, and then being given target continuum strain rates which are all zeros, rather than the true values. (Our previous discussion of parameter μ only concerned changing the uncertainties on some continuum strain rates, not changing their targets of zero.) At present there seems to be no easy, widely applicable method to measure these ancient continuum strains in the field.

However, when assessing the seriousness of the general underfitting of fault offsets, it is important to remember that most offset target values from the geologic literature are maxima along the fault trace, but our algorithm and our posterior analysis compare these to the mean offsets along each trace, as estimated by our reconstruction model. We noted in Data Collection that this overestimation of offsets in our input data set probably ranges up to 21%, with a mean likely to be ~10%. Accordingly, one could argue that about half of our model fault offsets are roughly consistent with geologic reality, while half are significantly underfit.

It has also been suggested that this bias could be eliminated by reporting maximum model offsets (along each trace) instead of mean offsets. We decline this suggestion because each maximum offset would come from a single finite-element, and therefore would be noisy and irreproducible (with slightly different F-E grids). This issue of the noisiness of single-element estimates of fault offset was illustrated in figure 13 of Bird (2009).

Elongations and/or contractions along restored cross sections were also underfit. Among the 33 sections whose initial ages are <48 Ma, a regression analysis of model versus target elongations (Fig. 12) gave a slope of only 58% average achievement, but a good correlation coefficient of 0.962. All model elongations and/or contractions had the correct sign. One likely cause of this underfitting is that some of the faults shown on the published cross sections are not included in our map-view data set of fault traces.

To characterize the success of the enforcement of the minimized-continuum-strain rate prior assumption (relative to its spatially varying uncertainty μ), we can use summary statistics. For example, during the run of the reference NI model from 0 Ma back to 6 Ma, our three dimensionless misfit measures were: N₀ = 0.015, N₁ = 0.282, and N₂ = 0.606. Here, N₀ is the fraction of elements (weighted by their areas) where the mean (over timesteps) ratio of scalar strain rate to uncertainty μ exceeded 2; its preferred value is 0.050 or less. N₁ is the mean (over area-weighted elements, and over timesteps) of the ratio of (positive) scalar strain rate to μ; its preferred value is 1.0 or less. N₂ is the RMS average (over area-weighted elements, and over timesteps) of the ratio of scalar strain rate to μ; its preferred value is also 1.0 or less. However, these low relative misfits on an area-weighted basis can hide larger continuum strain rates that occur around fault bends and intersections. Thus, we also monitor continuum strain rates in map-view, plotting each strain rate tensor as an equivalent pair of conjugate faults, as in Figure 13. We can also integrate model continuum strains over several million years of reconstruction, as in Figure 14. Wherever rocks of appropriate age outcrop, there are opportunities to test these predictions.

There are 162 paleomagnetic poles and/or VGPs in our database with magnetization ages less than 48 Ma. Among these, a regression analysis of model versus target rotations gave a slope of 72.8% average achievement and a correlation coefficient of 0.894 (Fig. 15). (In the Neotectonic Calibration section, we discuss two attempts to reduce the errors in model rotations, and the reasons why we consider those attempts unsuccessful.) In all cases where the target rotation had magnitude greater than 14°, the model gave a rotation of the correct sense.

The one class of data that our model failed to fit (even approximately) was paleolatitude anomalies (or latitude changes) derived from inclinations of paleomagnetism. The correlation coefficient between target and model latitude changes is only 0.188, and the regression slope shows only 3% of target latitude changes achieved. Also, the ranges are very different. In our model, most inland localities had latitude changes less than 1°, and latitude changes in the coastal terranes were no more than 3.8° (all northward). In contrast, the target latitude changes based on paleomagnetism range from 30° northward to 12° southward, and fully 17% of the targets have amplitudes of greater than ±5°. It appears that paleomagnetic estimates of latitude change in the western United States are basically uncorrelated with latitude changes based on restoring fault offsets. This situation is similar to that in western Canada, where the “Baja British Columbia” controversy (e.g., Cowan, 1994; Butler et al., 2001) arose from disagreements between the same two classes of data.

We created two series of paleogeologic maps, extending from 6 Ma to 48 Ma in steps of 6 m.y. Both are based on retrodeformation of the Geologic Map of North America by Reed et al. (2005) and Garrity and Soller (2009). Our two series of maps are almost identical except in scale and coverage. One series displays the whole model area at scale 1:4,800,000 with polyconic projection centered on the 110°W meridian. At this scale, it is impractical to display the unit symbol on each outcrop area; however, these can be identified by comparing outcrop shapes and colors to the published maps. The other series presents the greater southern California and Baja California area (which is the area of greatest strain) in detail at scale 1:1,875,000 with polyconic projection centered on the 116°W meridian; this series includes unit symbols on almost every outcrop. Bitmap-image versions of these two series are attached to this paper as Supplemental Material (see ^{footnote 1}). Object-oriented paleogeologic maps in Adobe Illustrator .AI format are also included in the Supplement Material.

To prepare these maps, we first converted the outcrop, contact, and fault layers of Garrity and Soller (2009) from ESRI ShapeFile (.SHP) format to our .DIG format. Then, we restored most digitized point locations to their former positions with program Restore. However, points that fell within the same finite element as any plate-boundary fault with the “throughgoing master fault” or “symmetric spreading system” designation (see  Appendix) were not restored and not displayed, because their polylines would have been hopelessly distorted and misleading. Restored positions of all other polylines were plotted with our graphical program RetroMap4. In addition, we overlaid the restored traces of all faults in our database back to their estimated initiation ages, using colored lines to indicate rake and dip.

The International Stratigraphic Chart (Cohen et al., 2013) was used to define the age-range of the unit symbol associated with each outcrop polyline. Then, colored outcrop areas that had not yet formed were omitted from each paleogeologic map. However, the fine contact lines (black, 0.28-point) separating these areas were never deleted, for two reasons. First, the contact lines published by Garrity and Soller (2009) have no metadata, making such deletion difficult to implement. Second, we feel that the contact lines provide useful location references when comparing paleogeologic maps to the present geologic map.

Similarly, the fault traces in heavy (black, 1-point) lines from Garrity and Soller (2009) were never deleted. First, they also have no metadata. Second, defining the age of initiation of any fault is difficult, unless it forms in an area of continuous deposition and then is fortuitously excavated by later geologic events. These fault traces in black are often redundant with our modeled fault traces plotted in color (Fig. 16), although they do not align very precisely. This discrepancy was also noted in a 0 Ma test map and therefore does not imply imprecision in our reconstructions. In principle, the black fault traces should be useful to display the locations of fracture zones in the Gulf of California/Mar de Cortés because our software does not display abandoned fracture zones shed from a symmetric-spreading system.

We note several known deficiencies in our paleogeologic maps: (1) There are many areas of white, where the paleogeology is unknown due to younger covering units. (2) Restored outcrops based on present outcrops may actually have been buried by covering strata in the past. (3) Black fault traces from the Geologic Map of North America appear on older maps, before the time those faults actually formed. (4) In detail, it can be noted that truncations and offsets of outcrop areas at minor fault traces remain the same as they are on the 0 Ma geologic map; this is because the deformation that we restored along each minor fault trace was distributed over a narrow (e.g., 4 km) fault corridor, not localized in a discontinuity at the fault trace.

The geographic reference frame for all parallels and meridians on paleogeologic maps is stable eastern North America. Because of North American apparent polar wander, the latitudes shown are not directly comparable with paleolatitudes that might be obtained from paleomagnetic inclination data.

As explained in the Introduction, our restoration method is distinctive in two ways: (1) It does not assume rigidity of microplates bounded by modeled faults. (2) It uses nearly all available data within southwestern North America (with natural weighting according to their uncertainties) rather than being dominated by a few favored constraints or by external boundary conditions. Naturally, our results are quite different from previous models.

Figure 17 shows the southern California part of our reconstruction at 26.6 Ma and compares it to Figure 7A of Wilson et al. (2005). There are many important differences: (1) Wilson et al. (2005) restore both the Sierra Nevada and future Channel Islands ~190 km farther southeast (relative to stable North America). (2) They restore the Channel Islands ~80 km farther northeast than we do, relative to the Sierra Nevada and San Andreas fault. (3) They reconstruct Santa Catalina Island in the Borderland, whereas we reconstruct it in the subsurface below the hanging wall of the Oceanside detachment fault. (4) They back-rotate the San Gabriel Mountains block by ~45° more than we do. Thus, there is no simple Eulerian rotation that could reconcile their 26.5 Ma reconstruction with ours. One major cause of these discrepancies is that they relied on a global plate circuit with assumed rigid plates to restore the Pacific and Farallon plates and then used those as references to restore microplates of North American affinity.

The reconstruction of McQuarrie and Wernicke (2005) extended back to 36 Ma. Figure 18 compares our reference model NI reconstruction at 36 Ma, using the same map projection. We find rough agreement on the reconstructed positions of the Sierra Nevada block and other features in the western United States. However, their reconstructed positions for Baja California and Baja California Sur lie far east of ours, by distances ranging from 154 to 363 km. Consequently, their reconstruction of the future western Transverse Ranges of California is also far east of ours, lying where we place the future Peninsular Ranges. The principal reason for these differences lies in the different ways that we handled uncertainty about the amount of extension in the Mexico part of the Basin and Range Province. The primary constraint on their reconstruction of Baja California and Baja California Sur with respect to stable North America was a global plate circuit (Atwater and Stock, 1998) assuming rigid plates. In contrast, we compiled extensional faults (both high-angle normal faults and low-angle detachment faults) from the map of Muehlberger (1992) and from regional studies, and then assumed that unknown offsets on these extensional faults were similar to those that have been measured on other faults of the same class (Table 1). Then we essentially integrated the motion of Baja California by adding fault heaves along paths going southwestward from stable North America.

This discussion suggests that it will be important to test whether and how our reconstructions might be reconciled with global plate circuits. This subject is complex because, in addition to the usual random errors in reconstructed plate positions, it will be necessary to consider possible systematic errors due to distributed strain in each of the plates. We hope to address this in a future study.

We support the conclusion of Langenheim and Powell (2009) that the 18 km left step of the San Andreas trace north of San Gorgonio Pass is a geologically young feature, developed in conjunction with left slip on the Pinto Mountain fault. (See Ingersoll and Coffey, 2017, for additional discussion.) We also suggest that Pliocene thrust offset on the Dillon fault may have contributed to forming this tectonic knot. Unfortunately, our present finite-element grids are too coarse to resolve the details of this process; also, our use of automatic trace smoothing along active strike-slip faults during restoration is potentially a complicating element.

At 6 Ma, we successfully restore the Pelona Schist in the San Gabriel Mountains block opposite the Orocopia Schist in the Orocopia Mountains (Fig. 9); these exposures had a common origin in the Chocolate Mountains anticlinorium (Ingersoll et al., 2014; Coffey et al., 2019). Thus, we accept 6 Ma as the likely origination time for this branch of the San Andreas system. Following Ingersoll and Rumelhart (1999), we note that during 12–6 Ma, the most active branch in southern California was the San Gabriel fault with dextral offset of ~43 km at its northwest end. During 18–12 Ma, the San Gabriel fault may have had an additional 15 km of earlier dextral offset, but the dominant branch was the San Gabriel–Canton fault with ~27 km of dextral offset. (See Crowell, 2003, for additional discussion.) During both epochs, much of this dextral offset eventually reached the Cristianitos fault, which had ~32 km of right slip and ~10 km of orthogonal extension, which initiated subsidence of the Los Angeles basin (Ingersoll and Rumelhart, 1999). Evidence for some of the offset along the San Gabriel–Canton fault may be hidden beneath the younger strata of the Los Angeles basin. The Cristianitos fault merged southward with the dextral-transtensional Oceanside detachment system.

In the more linear part of the San Andreas fault system to the northwest of the modern Transverse Ranges, our reconstructions (Fig. 19) are consistent with realignment of the Pinnacles Volcanics against the Neenach Volcanics at ca. 21 Ma (Matthews, 1976; Burnham, 2009). However, this is not the total offset or the origination age for the central California segments of the San Andreas fault because Sharman et al. (2013) showed that an additional 50–100 km of earlier dextral slip is needed to restore the middle Eocene Butano Sandstone opposite its likely sources in the Mojave Desert region.

Since the pioneering plate-tectonic reconstructions of Atwater (1970), it has often been inferred that the San Andreas dextral fault system had its origin when the Pacific and North American plate first made contact with each other at ca. 28 Ma. This may be correct, but San Andreas origination at 28 Ma is not shown in our maps, because we chose not to let global plate-tectonic models (based on assumed rigidity of plates) influence our summaries of fault history; instead, we assigned the first chapter of San Andreas history to begin at 37 Ma, because that is the minimum age of the offset Butano Sandstone. This assignment does not indicate that slip on the San Andreas fault necessarily began at 37 Ma; it only dictates that slip could not have begun prior to or during deposition of the Butano Sandstone.

In fact, the primary change in tectonic style that we note, beginning ca. 24 Ma, is SE-NW extension, as documented in the Plush Ranch–Vasquez-Diligencia (PR-V-D) basins (Hendrix and Ingersoll, 1987; Tennyson, 1989; Ingersoll et al., 2014); soon thereafter, SW-NE extension dominated along the linked detachment faults of the Colorado River extensional corridor. This extension became so robust by 21 Ma that the western edge of our model locally had a westward velocity of 8 mm/a (Fig. 20) relative to stable North America and the future Colorado Plateau. The initial cause of extension may have been triple-junction instability (Ingersoll, 1982) at the northwestern (Mendocino) triple junction, which may have been nearby northwest of the PR-V-D basins. A likely contributing and accelerating factor was the post–28 Ma development of ever-widening slab windows within the subducted Farallon plate (Dickinson and Snyder, 1979; Wilson et al., 2005) that allowed hot asthenosphere to directly contact the base of North American crust, resulting in uplift, thermal weakening and volcanism (Bird, 1979).

A similar clarification is needed concerning the origination age of dextral shear in the proto–Gulf of California/Mar de Cortés. By imposing the common offset history shown in Table 2 to all component faults along the plate boundary, we achieve restoration of both the Jolla Vieja and Poway conglomerates to the southwest of their Sonoran source (Abbott and Smith, 1989) at 48–42 Ma. This restoration is shown in the supplemental file GMUM-S_CA_NI_0–48Ma_restoration.pptx (see ^{footnote 1}). However, our model history of proto-Gulf dextral offset at a steady rate of 2.1 mm/a during 46–12.6 Ma is not unique. Hypothetically, offset could have started later (e.g., 28 Ma) and proceeded faster (e.g., 4.6 mm/a) to achieve the same result.

The western Transverse Ranges of California (Santa Ynez, Santa Susana, and Santa Monica Mountains, and the Northern Channel Islands) have rotated clockwise since 18 Ma, as shown by many paleomagnetic studies. However, our model rotation is only ~70° (as in Hornafius et al., 1986, Wilson et al., 2005, and Hillhouse, 2010) rather than ~120° (as in Nicholson et al., 1994, and Ingersoll and Rumelhart, 1999). Given our algorithm's tendency to underfit its deformation targets, it might be natural to assume that this is another example, and that greater weight (or lower uncertainties) on paleomagnetic data would have given a different result with more rotation. However, this is not the case. Rotation of the western Transverse Ranges was only possible because of the extremely large extension on the Oceanside detachment fault during 18–6 Ma (Crouch and Suppe, 1993; ten Brink et al., 2000; De Hoogh et al., 2019), and rotation was limited by the amount of detachment slip. During our reconstruction, we inferred an upper limit of 140 km on the extensional component of Oceanside detachment slip, based on San Nicolas Island geology (flat-lying unmetamorphosed Eocene sandstone; Vedder et al., 1974), which cannot have been covered by the hanging wall of the Oceanside detachment fault. We set the Oceanside extension target to this upper limit and attained it, restoring present San Nicolas Island just west of present Rosarito, B.C., where there are similar Eocene alluvial deposits (Fig. 19). This restoration is highlighted in supplemental file GMUM-S_CA_NI_0–48Ma_restoration.pptx (see ^{footnote 1}). Therefore, our model provides the maximum possible rotation of the western Transverse Ranges.

This reconstruction implies that the sinistral faults now defining the northern and southern margins of the western Transverse Ranges (Santa Ynez, Mission Ridge–Arroyo Parida–Santa Ana and Agua Blanca faults on the north; Santa Rosa Island, Santa Cruz Island–Point Dume–Santa Monica–Hollywood and Raymond faults on the south) formed at azimuths of ~020, consistent with principal stress directions rotated just slightly counterclockwise from present. Therefore, we predict that those faults will be found to have had monotonically sinistral strike-slip during their histories, instead of forming as dextral faults and later reversing their slip during rotation.

Along with the primary restored features discussed above, slip on innumerable minor faults and deformation of many local geologic units have been restored in order to provide temporally and spatially coherent paleogeologic maps at 6 Ma intervals. These maps provide templates with which paleotectonic and paleogeographic models may be tested. In addition, our maps display unanticipated kinematic patterns that may inspire new analyses of previous models.

In utilizing Restore's maps, researchers are cautioned to remember the following qualifications. On average, specific fault offsets, cross-section reconstructions, and paleomagnetically determined rotations of the model underfit the target data. As demonstrated in the Neotectonic Calibration section, this systematic underfitting is partially compensated at regional scale by continuum strain, producing total crustal velocities that approximately match GPS velocities. As a result, we have confidence in the regional kinematics of the model, but many local reconstructed areas may not be consistent with local geological constraints. This characteristic is especially notable with vertical-axis rotational data; individual finite elements rotate according to paleomagnetic data within them, but larger crustal blocks consisting of many finite elements generally rotate significantly less than the data-constrained finite elements within them. Where geologic constraints dictate that crustal blocks rotated rigidly, and paleomagnetic data are only available for local parts of these blocks, our maps could be manually altered to reflect this interpretation. Because our paleogeologic maps have been created as transparently and objectively as possible, we encourage researchers to examine all of our source data, or to construct their own maps using Restore in order to satisfy local constraints.

As an example of manual alteration of our results to more closely correspond to both paleomagnetic and geologic data, Figure 21 includes a modified paleogeologic map for 48 Ma based on the following: (1) the southernmost Sierra Nevada (Tehachapi block) has been restored counterclockwise 44° (Dickinson, 1996; Dickinson et al., 2005, and references therein), although counterclockwise rotation of the westernmost part by 59° might be appropriate (McWilliams and Li, 1985); (2) adjoining Mesozoic plutonic rocks of the Salinian block (southwest of the future San Andreas fault) have also been rotated counterclockwise and translated eastward to be more directly south of the southernmost Sierra Nevada; this translation is justified because of the general underfitting of Miocene east-west extension in the Mojave–Colorado River extensional domains; (3) coastal areas west of the Salinian block have also been moved southeastward along with Salinia; (4) the Chocolate Mountains anticlinorium has been rotated counterclockwise and the San Gabriel Mountains block has been translated eastward in order to better align offset depositional, structural and basement rocks trends along the anticlinorium (Ingersoll et al., 2014; Coffey, 2019; Coffey et al., 2019; and references therein); and (5) the western Transverse Ranges and adjoining areas to the west have been translated northward to lessen the gaps between these terranes and Salinia. Each of these manual changes to the 48 Ma paleogeologic map is based on underfit paleomagnetic data and suspected unaccounted-for crustal extension, and/or unknown magnitudes of oblique slip along large strike-slip faults. The rigor and objectivity of the Restore paleogeologic map are retained at small scale; the manual modifications, which are based on multicomponent geologic reasoning, apply at large scale, where dictated by local relationships.

Apart from local inconsistencies, as illustrated above, these Restore reconstructions achieve our overall goals of transparency, objectivity, and internal consistency. Primary achievements include: (1) reversal of dextral slip on the San Andreas fault system, including the San Gregorio–Hosgri, Hayward, Calaveras, Rinconada, San Gabriel, Newport-Inglewood, Elsinore, and San Jacinto faults; (2) reversal of clockwise rotations of segments of the Transverse Ranges by different amounts at different times; (3) closure of the Inner Borderland of southern California by reversal of extension along the Oceanside detachment, seafloor spreading, and dextral slip; (4) closure of the Gulf of California/Mar de Cortés; and (5) reversal of crustal extension in the Basin and Range province of the United States and northern Mexico, including the Rio Grande rift. Our 48 Ma paleogeologic map is the most robust palinspastic reconstruction for southwestern North America at the beginning of Basin and Range extension, and provides a good base map for both kinematic and dynamic studies of the earlier Laramide and Sevier orogenies. However, one problem remaining for future resolution is whether these reconstructions can be reconciled with global-plate-circuit reconstructions of the Pacific plate, when allowing for distributed strain in all plates. We will reconstruct the kinematics and paleogeology of those orogenies in a companion paper extending back to 90 Ma.

This appendix explains expected changes through time in the maps of fault traces, and how they are approximated in version 4 of program Restore. Most examples are taken from the western conterminous United States or Gulf of California/Mar de Cortés.

Some faults have both dip-slip and strike-slip components, either in different geologic epochs, or simultaneously. The Restore database structure and the Restore algorithm keep these two components of offset conceptually distinct, and model them with distinct sets of data and code.

Such faults often appear in groups with little variance in rake (e.g., all normal, all thrust, or all strike slip within one tectonic province). Contrary to plate-tectonic principles, the majority of these fault traces are short, and often unconnected at their ends. Specifically, many dip-slip faults do not end at connections to transform faults. Therefore, these faults often have variable slips and offsets, which go to zero at each trace end. We do not try to represent this internal variation in our input database of fault offsets; it is understood that offset targets are maximum offsets, near the center of each trace. If the finite elements used in a Restore computation are locally smaller than the lengths of these fault traces, then an elliptical pattern of offset can be expected to appear naturally in output maps of velocity, offset rate, and total offset. This elliptical variation arises from the conflict between the (uniform) desired offset target and the opposing target of minimizing continuum deformation in the crust surrounding each fault.

In theory, the fault trace of any dip-slip fault should move with the hanging wall. Therefore, as time passes, a thrust fault buries and hides deformed rocks of the footwall. Likewise, normal faults expose extensive areas of deformed footwalls that may be capped by mylonites and/or slickensides. It is clear that the active fault trace of a detachment does not remain at the eroding breakaway trace, but lies where active hanging wall touches active footwall at Earth's surface, which is also where the fault dives into the mid-crust.

There are three problems with implementing this natural “hanging-wall” rule in Restore:

1. When faults are closely spaced (relative to finite-element sizes, so that there are several traces in each finite element), there may not be any adjacent unfaulted elements to define the hanging-wall velocity.

2. During retro-deformation, normal faults would subduct their footwalls, requiring any other faults (and paleomagnetic, stress-direction, and cross-section data) located there to be identified and removed from the database for all earlier timesteps. Likewise, retro-deformation of thrusts would educt their footwalls, requiring other (crossing) fault traces to be extended before use in earlier timesteps. This would lead to very complex code, with possible artifacts.

3. In cases where dip-slip faults do connect to strike-slip faults at their ends (e.g., Las Vegas Valley dextral shear zone; Duebendorfer and Wallin, 1991), the extent and activity of the strike-slip fault would have to be adjusted every timestep.

All of these complications resulting from an exact “hanging-wall rule” are too difficult to implement, so we opt for a simpler alternative: continental plate-interior fault traces with dip-slip motion are embedded in the deforming grid, so that they move and strain but do not change their topologic relations (including their connections) or their desired offset target histories over time. The only changes in the topology of the map of continental plate-interior dip-slip fault traces occur when a restoration calculation reaches geologic epochs at which certain faults had not yet formed, and their traces are automatically removed from all maps and calculations. (If one is uncertain about the initiation age of any fault, and wishes to keep it in the maps for early epochs, this can be achieved by assigning it a long-lasting early chapter of slip with a target offset of zero and a small uncertainty.)

For an important detachment fault of really large offset (e.g., Oceanside detachment in southern California), one can approximately simulate the desired motion of the trace with the hanging wall by defining wide fault-corridor elements that are placed so that the fault trace runs along the hanging-wall side of the corridor. Then, this trace will move only slightly and incidentally with respect to its hanging wall. We used this method to represent the Oceanside detachment during the 18–6 Ma part of our model.

Continental plate-interior faults with strike-slip offset are initially modeled in Restore by the same rule as dip-slip faults (discussed above): the digitized points that define their traces are embedded in the surrounding crustal continuum and move with it over time. That means that fault azimuths may change if the surrounding crust is rotating. Also, the length of the fault may increase if the surrounding continuum is extending or decrease if the surrounding continuum is shortening. But the topology of fault connections and close relationships (e.g., en-echelon) will not change.

Experience has shown that two additional algorithmic steps (and one data-preparation rule) are beneficial for reasonable handling of all strike-slip cases.

First, it is necessary to preserve the smoothness (i.e., continuous variation of strike along the trace) of strike-slip fault traces as we retro-deform the fault network. If kinks in these traces are not removed, they tend to become more extreme (and less plausible) during retro-deformation. This makes strike-slip more difficult and may cause the fault's offset target to be missed. Also, surrounding crust (and other adjacent faults) can become involved in this non-physical model artifact. To prevent this, we apply a small amount of diffusive smoothing to the locations of all strike-slip fault traces after each deformation step in which that fault was active. Taking some small, dimensionless coefficient ε (e.g., ε = 0.01), we replace the longitude and latitude coordinates of all interior fault-trace points along a strike-slip fault with the sum of (1–ε) times their previous values, plus ε times the average of the coordinates of the two surrounding trace points. This operation is not applied to the two end points of a strike-slip trace, so they are not moved. The asymptotic limit of diffusive smoothing (if carried too far) would be a strike-slip fault trace that would be an arc of a great circle. In practice, this would not be noticeably different from the arc of a small circle that is expected from plate theory for rigid microplates.

Another possible problem is presented by fault traces that cross in an “X” pattern, where at least one of the two traces has strike-slip offset. In plate-tectonic theory, and also in reality, we expect that such “quadruple junctions” or “X” intersections are usually unstable, and quickly break down into two triple junctions (“T+T”), which then move apart. A reasonable expectation is that the strike-slip fault will divide the other fault trace into two disconnected pieces (especially if it has higher offset rate). The difficulty arises because Restore does not have any special code to simulate quadruple junctions or triple junctions or to break up fault traces, so it would simply keep both fault traces internally connected, and the fast strike-slip fault trace would cause the other intersecting fault trace to be bent into a “Z” or “S” shape. This is especially damaging if the second (bent) trace also has strike-slip offset. To prevent this, we observe a data-preparation rule: when faults of large (or fast) offset cross, at least one of the faults (usually the slower one) should be divided into two segments (with distinct names, trace numbers and offset histories) during data preparation. Thus, we preemptively change the topology of most “X” intersections into “T+T” intersections that are easier to handle.

The last problem is that Restore does not model strike-slip faulting with an actual velocity discontinuity; instead, it approximates this with a high rate of simple shear distributed across each finite element that contains the active strike-slip fault. The union of all such elements is called a “fault corridor.” Where other fault traces enter this fault corridor (perhaps approaching a “T” intersection with the strike-slip fault, as described above), their ends would be locally bent by the strong simple shear in the corridor. To avoid this, we make a local and specific exception to the general principle that continental plate-interior fault traces are embedded in surrounding crust. Whenever possible, the Restore code extrapolates the model displacement of the last (or last few) trace points on the intersecting fault based on the motion of adjacent elements that are just outside the fault corridor.

By “spreading system,” we mean a chain of plate-boundary faults that alternates between spreading-rise segments (which contain both normal faults and dike swarms) and transform (strike-slip) faults. We consider a “spreading system” to be “symmetric” if each spreading rise contributes new crust equally quickly to both sides. As the system operates, it generates new oceanic crust imprinted with an orthorhombic texture of rise-parallel inactive dikes and normal faults and transform-parallel inactive fracture zones.

All of the fault traces comprising a symmetric-spreading system are in motion with respect to each of the two bounding plates. However, they share a common horizontal velocity, which is the average of the horizontal velocities of the two bounding plates. Therefore, a symmetric-spreading system can operate for millions of years without changing its form in map view, if watched from this average-velocity reference frame.

In program Restore, faults with the “symmetric-spreading-system” attribute do not move as if embedded in either of the adjacent plates. Instead, they move (together) at the mean horizontal plate velocity. To determine these velocities, the program looks 100 km away, into the interiors of the bounding plates, to sample the horizontal velocities at those points. At ridge-transform corner points, which are joins between adjacent faults in the system, the “look direction” bisects the angles between faults. The advantage of this algorithm is that it operates locally, always producing a reasonable result, and independent of any need for a global plate-motion model or defined plate names or explicit Euler poles. If the bounding plates are not rigid, but are deforming (because of numerous other faults which are simultaneously active), this algorithm will produce a similar mild and regional slow straining of the symmetric-spreading system.

Faults that are to be part of a “symmetric-spreading system” are so designated in the fault-traces (.DIG) file. Either end of a strike-slip fault (or both ends) can receive this designation. However, where a symmetric-spreading system ends, and connects to more irregular fault networks (of a continental interior), the outer end of the last transform might have the designation of “throughgoing master fault,” as described below.

We observe one data-preparation rule for strike-slip transform faults within a “symmetric-spreading system”: they are usually defined by only two digitized fault-trace points, one at each end. Therefore, they remain geometrically straight at all times, and the “diffusive smoothing of strike-slip fault traces” discussed in other sections of this Appendix does not apply to them.

It is often possible to estimate neotectonic slip rates for a long strike-slip fault (using dated Quaternary offsets, interseismic GPS velocities, and kinematic modeling programs such as NeoKinema). This is routinely done as part of seismic-hazard estimation projects.

However, total fault offset is NOT a simple integral of slip rate over time! (Or, more precisely, the total offset depends on exactly how the time integration is conducted and is therefore nonunique.) Even the initiation age of the fault may be undefined in many cases. For example, consider one of the northernmost segments of the San Andreas fault in California:

F7001R San Andreas 01 (offshore segment of 2011 CFM) dextral fault, CA.

This segment extends from (40.207°N, 124.343°W) in the southern part of Cape Mendocino southeast to an offshore point (39.148°N, 123.795°W) that is WNW of Elk CA. Its length is 131 km, with its NW end lying only 35 km from the Mendocino triple-junction (NA-PA-JF\NA), and its SE end is located 166 km along the trace from this same triple-junction. Baldwin (1996) reported dextral offset of 6.1–13.1 m on this segment since 547–367 years ago, for a rate of 22 ± 8 mm/a (which, however, may only sample 1–2 earthquake ruptures). Prentice et al. (2000) reported offset of 180 m since later than 13 ka, for a minimum rate of 14 mm/a. Using these data (and also local GPS benchmark velocities), the NeoKinema modeling program estimated a neotectonic dextral rate of only 7.5 mm/a (Field et al., 2013). Even though the proper neotectonic rate for this segment is uncertain, let us stipulate (for purposes of argument) that the real neotectonic rate is 14 mm/a (Prentice et al., 2000).

An observer standing on the northeastern side of the fault (in the Coast Ranges) and looking SW might argue that this segment formed sequentially as the Mendocino triple-junction passed by, with relative velocity of 14 mm/a. Therefore, the inferred age of this fault segment would vary from 11.9 Ma (on its SE end) down to 2.5 Ma (on its NW end), and its expected total slip would vary from 166 km (on its SE end) down to 35 km (on its NW end).

Yet an observer attached to the NE corner of the Pacific plate (e.g., in the Spanish Submarine Canyon area to the SW of Cape Mendocino) might argue that this part of the Pacific plate has been in contact with North America plate since 28 Ma (Nicholson et al., 1994; Atwater and Stock, 1998), so the age of initiation of this transform fault segment should be 28 Ma. Also, its offset since that time is that of the Pacific plate, which rotated 22.4° since 28 Ma about an Euler pole at (58°N, 287°E) with respect to North America (Engebretson et al., 1985, Table 3). This implies dextral fault slip of up to 1490 km (which should be reduced by the amount of relative movement between the Coast Range and North American interior; this correction is probably no more than 500 km).

Another point of view is possible: Burnham (2009) summarized compelling evidence for dextral offsets of 205 ± 30 km since 16.6 Ma (implying mean rate of 12.3 ± 1.8 mm/a) based on 13 offset geologic features that occur farther to the southeast along connected branches (e.g., the Pilarcitos, Peninsular San Andreas, and San Gregorio faults).

This is the paradox: When strike-slip faults have offsets comparable to (or exceeding) their lengths, then their initiation ages, total offsets, and slip histories are all unclear because there is no unique choice of reference frame for such calculations! Without some additional specific conventions, the best we can say is that these parameters lie between the bounds that would be estimated by observers on the left side, and those on the right side of the trace.

Segmentation of these faults in the geologic past is also problematic. In neotectonic studies, segmentation may be introduced because:

1. slip rate changes at a branch point in the map of fault traces;

2. rake angle changes (in concert with fault strike) at a fault bend;

3. seismic coupling changes (e.g., in the “creeping section” of the San Andreas); or

4. historic earthquake ruptures ended at certain points along the trace.

However, if segmentation is defined by branch points, then these have relative motion, and they may change their order and/or topology (or even appear and disappear) over geologic time. Fault strikes and rakes may change during retrodeformation. And seismicity parameters are unknowable in the geologic past. All of this means that there is no “natural” or universal definition of where neotectonic fault segments should be located in the geologic past; there can only be somewhat artificial, algorithmic definitions.

In Restore, any such fault is usually designated as a “throughgoing_master_fault” (in the fault-traces .DIG file). Such faults may be composed of multiple separately digitized segments, and thus may change their names (and trace #s) along strike, as long as there is no ambiguity about which two adjacent fault traces are meant to be connected at each join at each time. (That is, a “throughgoing master fault” is not allowed to branch, and results are unpredictable if the Restore user tries to create such a branch. Instead, branch faults with lesser offset rates should be treated as ordinary continental plate-interior faults, of type #1, as described above.)

Like continental plate-interior faults, throughgoing master faults are embedded in surrounding continuum, and initially move with it. Specifically, the digitized-coordinate points that define the fault trace move with the continuum. As narrow fault-corridor finite elements shear, these digitized points along the fault trace move at velocities which are between those of the left-wall block and those of the right-wall block.

Whenever strain becomes excessive (e.g., every ~0.8 m.y.), these fault-corridor finite elements are manually re-defined (by selecting new trios of nodes to define their corners), and new internal coordinates are computed for all fault-trace points within these new elements. Thus, the evolution of the master fault trace is unaffected (to first order) by this local re-gridding.

Some geologists might expect that evolving strike-slip faults tend to become straighter and smoother as slip increases. Therefore, during retrodeformation, one might expect them to become slightly rougher and less smooth. However, without special care (and coding), this process becomes an artificial instability in program Restore, with implausible results. To correct this, all active strike-slip fault traces are subjected to diffusive smoothing (as described above in section #1), proportional to their slip rates, and applied in small increments after their translation within each Restore timestep for which the fault is active.

Experience with Restore also shows that “matching” ends of adjacent segments of a throughgoing master strike-slip fault can become separated, which would lead to artificial hardening of the fault zone and other undesirable artifacts. Therefore, special code is added to average these segment ends to a common location after each translation timestep, preventing any separation at the join.

Experience also shows that segment-boundary joins may tend to become sharp fault bends, which also causes artificial resistance to slip and other unwanted artifacts. Therefore, additional special code is added to replace the (common) join location with the average of the locations of the two adjacent digitized points in the joined segments. This tames the bending instability and keeps the fault trace reasonably smooth in map view.

The collective result of these (somewhat arbitrary, but specific) choices in our kinematic algorithm is that the segments of a throughgoing master fault (and their associated offset targets and histories) remain well defined during any Restore computation. However, the user must anticipate that the segments are not fixed with respect to either the left-wall microplate or the right-wall microplate, but drift with respect to each of these, at intermediate velocities. Thus, their behavior is qualitatively (but not quantitatively) similar to the behavior of symmetric-spreading-system faults described in #2 above.

To return to our introductory paradox about the northern San Andreas fault near the Mendocino triple junction: This complex and unstable triple junction is left just outside the domain of our model, but is close to the northeast corner of the Pacific plate, which is southwest of the San Andreas fault and south of the Mendocino transform fault. Therefore, this plate corner and the implied nearby triple junction are restored southeastward at velocities which could be as high as the full PA-NA relative velocity. This progressively exposes the northeast-wall block of the northern San Andreas fault as the former continental margin, which was then part of the accretionary forearc of the Cascadia subduction zone. The digitized trace points that define the northernmost segment of the San Andreas fault are sheared (along with their enclosing fault-corridor elements) and translated southeastward at a slower velocity (about half of its slip rate). However, this northernmost segment progressively loses its southwestern-wall block, and thus becomes part of the perimeter of the model. Thus, this segment has only a short “life” in terms of its effects on the rest of the model, regardless of whether the user decides to keep or delete the associated fault-corridor elements. Meanwhile, the northeastern corner of the Pacific plate “overtakes” progressively older and more southeastern segments of the San Andreas fault during retrodeformation, so it never loses the ability to shear and move with respect to inland North America (until the time which the data editor has chosen as the initiation age of the whole San Andreas fault). This resolves the paradox, even if the resolution is not intellectually satisfying.

The University of California, Los Angeles, supported the first author with sabbatical research time in 1983, 1990, 1997, and 2004, as well as other research time and library and computational facilities. This work was supported by the National Science Foundation under grants EAR-9316169 and EAR-9614263, and by the U.S. Geological Survey (USGS), Department of the Interior, under USGS award number 01HQGR0021 to the University of California. The views and conclusions contained in this document are those of the authors and should not be interpreted as necessarily representing the official policies, either express or implied, of the U.S. Government. Bruce Luyendyk contributed a data set of selected paleomagnetic results for the southern California region. We thank John Fletcher, John Geissman, Craig Jones, and Mike Oskin for thorough reviews.