Abstract

Geologic slip rate determinations are critical to both tectonic modeling and seismic hazard evaluation. Because the slip rates of seismic faults are highly variable, a better target for statistical estimation is the long-term offset rate, which can be defined as the rate of one component of the slip that would be measured between any two different times when fault-plane shear tractions are equal. The probability density function for the long-term offset since a particular geologic event is broadened by uncertainties about changes in elastic strain between that event and the present that are estimated from the sizes of historic earthquake offsets on other faults of similar type. The probability density function for the age of a particular geologic event may be non-Gaussian, especially if it is determined from crosscutting relations, or from radiocarbon or cosmogenic nuclide ages containing inheritance. Two alternate convolution formulas relating the distributions for offset and age give the probability density function for long-term offset rate; these are computed for most published cases of dated offset features along active faults in California and other western states. After defining a probabilistic measure of disagreement between two long-term offset rate distributions measured on the same fault train (a contiguous piece of the trace of a fault system along which our knowledge of fault geometry permits the null hypothesis of uniform long-term offset rate), I investigate how disagreement varies with geologic time (difference in age of the offset features) and with publication type (primary, secondary, or tertiary). Patterns of disagreement suggest that at least 4%–5% of offset rates in primary literature are fundamentally incorrect (due to, for example, failure to span the whole fault, undetected complex initial shapes of offset features, or faulty correlation in space or in geologic time) or unrepresentative (due to variations in offset rate along the trace). Third-hand (tertiary) literature sources have a higher error rate of ∼15%. In the western United States, it appears that rates from offset features as old as 3 Ma can be averaged without introducing age-dependent bias. Offsets of older features can and should be used as well, but it is necessary to make allowance for the increased risk, rising rapidly to ∼50%, that they are inapplicable (to neotectonics). Based on these results, best estimate combined probability density functions are computed for the long-term offset rates of all active faults in California and other conterminous western states, and described in tables using several scalar measures. Among 849 active and potentially active fault trains in the conterminous western United States, only 48 are well constrained (having combined probability density functions for long-term offset rate in which the width of the 95% confidence range is smaller than the median). Among 198 active fault sections in California, only 30 have well-constrained rates. It appears to require approximately four offset features to give an even chance of achieving a well-constrained combined rate, and at least seven offset features to guarantee it.

INTRODUCTION

For about half a century, geologists have attempted to measure the slip rates of active faults by finding, documenting, and dating offset features. Such studies helped to test the theory of plate tectonics, and they continue to provide ground truth for those who model continental tectonics. These rates are also critical to estimates of future seismicity, which lead in turn to estimates of seismic hazard and risk, and in some cases to revisions of building codes. Considering that lives and large investments are at risk, the treatment of the uncertainties in these data has often been surprisingly casual. Government agencies that have tabulated slip rates and uncertainties have rarely specified the nature of the distribution for which they report the standard deviation (or other undefined scalar measure of uncertainty). They have often arrived at their preferred rates and uncertainties by deliberation in a committee of experts, which is an undocumented and irreproducible process. A disturbingly large fraction of the rates have been quoted as having uncertainties of exactly ¼ or ½ of the preferred offset rate, suggesting that the subject did not receive very serious consideration.

It might seem desirable to answer such questions with very intensive resampling of slip rates on a few representative faults. For example, the National Earthquake Hazard Reduction Program has funded multiple investigations of the San Andreas fault. However, there are two obstacles: first, exposed geology typically only provides a limited number of datable offset features (if any) along any given fault trace. Second, slip rates definitely vary in time (e.g., during earthquakes), and may also vary in space (along the trace), which makes it very difficult to conclusively falsify any single measurement with another single measurement. Many authors prefer to resolve discrepancies by increasing the number of free parameters in the slip history of the fault.

Consequently, a purely frequentist approach to determining the uncertainty in geologic slip rates is not practical or correct. We must rely on a combination of: (1) redefinition of the objective, from a complex slip history to a single long-term slip rate; (2) Bayesian approaches, in which prior assumptions about the shapes of distributions substitute for multiple sampling; and/or (3) bootstrap methods, which use properties of the distributions of offset distance, offset age, or offset rate for all active faults in some class to help estimate the corresponding distribution of a particular fault. One objective of this paper is to begin a formal discussion about which redefinitions, prior assumptions, bootstrap methods, and computational paths the community of geologistsand geophysicists might choose to support as a standard. Standardization would be helpful because it would: (1) provide guidance to investigators; (2) increase reproducibility of conclusions, and (3) permit automated retrospective revision of geologic offsetrate databases if prior assumptions should need to be changed, or if systematic errors in the geologic time scale should be discovered in the future.

One of the challenges we face is to mitigate the adverse effects of three kinds of misleading data, which cannot usually be identified in isolation or in advance: fundamentally incorrect offset rates (due to, for example, failure to span the whole fault, undetected complex initial shapes of offset features, or faulty correlation in space or in geologic time); inapplicable offset rates (which are correct as averages over their lengthy time windows, but misleading when considered as neotectonic rates); and unrepresentative offset rates (very small or zero rates measured at or beyond the extent of the fault trace which is active in the neotectonic era). I propose that the first two problems can be handled by developing bootstrap estimates of their frequency, and then merging comparable offset rates (along one fault) with a formula built to reflect these probabilities. I suggest handling the third issue by using the small- or zero-offset data to redefine the length of the active fault.

The method advocated here has the following nine steps.

  • (1) Estimate the probability density function of one scalar component of the far-field cumulative offset since a particular geologic event, including uncertainty due to plausible but invisible elastic relative displacements which leave no geologic record.

  • (2) Estimate the probability density function of the age of the geologic event associated with this offset (which will frequently be a smoothed-boxcar or other non-Gaussian distribution).

  • (3) Convolve these two distributions to obtain the probability density function for the long-term offset rate for this particular pair of offset features.

  • (4) Define a scalar, dimensionless measure of the disagreement between two long-term offset rate distributions determined for the same component of slip on the same fault train (defined below).

  • (5) Identify strong disagreements between multiple rate estimates for a single fault train, and calculate how their frequency in a given large region varies with age of the offset features and with the type of literature source

  • (6) Estimate the fractions of published geologic offset rates which are incorrect or unrepresentative (for each type of literature source) when the offset feature is young.

  • (7) Estimate the additional fraction of published offset rates that are inapplicable to neotectonics, as a function of the age of the offset feature.

  • (8) Considering results of steps 5–7, merge all offset-rate probability density functions from one fault train (possibly including incorrect, unrepresentative, and inapplicable ones) to provide the best combined estimate of the long-term neotectonic offset rate (under the assumption that it does not vary along the trace).

  • (9) Express this distribution in terms of simple measures (mode, median, mean, lower and upper 95% confidence limits, and standard deviation) and present these in tables.

The Bayesian aspects of the program are most apparent in steps 1 and 2, while the bootstrap aspects are most apparent in steps 5–7. Steps 3–4 and 8–9 are purely mathematical. However, at several points it will be necessary to deal with incomplete information, as when only a lower limit on an offset, or only an upper limit on its age, is available. In these cases, I rely on a prior estimate of the probability density function for offset rate that is determined by bootstrap estimation based on similar faults in the region. The assumptions necessary to justify this are described in the next section.

This paper has a rather narrow focus on determining long-term offset rates of faults only by use of offset geologic features (and groups of features). Several other valid approaches are available: (1) use of geodetic estimates of differences in benchmark velocities across faults; (2) use of plate tectonic constraints on the total of the vector heave rates for a system of faults composing a plate boundary; (3) local kinematic-consistency arguments that extend a known slip rate from one fault to those that connect to it; and (4) use of instrumental and historical seismicity and/or paleoseismicity. There is extensive literature on numerical schemes for merging these approaches, and it is definitely true that using a variety of schemes will reduce uncertainty. I am also involved in such studies, in which we use our kinematic finite-element code NeoKinema (e.g., Bird and Liu, 2007). However, merging geologic, geodetic, plate tectonic, and perhaps seismicity data should only be attempted after properly characterizing the uncertainties in each. This paper addresses the more limited task of determining best-estimate offset rates and their uncertainties purely from offset geologic features.

BASIC ASSUMPTIONS

1. I assume that geologists can accurately distinguish tectonic faults (those which cut deeply into the lithosphere) from surficial faults surrounding landslides (which merge at a shallow depth) and from surficial faults associated with sediment compaction, groundwater withdrawal and recharge, or magma chamber inflation and deflation (which root in a subsurface region of volume change). I assume that only tectonic faults are included in the database.

2. I assume that the senses of predominant offset reported for an active fault (dextral or sinistral and/or normal or thrust) are correct. The reporting geologist typically has access to many scarps and offset features that are not datable, but can be described as “young” with high confidence. Also, regional fault-plane solutions and/or other stress-direction indicators provide guidance.

3. I assume that tectonic faults have motion that is geometrically and physically related to relative plate velocities, so that long-term offset rates that are orders of magnitude faster or infinite are not plausible.

The implication of these assumptions is that long-term offset rates of tectonic faults are defined to be non-negative, and that a prior probability density function for long-term offset rates of active tectonic faults can be estimated from research in a particular plate-boundary region.

For procedural reasons, it may be desirable to add another qualifying assumption.

4. Data to be processed by this proposed algorithm should ideally come from the peer-reviewed scientific literature (including theses, but excluding abstracts). This helps to protect data quality in several ways. (1) The numbers describing offset, age, and uncertainty will be the authors' final best estimates, recorded in an archived source. (2) Each entry will have been screened for elementary errors in logic. (3) Data will be certified as having significant information content. As a counterexample, imagine that someone wished to treat a recent snowfall as an overlap assemblage, and entered hundreds of “data” showing that each of the faults in California was motionless in that particular week. A few such “data” do no harm, as this algorithm will be designed to return the prior distribution of long-term offset rates in such cases. However, it would be very undesirable for the combined long-term offset rates of faults to be too heavily weighted toward this inherited prior distribution, when other data of higher information content are available. One opportunity for future improvement in this algorithm would be to define and quantify the subjective notion of information content and then to derive its proper role in weighting.

OFFSETS AS COMPONENTS OF SLIP

Assume that one side of a fault (e.g., footwall) is taken to be fixed, to provide a temporary and local reference frame for displacements and velocities. The displacement of the moving side is the slip vector. Slip can be considered as consisting of a vertical component (throw) and a horizontal component (heave). The heave is a two-component horizontal vector, which can be further subdivided into a fault-strike-parallel component (strike slip) and a perpendicular horizontal component (closing or opening, i.e., or convergence or divergence). Here I use the generic word offset to stand for any of these scalar components of fault slip. It is only meaningful to compare long-term offset rates that refer to the same component of slip rate. In practice this component is usually the strike slip (for predominantly dextral and sinistral faults) or the throw (for thrust faults and high-angle normal faults) or the divergence (for low-angle normal faults and/or magmatic centers of crustal spreading). Occasionally the convergence component can be estimated for low-angle thrusts by use of borehole, tunnel, and/or seismic reflection data. Data on divergence or convergence can be compared with data on throw by measuring or assuming the mean dip of the fault within the seismogenic layer. In this study, cases of oblique slip (mixed strike slip and dip slip) are usually handled by treating these two components as separate processes (which happen to occur simultaneously on the same fault trace).

In order to determine slip or offset, it is theoretically best to map (or excavate) a pair of offset piercing points that were separated when fault slip disrupted an originally continuous piercing line. Piercing lines include rims of impact craters, shorelines (formed during a brief highstand), terminal glacial moraines (formed during a brief glacial maximum), and lava flows and debris flows that were confined to narrow straight valleys crossing the fault trace. If the topography already included fault scarps at the time of the formation of any of these features, they may have formed with initial kinks. Then, there is danger of erroneous interpretation.

A different kind of piercing line may be defined as the intersection of two planar features of different ages, if the same pair of features is found on both sides of the fault. For example, flat-lying sediments may be intruded by a vertical dike, with both planar features cut by the fault. Or, gently dipping alluvial fan surfaces may be truncated laterally by steep cut banks on the outside of channel meanders. Such intersection piercing lines present the risk of another type of misinterpretation: the throw will be recorded beginning with the formation of the quasi-horizontal feature, but the strike slip will be recorded only after the formation of the quasi-vertical feature. In such cases, it is best to treat the throw and the strike-slip components of the slip as separate problems, each with their own constraining data and resulting model. Then, separate ages may be assigned for the critical geologic events that initiated recording (e.g., sedimentary bed age versus dike or cutbank age). Or, the same age may appear as a best estimate for one offset component, and as an upper limit for the other offset component (e.g., cosmogenic nuclide age of a fan surface truncated by a cutbank).

Where the geology does not present piercing points created by the rupture of a piercing line, it is still possible to measure some kinds of offsets using pairs of fault-surface trace lines, which are the intersections of originally continuous planes with the fault surface. Offset of a quasi-vertical and fault-perpendicular dike, vein, or older inactive fault can provide an estimate of strike slip. Offset of a quasi-horizontal sedimentary bed or erosion surface can provide an estimate of throw.

One special note is necessary concerning throws of thrust faults. New thrust faults typically propagate upward, and prior to break-out they may be blind thrusts expressed as monoclines at the surface (especially if sediments are present). Even after breakout, the vertical height of the monocline is often greater than the throw measured directly across the fault trace, because the height of the monocline incorporates the entire zone of anelastic deformation. If the purpose of the database is to support estimates of seismic hazard, then it is probably prudent to assume that the throw of thrust faults at several kilometers depth (where large stress drops may occur during earthquakes) is given by the height of the monocline, not by the smaller throw seen directly across the trace in outcrop. This phenomenon of distributed anelastic strain around the shallow portions of thrust faults is conceptually separate from the following discussion of distributed elastic strain around all seismic faults.

SHORT-TERM VERSUS LONG-TERM OFFSET RATES

A majority of active continental faults are seismic, so they slip as much as 15 m in a few seconds, and may then not slip at all (at the surface) for hundreds to thousands of years. The modeling methods of Bennett (2007) deal with time-dependent slip rates, but work best in cases where the variation in slip rate has been smooth, so that it can be described by approximately as many parameters as there are offset features. Also, Bennett's derivation assumes that unmodeled effects have Gaussian distributions, which is not a good description of slip-rate variations within one seismic cycle. Therefore, this study takes a different approach: instead of attempting to recover the history of slip rate on a fault, I attempt to recover only a single long-term rate (defined below) which is presumed to be steady over time scales of a few million years. This less ambitious objective is actually difficult to reach.

In many cases, geodesy has shown that benchmarks more than 10–30 km from active faults have relative velocities that are steady through time, and that differences between far-field and near-field displacements are accommodated by changes in elastic strain. Let F(t) be some offset (as a function of time t) measured at the fault trace (or from fold amplitude, in the case of thrust monoclines), while D(t) is a distant (far field) relative displacement component expressed as the same kind of offset, and E(t) is the cryptic distributed offset due to changes in elastic strain along the same fault (but excluding any displacement due to strain fields of other faults). In this notation,  
formula
where D(t) is monotonically increasing and curvilinear in time, F(t) is non-decreasing but highly nonuniform, and E(t) is an irregularly oscillating “saw tooth” (or “fractal saw tooth”) function that is bounded within a range of ∼15 m or less (Fig. 1). In this paper, short-term offset rate (S = ΔFt, for some interval of time Δt that usually ends at the present) is used for the quantity that has typically been derived from offset features in the literature. Long-term offset rate L = ΔDt will be the goal of statistical estimation because it has many fewer parameters, and can perhaps even be approximated as independent of time within the Quaternary history of the western United States (although that remains to be demonstrated). Then,  
formula

Because long-term offset rate L is central to this paper, but rather abstract, the following are some idealized cases in which L has a simple physical meaning.

  1. If a fault is aseismic and creeps very regularly (e.g., due to strain-rate-strengthening rheology in the gouge), then ΔE = 0 and L = S.

  2. If the active fault is planar, has a uniform slip-rate vector (along its trace), and is sufficiently isolated from other active faults, then one can use geodetic benchmarks far from the fault (where ΔE ≅ 0 according to dislocation modeling) to measure ΔD and thus L directly.

  3. If we were willing to install stress meters beside creep meters and record the history of shear stress changes and fault movement for several centuries, we could then choose any two different times with equal shear stress, assume that ΔE = 0 for that particular Δt, and obtain L directly from the ΔF measured by creep meters in the same time window.

Fortunately, the maximum possible difference between short-term offset rate and long-term offset rate tends to decrease with increasing measurement duration Δt:  
formula
where operator ‖ ‖ indicates an absolute value, and sup (‖ΔE‖) (the maximum of the absolute value of the change in elastic displacement) can be estimated based on historic earthquakes on other faults of similar style.

PROBABILITY DENSITY FUNCTIONS AND FUNCTIONAL COMBINATIONS

We now begin to treat ΔF, ΔE, ΔD, and Δt as random variables, which all refer to one particular offset feature on one particular fault. They will not always be random in the frequentist sense that they display variation in repeated trials; many are random in the Bayesian sense that a range of values is plausible with varying degrees of belief (in a community that shares certain assumptions and a common pool of literature) (Cowan, 2007).

Each variable (taking ΔF as an example to demonstrate notation) will be characterized by a cumulative distribution function PFf), where P() indicates a dimensionless probability, and by an associated probability density function (pdf) denoted by pΔF(f):  
formula
 
formula
 
formula
 
formula
In this paper I will employ some ad hoc functions that are almost pdfs, in the sense that their integrals are dimensionless numbers of order unity. For compact notation, it is very convenient to convert these to pdfs by defining a normalizing operator N() for any function that has the physical units of a pdf, and that adjusts this function to an integral of unity, as required by equation 4d. That is, for any function y(x) that has a finite integral, and whose integral is nonzero, and assuming that the physical units of y and of x are respectively reciprocal,  
formula

Note that this normalizing operator does not cause any change in the physical units of its argument function, because the denominator in equation 5 is restricted to be dimensionless.

We also need to combine random variables by elementary operations like addition, subtraction, and division, so it is useful to review how their pdfs will combine. Let X and Y be two independent random variables, with associated pdfs of pX (x') and pY (y'). Let H(X,Y) be an algebraic function that is both differentiable and monotonic (hence, invertible) with respect to each of X and Y. Then, the pdf of H, symbolized by pH (h'), can be expressed as either of two weighted convolutions:  
formula
or  
formula
where all three terms inside each integral are uniformly non-negative. See the Appendix for proof.

OFFSET AT THE FAULT TRACE, ΔF

Stated Offset is a Best Estimate

The measurement of ΔF is made with a ruler (either on the ground, or on a map or section). This is a well-studied type of measurement, for which it is traditional to assume a Gaussian (normal) distribution. However, under the assumptions of this project, ΔF is not permitted to be negative, while the long tails of a Gaussian distribution would permit that possibility. Therefore I assume a truncated and renormalized Gaussian distribution:  
formula
where ΔFn is the nominal or most likely value and σΔF is its standard deviation.

In most cases, the uncertainty in the offset (σΔF) is dominated by the somewhat subjective decision of where to place the ruler, rather than by the mechanical operation of reading it. This is especially an issue when the offset ends of the piercing line no longer reach the fault trace because of erosion. It would be desirable for investigators to ask coworkers and colleagues to make independent measurements of offset “at” (i.e., “subjectively projected to”) the fault trace, so that ΔFn and σΔF could be estimated from a sample population of more than one (a frequentist approach).

When the investigation is no longer active, and the nominal offset ΔFn is obtained from the literature or from a database, there are several alternatives for estimating σΔF. If the authors cite a standard deviation explicitly, it should be honored. If the authors cite a range for the offset (e.g., 5–9 m, or equivalently 7 ± 2 m), it is important to know what level of confidence was intended. For example, if the range was a 95% confidence range, then σΔF= (ΔFmax – ΔFmin)/4 (e.g., 1 m in this case). Unfortunately, many authors do not state the level of confidence, and another possibility is that their intention was to present “ΔFn ± σΔF,” in which case σΔF = (ΔFmax– ΔFmin)/2 (e.g., 2 m in this case). In cases where no clarification of intent is available, I adopt a compromise rule that assigns σΔF = (ΔFmax– ΔFmin)/3; that is, the quoted range will be interpreted as an 87% confidence range.

If no standard deviation and no range are quoted (or if this information has been lost during database compilation), then the only recourse is infer an implied standard deviation from the position of the least-significant nonzero digit in the quoted nominal offset ΔFn. Implied standard deviations can be defined as one-half of the increment caused by increasing the right-most nonzero digit by unity. For example, an offset stated as either 1 km or 1000 m would be assigned an implied standard deviation of 500 m, while an offset of 1.3 km or 1300 m would be assigned an implied standard deviation of 50 m. Implied standard deviations always have magnitudes of 5 × 10i of the current length units, where i is an integer.

Occasionally, a study will suggest two alternate pairings of piercing points as distinct possibilities. For example, either of two headless alluvial fans (of similar age) on the downstream side of the fault trace might be candidates for pairing with a canyon on the upstream side of the fault. One might consider using a bimodal probability density function pΔF in such cases, and the equations proposed below are general enough to handle such a model. However, when two possible pairings are suggested, we must also worry that these two cases may not be exhaustive, and that these two possible pairings may have been selected from among many, based on prior estimates of the offset rate for that fault. Such circularity would be impossible for regional modelers to unravel, and with the involvement of naive intermediaries it could easily grade into pathological science. To prevent this, it may be wiser to simply disregard any offsets for which bimodal (or multimodal) pdfs have been suggested.

A deep question that arises in this context is, how often have authors selected the offset features that they chose to report in the literature based on conformity with prior beliefs about the expected slip rate? Such selection might be completely unconscious if the geologist knew the age of the potential offset feature at the time of mapping. Fortunately, it is more typical that the age is determined only by time-consuming lab work. Then, any selection is at least conscious. Whether any investigator would be willing to report applying conscious selection is also problematic.

Stated Offset is an Upper Limit

Offsets are sometimes reported only as upper limits, with no best estimate. This occurs when the offset feature spans more than one active fault, but is not seen between the faults (due to erosion or burial). It is reasonable to assume that quasi-parallel faults of the same style have the same sense of slip, so the total offset serves as an upper limit ΔFmax for each fault individually. According to the assumptions of this study, zero is the default lower limit (because we assume that we know the sense of slip of each active fault). Rather than assuming a rectangular uniform (boxcar) pdf for such cases, I recognize that the upper limit is at least slightly flexible (subject to uncertainty), and use a smoothed-boxcar distribution with a complementary-error-function (erfc) shape around the upper limit:  
formula
where the standard deviation of the upper limit, σΔFmax, is either stated by the authors or estimated by one of the methods listed above. If the authors have indicated a range for the upper limit, its standard deviation is one-third of the quoted range. Otherwise, an implied standard deviation is determined from the position of the least significant nonzero digit in the number.

Stated Offset is a Lower Limit

Offsets for which only lower limits ΔFmin are reported present a greater challenge. This happens most often with high-angle normal faults, where the topographic height of the scarp at the mountain front gives a lower limit on the throw. Actual throw would typically be greater if we could measure the erosion of the mountain footwall and the sediment thickness on the valley hanging wall, both of which reduced the scarp height before it could be measured. Lacking an upper limit, we cannot write a one-sided pdf analogous to equation 9, because the function would have an infinite integral, and the normalizing operator could not be applied. This is a situation in which I appeal to prior knowledge to put loose but reasonable bounds on an otherwise open ended distribution. Assuming that we have a prior estimate
\(\mathit{p}_{\mathit{L}}^{prior}\)
of the pdf of long-term offset rates for this type of fault in the same region, I create a generic prior estimate of the offset. If Δt were precisely known, this would be:  
formula
The elastic term is dropped for simplicity, because it has zero mean, and because it is no larger than ∼6 m, which is not significant relative to the height of mountain ranges. Now, given that both Lprior and Δt can be represented by pdfs, and applying equation 6 to equation 10, the pdf of the prior fault-trace offset is:  
formula
[ where functional definition of the pΔt(t'), the pdf for the time interval Δt, is deferred to the section below]. Then, after application of the lower limit, the recommended form is (erf is error function):  
formula

Stated Offset is Zero

Observations that an active fault is covered by an unfaulted overlap formation of known age can be quite valuable. This is one of the few ways to directly falsify incorrect offset rates from imagined offset features adjacent on the same fault section. Or, the point of no offset may serve to document variation in long-term offset rate with distance along a fault trace. In tables and algorithms (like this one) that are primarily concerned with offset features, overlap assemblages might appear as offsets of 0 m, or (more often) as positive offsets that have a constrained ending time before the present. The associated pΔF for the fault-trace offset during post-overlap time is a delta function at the origin. In numerical computations a delta function is inconvenient, and I have approximated it with one of several distributions of finite width. If a standard deviation is explicitly stated for the zero offset, then a truncated and normalized Gaussian distribution is used (equation 8, with ΔFn= 0). If an upper limit on the nominally zero offset is stated, then that is used to create a uniform (boxcar) distribution extending from zero to the upper limit. If neither standard deviation nor upper limit is stated, then I infer a minimal amount of geologic imprecision:  
formula
where ϵ = 0.01 m. The basis for this arbitrary value is that an overlap assemblage might easily have ∼100 microcracks/transect of ∼100 μm offset each, and yet appear unfaulted to a field geologist. As discussed in a previous section, observations of zero offset could potentially become very numerous if they included every paved road crossing a fault, and there are practical reasons for limiting such data to only those with Δt of at least 100–1000 yr, data that add meaningful amounts of information.

CRYPTIC ELASTIC OFFSET, ΔE

Unlike ΔF and ΔD, cryptic elastic offset ΔE is not assumed to be positive. It will often be negative, whenever an earthquake (or multiple earthquakes) during the measurement time window caused fault-trace offsets exceeding the long-term expectation (ΔF > LΔt), so that regional elastic strain changes were net unloading rather than net loading. If E(t) is bounded by upper and lower limits (related to physical limits on the magnitude of shear stresses, and the geometry of the fault and elastic lithosphere), and if ΔE = E(t2) – E(t1) is measured between fiducial times that are both randomly selected, then the expectation (mean, first moment) of ΔE is zero.

Another question is whether the fiducial times have really been selected without regard to earthquake times. One can argue that some debris flows were caused by earthquakes on the fault, which they then crossed to create a post-seismic piercing line. One can also argue that faults are more likely to be studied just after earthquakes. I have no solution to offer, except to make a hopeful assumption that any biases caused by these post-seismic selection pressures on t1 and t2 will tend to cancel.

For long time windows associated with large offsets, meaning those that are likely to have encompassed multiple earthquakes, there is little likelihood of correlations in phase relative to the “earthquake cycle” (if such a thing even exists) at times t1 and t2. Therefore, the chances of having a net gain in cryptic elastic displacement are equal to the chances of having a loss of the same size. In this long-time limit, the probability density function of cryptic elastic offset (pΔE) can be estimated by mirroring-and-halving the probability density function of coseismic offsets in historic ground-breaking earthquakes (pC) on faults of the same type (e.g., only strike slip, only thrust, or only normal):

There are some issues raised by this approximation. First, the conditioning event that causes a coseismic offset to be measured and published has not been clearly defined. This could cause our estimate of pC to be inaccurate near the origin (small offsets). However, if we focus on the importance of estimating the tails of pC and pΔE, and note that the most dramatic scarps and largest earthquakes on land are always studied, this issue seems less serious. It is also necessary to assume that the historic period has been long enough to display the largest earthquakes that can occur in each tectonic setting. Bird and Kagan (2004) argued that this seems to be true, at least when earthquakes are compiled globally for a century.

Wells and Coppersmith (1994) created a global compilation of slip in historic ground-breaking continental earthquakes (excluding subduction events). I use maximum slips (rather than average slips) from their Table 1 (Wells and Coppersmith, 1994) because these are more frequently reported, and because they are more relevant to determining the tails of the pC and pΔE' distributions. I divide them into three populations of dominantly strike-slip, thrust, or normal-faulting earthquakes. When each population is plotted as a cumulative distribution (Fig. 2) they seem to have an exponential character:  
formula
 
formula

The median of the maximum slips is probably the most stable measure of these small samples: 2.0 m for strike slip, 1.3 m for thrust, and 1.15 m for normal earthquakes. Using these to determine the λ values results in 0.347 m−1 for strike slip, 0.533 m−1, for thrust, and 0.60 m−1 for normal. If the offset being studied is not the full slip, then these values would be geometrically transformed. For example, if the offset in question is the throw of a thrust fault, and the assumed dip of thrust faults is 20°, then the median maximum offset in historic thrust earthquakes would become a median throw of 1.3 m × sin(20°) = 0.44 m, and λ would correspondingly increase to 1.56 m−1.

Simple symmetry of pΔE breaks down for short time windows Δt. When the time window is so short that the number of expected earthquakes is <1, the positive (loading) side of pΔE is peaked around modest values of LΔt, but the negative side of the distribution has lower probability densities extending further from the origin and representing the small chances of large coseismic unloading of elastic strain. (The two halves of the pdf still balance so that the expectation is zero.) In this limit, I use the prior distribution of long-term slip rates
\(\mathit{p}_{\mathit{L}}^{prior}\)
(on faults of the same type) as a proxy for the unknown long-term rate of the fault in question, and convolve this with the pdf of the time available for strain accumulation, so that:  
formula
where k is a constant numerically determined to give pΔE an expectation of zero. Note that on the negative side (e' < 0) this pdf has the same form as in equation 14, which expressed the long-term limit.
For time windows Δt that are neither very short nor very long, a gradual transition in the positive side of the pdf is desirable. I propose as a working hypothesis that the complement of the cumulative probability of a certain positive elastic offset (which is the probability of an even larger elastic offset) decreases with offset as the product of (1) the estimated probability that such an offset could be accumulated in the time available (based on the prior distribution of long-term rates), with (2) the probability that the future seismic conversion of such an offset (ΔE → ΔF) would produce a coseismic offset consistent with the historic record:  
formula
The form of the pdf for positive elastic offsets can then be determined by taking the derivative (which becomes fairly complex because Lprior, Δt, and C are all random variables):  
formula
Finally, the positive and negative sides of the pdf are properly weighted and jointly normalized:  
formula
where k is computed (as in equation 16) to give an expectation of zero. As Δt increases there is a gradual transition from the growing distance scale of LΔt to the limiting distance scales of pC on the positive side of pΔE. The transition will take place (for any given small positive e') over a range of Δt because of the uncertainty in L and the fact that it has been replaced by the prior distribution as a proxy. However, for very large positive e', greater than the offsets in historic earthquakes, pΔE will be nearly zero regardless of Δt, because of the first factor in the last line of equation 18.

DISTANT (FAR FIELD) OFFSET, ΔD

Because distant (far field) offset ΔD is the sum of ΔE and ΔF (equation 1), its pdf is obtained from the convolution (equation 7) of pΔE with pΔF However, there is also a requirement (by initial assumption) that ΔD is positive. Thus, the convolution is truncated and normalized:  
formula

The complete generality of this convolution equation permits any of the different pdfs that were previously proposed (for different cases of available information) to be incorporated. A real example from the Teton normal fault of Wyoming is shown in Figure 3. Note that in other cases where the length scale of ΔF is much greater than the length scale of ΔE, the form of the convolution based on equation 7 is numerically superior to the alternative form implied by equation 6, which would be inconveniently stiff and difficult to compute precisely.

AGE OF THE OFFSET FEATURE, Δt

Stated Age is a Best Estimate

The offset features that are simplest to date are igneous dikes, lava flows, and pyroclastic flows, because they contain minerals that formed, and isotopic systems that blocked, at the time when the piercing line was established. I disregard as insignificant the small span(s) of years between laboratory fiducial time, laboratory measurement, publication of the measurement, and analysis of the published data. In these cases, the pdf for the age of the offset feature [pΔt(t')] can be approximated as equal to the pdf for the laboratory determination of radiometric age [
\(\mathit{p}_{{\Delta}\mathit{t}}^{lab}\)
(t′)]. These are estimated to be Gaussian, and characterized by a nominal or best-estimate (mode = median = mean) age [
\({\Delta}\mathit{t}_{n}^{lab}\)
] and by a standard deviation, which will depend on repeatability, instrument calibration, blank level, and any complications in the sample history (such as Ar loss or Pb loss) that can be handled by established methods. I truncate the age distribution so that it does not extend to negative ages:  
formula

If the authors quote

\({\sigma}_{{\Delta}\mathit{t}}^{lab}\)
explicitly, it should be honored. If they only quote a range of ages, then (as I did in the case of range of offsets) I will use the rule
\({\sigma}_{{\Delta}\mathit{t}}^{lab}\)
= Δtmaxtmin)/3. If only a bare age is stated, then the implied standard deviation will be defined to be one-half of the value of an increment of one in the rightmost nonzero digit; for example, quoted ages of either “1.3 Ma” or “1300 ka” would have an implied standard deviation of 50 k.y.

In other cases, volcanic strata are not directly dated, but are correlated with strata elsewhere that have been extensively studied and accurately dated (e.g., Bishop tuff or Peach Springs tuff in eastern California). While such correlations are not always correct, this is not a situation that is easily described by a probability density function for the age of the bed. Instead, I consider herein the possibility of faulty correlations in geologic time as a potential cause of fundamentally incorrect offset rates.

Very young ages (younger than 100 ka) of sedimentary layers based on 14C or cosmogenic nuclides tend to be biased toward the old side. This is because the event being dated (fixation of atmospheric 14C into wood, or initial exposure of a future rock clast on a weathering hillside) preceded the deposition of these samples in the fault-covering stratum. This upward bias in the age could create a downward bias in offset rates, if not recognized. In well-funded studies, this problem is sometimes mitigated by quoting the youngest age found in a group of analyses of samples from a single stratum, or by more sophisticated methods taking into account prior knowledge (e.g., Biasi and Weldon, 1994). Even then, some uncertainty and bias remain. This is similar to the problem of determining extinction ages from a finite distribution of discrete fossils (Weiss and Marshall, 1999).

Often this effect is uncontrolled, and an uncertain amount of inheritance (Δtinheritance) occurred between the beginning of time recording and the formation of the piercing line. As  
formula
it follows that we should broaden the pdfs of many laboratory 14C and cosmogenic nuclide  
formula
I assume a simple one-parameter exponential model for the distribution of inheritance:  
formula

The study of Biasi and Weldon (1994) is an excellent resource for estimating the ω of 14C inheritances in California, because they compared 106 laboratory ages (all younger than 2 ka) to posterior age distributions constrained by stratigraphy, sedimentation rates, and one historic earthquake. In the most extreme case (W125ac) the inheritance of one sample appeared to be ∼440 yr. However, the median apparent inheritance of strata (for those strata where it was positive) was only 51 yr. I use this value to estimate ω ≅ 0.0136 a−1 for 14C in California.

This issue is more serious for cosmogenic nuclides, because erosion is a slow process, and rock clasts are more durable on the surface than charcoal. Bierman (1994, p. 13,893) mentioned in a review, “In every study to date, exposure age variability on a single geomorphic surface is significant (usually >10%) and in most cases greater than expected from analytical uncertainty.” Furthermore, most of these early studies were conducted on eroding surfaces; for rock clasts on an offset (inactive) alluvial fan, the inheritance includes some fraction of the transport time as well as the time needed for erosion. It is often thought that boulders are the best sampling targets on fans because their transit times were the least (e.g., Blumentritt et al., 2005a). However two methodological studies by Blumentritt et al. (2005b) and Perg et al. (2005) of alluvial fans in the Mojave Desert of California established that ages of desert pavement are younger and more reproducible than ages of boulders; but they still show inheritance, which can be estimated by sampling to 2–3 m depths in pits (Machette et al., 2003). The correction for inheritance in arid Death Valley is typically 20–130 ka (Machette et al., 2003), but on the forested northern side of the San Gabriel Mountains it is only 1–20 ka for boulders and 3–8 ka for fine sediment (Matmon et al., 2005). In cases of cosmogenic nuclide ages that were not already corrected for inheritance, I use equations 22–24 to broaden the age range toward the young side, taking ω = 3.5 × 10−5 a−1, corresponding to a median inheritance of 20 ka Unfortunately, this is a very crude allowance for an effect that probably depends strongly on basement rock type, topographic slope in the eroding region, slope direction, and the many other variables that can be loosely termed “climate.”

Stated Age is Bounded by Two Limits

Many offset features are not datable at all, but have ages constrained by stratigraphic succession and/or crosscutting relations with other beds that have been dated. For example, an offset feature in rocks containing Miocene fossils is younger than the Oligocene-Miocene faunal transition at 23.8 Ma, but older than the Miocene-Pliocene faunal transition at 5.3 Ma. These bounding ages are each subject to their own uncertainties (as shown by past revisions of the geologic time scale; Harland et al., 1990). The appropriate pdf for such cases is a uniform distribution with smoothing of the left and right discontinuities to reflect the uncertainty in the upper and lower limits:  
formula

The comments listed after equation 21 about the selection of standard deviations also apply here.

Stated Age is a Lower Limit

Tabulated data also include some cases in which only a lower limit is provided for the age of the offset feature. Examples include uranium-series (closed system) and uranium-trend (open system) dating of soil carbonates, including carbonate rinds on clasts in alluvial fans; the events being dated postdate sediment deposition by an unknown amount of time.

A pdf cannot be defined for a random variable that is potentially infinite. The age of the Earth's oldest surviving crust (ca. 4 Ga) can always be used as Δtmax in equation 25, but this only solves the problem in a formal sense. In California, I assume a smaller default maximum age of Δtmax± σΔtmax = 135 ± 70 Ma, based on the Jurassic–Cretaceous age of the prominent plutons, volcanics, and sedimentary wedges that make up a large fraction of the upper crust. Even this more modest age limit can lead to artifacts. For example, if it is known that the San Andreas fault slipped 750 m in >0.03 m.y. (Rasmussen, 1982), then the obvious conclusion is that the slip rate is <25 mm/yr. But if we assume a uniform probability distribution in 0.03–135 m.y. for the age of the offset feature, then the 95% confidence limits on the age become 3.4–131.6 Ma, and the 95% confidence limits on slip rate become 0.0057–0.22 mm/yr, values that are clearly too low. (In cases where the investigation is still active, a tighter constraint should be supplied by reviewing the regional geology and entering the age of the oldest formations exposed locally as the upper limit on age of offset.)

To avoid this artifact, I make an additional assumption that any stated minimum age also provides an order-of-magnitude estimate for the true age, which is very unlikely to be more than 10 times larger. Therefore, I use an additional exponential term to weight the probability density near the lower limit  
formula
where λ = –ln(&12frac;)/Δtmin This assumption gives a more reasonable result (i.e., a distribution of L extending from ∼15% of the expected upper limit, to the expected upper limit) in the processing of most published results. (However, it will lead to a new artifact of rate overestimation if offsets are published with minimum ages that are 10 times smaller than the true age. Such data have little information value, and should normally be screened out during the process of peer review.)

Stated Age is an Upper Limit

When only an upper limit on the age of the offset feature is available, the pdf is defined as a uniform distribution with one rounded shoulder, by simply dropping the erf() factor from equation 25:  
formula

LONG-TERM OFFSET RATE, L, FROM A SINGLE OFFSET FEATURE

Combining equation 2 with equation 6 or equation 7, we get two alternative convolutions for the probability density function (pL) of the long-term offset rate (L) of a particular offset feature on a particular fault:  
formula
and  
formula
where ΔD, Δt, and L are all non-negative. In most cases, equation 29 will be used for numerical computation, because it does not involve an integrand that is singular as ℓ' → 0.

One property of these pdfs for long-term rate that may be unexpected is that they are not symmetrical, even when the pdfs for far-field offset and for age are each symmetrical. This is a result of the division of the former by the latter, and the resulting partial-derivative terms (cf. equations 6 and 7) that appear in equations 28 and 29. The probability densities of rates near the lower limit are higher than the probability densities of rates near the upper limit, so that the four central measures of the rate distribution typically occur in the order: (mode) < (median) < (mean) < (mean of nominal lower and upper limits). An example is shown in Figure 4. Previous studies (including my own) that did not consider this effect probably had a subtle systematic bias toward artificially high long-term rates.

One necessary caution is that these equations only define distributions for the average offset rate since creation of a particular offset feature; there is no implication that the actual rates have always been equal to the average rate. Possible time dependence of rates is assessed later in this paper, by comparisons between different offset features.

FAULT TRAINS VERSUS FAULT SEGMENTS AND FAULT SECTIONS

In the following I contrast and then combine long-term offset rates from different offset features along the same fault: initially to learn about the prevalence of errors and rate changes, and later to determine the best possible combined rate. Such operations require deciding which rates are comparable, and the logic to follow requires that we select them on an a priori basis (using only the locations of the offset features relative to the fault network), not on an a posteriori basis (considering the rate values). So, it is useful to have a word or phrase to designate the geographic extent within which rates can be compared, and some rules for defining these geographic extents.

I use the phrase fault train to mean a contiguous piece of the trace of a fault system along which our knowledge of fault geometry permits the null hypothesis of uniform long-term offset rate. One basic rule for defining fault trains is that each one must end at a triple junction (or higher-order junction) of active faults, because the addition and/or subtraction of slip-rate vectors that occurs at a junction normally implies a change in all components of long-term offset rate. There may also be additional points along a trace that require a boundary between adjacent fault trains, but this depends on the particular offset rate in question, and on fault geometry. For example, a right-angle bend in a fault trace where a strike-slip fault (tear fault) connects to a dip-slip fault would be a point of expected change in all the offset measures typically catalogued (strike-slip rate, throw rate, and opening and/or closing rate), even though it would not be a point of change in the heave rate vector. Thus it should be a boundary between two fault trains unless the offset rate in question is the heave rate vector (a rare case). On the other hand, some dip-slip faults have complex traces with many bends, where there are also changes in dip (lateral ramps, mullion), but have parallel slip vectors of uniform trend and plunge. If the offset rate in question is the throw rate, such complexities of the trace would not require subdivision of the trace into different fault trains.

The subdivision of fault traces into trains should ideally be independent of seismicity, because long-term offset rates have already been defined as far-field rates, or alternatively as averages over many earthquakes. In contrast, the fault segments mapped by previous Working Group(s) on (Northern/Southern) California Earthquake Probability were recognized primarily on the basis of the extent of historic or paleo-seismic ruptures, and/or expert opinion about the likely extents of future ruptures. Because fault intersections and fault bends can be (or at least are often thought to be) barriers to rupture, there are many fault segment boundaries that are also mandatory fault train boundaries, but the two concepts are still distinct. Some segment boundaries based on terminations of past ruptures are not required to be fault train boundaries, because there is no geometric basis for expecting a change in long-term offset rate at those points.

The newer fault sections mapped by the most recent Working Group on California Earthquake Probabilities are defined (in draft documents by C.J. Wills and coauthors) as “parts of faults that are defined by changes in geometry, slip rate, geomorphic expression, or earthquake rupture history.” Therefore, the set of fault section boundaries is the union of a minimal set of required fault train boundaries with the set of fault segment boundaries. So, fault sections are typically shorter and more numerous than either fault trains or fault segments. Any fault section is a valid fault train because it fits the definition above. However, if the fault section was limited by rupture extent, or by a boundary in geomorphic expression, it may be shorter than the longest possibly fault train that could be defined in that area.

Defining fault trains that are shorter than the longest possible fault train in the area is not an error in logic. Any valid fault train can be subdivided into shorter fault trains (which will have identical long-term offset rates). However, it is counterproductive, and therefore an error in scientific strategy. It adds parameters to the regional model that are not needed. At the same time, it reduces the number of pairs and multiples of offset features that can be contrasted to learn more about sources of error. It also reduces the amount of data redundancy in the computation of combined rates, thereby increasing their uncertainties.

DISAGREEMENT BETWEEN TWO LONG-TERM OFFSET RATES

According to the conceptual model of rigid microplate tectonics, offset rates should typically be nearly constant along any fault train. However, slip rates might still vary if the distance to the (typically unknown) Euler pole were less than several times the width of the smaller microplate. Other tectonic paradigms with more distributed deformation away from (mapped) faults allow more parameters for offset rate variation along the trace. In addition, if one or both of the offset features is many millions of years old, the long-term offset rate for that fault may have changed through geologic time, causing a disagreement between rates averaged across different time windows. Long-term offset rates may also disagree in cases of fundamentally incorrect data, including failure of the piercing-point pair to span the whole fault zone, undetected complex initial shapes of offset features, or faulty correlations in space or in geologic time. Second- or third-hand data compilations may introduce more artificial disagreements through clerical error, overgeneralization, or oversimplification of the uncertainties.

In order to assess the importance of each of these effects, it is useful to have a quantitative measure of the level of disagreement between any two long-term offset rates (for the same component of slip, on the same fault train). I concentrate first on disagreements between pairs of features, because this is conceptually and mathematically simpler than an ensemble analysis. However, once I have estimated the fractions of incorrect and inapplicable data, I use these in ensemble formulas to estimate combined distributions of offset rates as accurately as possible.

Because the long tails of most individual pL(ℓ') distributions will almost always overlap to some extent, agreement or disagreement is not well represented by a purely binary classification. Instead, a probabilistic scalar measure of disagreement is appropriate.

One condition sufficient to demonstrate that two independent scalar random variables X and Z have different values is that their values are on opposite sides of some intermediate value y: X < yZ. Another sufficient condition is the mutually exclusive alternative Z < yX. Here I adopt as a disagreement measure the scalar quantity δ defined by:  
formula
in terms of the two cumulative probability functions (for X and Z both non-negative).

A test of two distributions that have no overlap, and disagree completely, gives δ = 1. When X and Z have Gaussian distributions that are identical, δ = &12frac;, which is the lowest possible result. The same result of δ = &12frac; is obtained when X and Z have Gaussian distributions with the same median but different standard deviations. This is an important property of the δ measure, which is not shared by some other popular metrics of differences between distributions (e.g., the maximum of the absolute value of the difference between the two cumulative probability functions). Use of δ is appropriate in Bayesian models like this project, where the distribution functions for long-term offset rates portray the plausibility of alternative rates given available information, rather than the frequency with which rates appear during repeated sampling. That is, δ is a conservative measure of differences between distributions, which tends to minimize the difference in cases of great uncertainty.

Whenever we compare two rate data from the same fault train, our null hypothesis is that the rates are the same. Low values of δ are consistent with the null hypothesis, but high values of δ indicate that this null hypothesis can be rejected with a high degree of confidence. To quantify this relationship, it is necessary to consider δ as a random variable, with randomness inherited from the measurement errors in the two data that were then expanded into model cumulative-probability functions

\(\mathit{P}_{\mathit{L}}^{X}\)
and
\(\mathit{P}_{\mathit{L}}^{Z}\)
and compared. I conducted Monte Carlo experiments on many pairs of simulated rate data sampled from a common Gaussian distribution. The model rate distribution for each of these simulated data was also Gaussian, but was centered on the sampled value rather than the true value. The standard deviation of each model distribution was the same as the standard deviation of the common parent distribution. The result was that values of δ > 0.846 occurred in only 5% of the cases: P(δ < 0.846) = 0.95, under all the conditions stated above, when the null hypothesis is true. However, the case of two equally broad distributions is not typical. In other Monte Carlo tests where the two simulated data were sampled from two Gaussian distributions with very different widths (but still sharing a common median, representing the true value), the result was the threshold δ ' for which P(δ < δ') = 0.95 rose to approach, but not equal, 0.95. Therefore, I propose the following.

Conjecture

When two data are samples of a common value, but measurement methods have introduced random errors, and when the model distribution functions for those two data accurately reflect the distributions of sizes of the measurement errors, then P(δ< 0.95) ≥ 0.95. That is, there is no more than 5% chance of observing δ > 0.95 when comparing two honest but imprecise measurements of a common value.

Once δ values have been computed for all the comparable pairs of long-term offset rates in a large data set, it is possible to begin to infer the main sources of disagreement. For example, the importance of changing slip rates over geologic time can be estimated by plotting δ as a function of the absolute value of the time gap between the ages of the two offset features. The importance of relative microplate rotation and distributed deformation might be estimated by plotting δ as a function of separation distance along the fault trace. The importance of clerical and oversimplification errors can be assessed by comparing δ values for pairs of close, young offset features in primary versus secondary or tertiary literature. It is probably reasonable to lump together δ values from all the different classes of offset measures (dextral and sinistral strike slip, normal, and thrust fault throw) when conducting these empirical studies. Such lumping makes it possible to collect enough comparisons for statistical inferences.

ESTIMATION OF THE FRACTIONS OF INCORRECT AND INAPPLICABLE DATA

Once physical causes (variations of offset rate in time and space) and clerical causes (second-hand reporting) for disagreements have been minimized, there is likely to be a residual population of disagreements that is largely due to fundamentally incorrect interpretation of the geologic history and/or miscorrelation of features or strata. This fraction of incorrect data can be estimated from the pattern of disagreements. Let us adopt δ ≥ 0.95 as the criterion for a strong disagreement, corresponding to at least 95% confidence that the two long-term offset rates are different. Then, <5% of comparisons between correct data (for the same component of slip, on the same fault train, in the same tectonic epoch, from primary sources) should lead to a strong disagreement. Also assume that incorrect data are sufficiently erratic in their values so that comparisons of incorrect data with correct data (or with other incorrect data) will always lead to strong disagreement. This assumption amounts to a restrictive or conservative definition of incorrect, which will, in turn, tend to minimize the fraction of incorrect data that we infer at the end of this study. If α is the fraction of incorrect data, then the fraction of comparisons that result in strong disagreements (β) should be  
formula
which is equivalent to the inverse relation  
formula
so that (for example) β = 0.1 → α ≥ 0.027; or β = 0.3 → α ≥ 0.14; or β = 0.5 → α ≥ 0.27.
The value of α obtained from such comparisons (equation 32) between data regarding young offset features (presumed to be from the same tectonic epoch) estimates the fraction of the data that are fundamentally incorrect; this may be labeled
\({\alpha}_{neotectonic}^{primary}\)
. However, if comparisons are extended further back in geologic time, we expect βprimary(t) and αprimary(t) and to increase, because in addition to incorrect data there will now be some data that are inapplicable to neotectonics. That is, even correct rates will begin to disagree because some rates have varied across different tectonic epochs. The increase in αprimary(t) from its neotectonic value provides an estimate of the fraction of inapplicable data. However, now the situation is less symmetric than before, so equations 31 and 32 no longer apply. If βprimary(t) increases when comparing rates from progressively older offset features with those from young offset features, this increase in β is due almost entirely in increased probability of inapplicability of the rate based on the older offset feature. Therefore,  
formula

For purposes of combining rate data into a single best-estimate neotectonic rate, it makes little difference whether a rate is misleading because it is incorrect, or because it is inapplicable. Both possibilities will be expressed by the single fraction α(t) in the following section.

MERGING MULTIPLE RATE DATA INTO ONE COMBINED PROBABILITY DENSITY FUNCTION

For many projects, it may be necessary or desirable to merge all available long-term offset rates (for the same component of slip, on a single fault train), producing one combined

\(\mathit{p}_{\mathit{L}}^{combined}\)
(ℓ′) For example, this step is required before simulating neotectonics with our kinematic finite-element program NeoKinema (Bird and Liu, 2007). Taking this step implies a working hypothesis that variations in long-term offset rate along one fault train (due to nearby Euler poles and/or distributed anelastic deformation) are typically minor in comparison to the apparent variations caused by measurement errors and/or real variations caused by time-dependence. Those who do not accept this working hypothesis should simply use rates from individual offset features as point estimates, and ignore the combined pdfs of fault trains derived here.

We cannot expect to accurately treat those faults whose long-term offset rates vary along their traces with a simple algorithm like the one developed in the following. Nevertheless, some attention to possible artifacts of a uniform offset-rate assumption can prevent egregiously misleading results. The most serious cases are those with at least one datum showing an older overlapping formation (indicating that the fault has been locked since that time), combined with at least one datum showing younger scarps or other positive offsets. Examples include the Hoback normal fault of Wyoming, the Amargosa–Death Valley detachment fault of California, and the San Gabriel–Vasquez Creek dextral fault of California. If we consider such logical contradictions as necessarily proving error in at least one of the data, then we can combine them with a standard procedure (detailed below) that takes account of the α of each datum. Experience shows that the fault-pinning overlap data will often “win” in the sense of producing a combined pdf for long-term offset rate that has a very small median and very small upper 95% confidence limit. However, a more plausible interpretation is that this particular fault had variable long-term offset rate along its trace, and that the pinning overlap formation(s) should be regarded as effectively reducing the length of the fault trace that is currently active. In the absence of information on the probabilities of slip rate variations (along the trace) of various magnitudes, it is not yet possible to weight these alternatives quantitatively. However, I would not want to communicate invalid assurances about the safety of possibly dangerous faults. Therefore, I adopt a somewhat arbitrary rule here: fault-pinning (zero rate) data of greater age will not be used in computing combined offset rates when there is at least one datum of younger age showing positive offset. Instead, the resulting combined pdf will be marked with an asterisk to indicate that the active fault length in the neotectonic era is less than what it was in some paleotectonic era, and that the fault-pinning data should be considered when defining a reduced active fault length.

With these caveats, I now merge multiple rates (for a single offset component on a single fault train) under a working hypothesis that the long-term offset rate is a single non-negative number. The merging of multiple rates should take into account the probability (αi) that any individual rate is incorrect or inapplicable. Even with a number (i = 1,2,...n) of rates, there is also a non-vanishing probability (graphic) that all the data are incorrect or inapplicable. Therefore, it is necessary to invoke a prior distribution to write a formula for the combined pdf. Fortunately, these prior distributions can be obtained by iterative bootstrap processing, as described below. The prior distribution also serves as estimate in cases of faults that have no associated data.

When there are no data concerning a particular fault train,  
formula
When there is one datum concerning a particular fault train,  
formula
Where there are two data concerning a particular fault train,  
formula
in which the superscripts on the probability density functions are identifying indices (not exponents), and Lscale ≡ 1 km/m.y. is introduced to correct the units of the product before normalization. In numerical computation this factor has no effect and may be omitted.
Where there are three data concerning a particular fault train,  
formula
The pattern is clear: the number of terms increases as 2n, while the computational work increases as n2n. This makes it impractical to write (or compute) the exact result for cases with n = 12–18, which occur along some trains of the San Andreas fault. I handle such cases with an approximate, iterated formula:  
formula
for all i (sequentially increasing) in the range 2 ≤ in, where
\(\mathit{p}_{\mathit{L}}^{uniform}\)
is defined as:  
formula
and the range in which all distributions are computed is 0 ≤ ℓ'< Lmax Approximate formula 38 preserves the same number of terms, and the same set of leading coefficients [permuted serial products of αi or (1–αi)] as exact formulas like equations 36 and 37. However, it is approximate for two reasons. (1) The first term of the expanded polynomial form of the solution (which is generally the smallest) is proportional to graphic rather than to graphic and (2) weights added by the normalizing operator N() are somewhat different because the operator is applied sequentially to each factor, instead of being applied once in each term of the result (except the leading term). This formula also has important advantages: (1) closed form for any value of n, allowing implementation in a program; (2) computational work just proportional to n; and (3) simple interpretation. Equation (38) can be seen as a serial string of identical signal processors (one per datum), each of which is a leaky filter. The leak term passes the previous signal through unchanged except for fractional gain of αi; this represents the case that datum i is incorrect or inapplicable, and therefore should have no effect on the result. The filter term multiplies the previous pdf by the pdf of datum i and applies gain (1 –αi); this represents the case that datum i is correct and applicable.

Tests were run on sets (n = 4–18) of real data, after the data were sorted either in order of increasing nominal rate, or in order of decreasing nominal rate. Computed medians and standard deviations of the combined rates were the same for either order of data, confirming the desired symmetry of this approximate algorithm, even in the presence of numerical roundoff error.

In another test, accurate equation 37 was compared to approximate equation 38 for a number of actual cases of faults with n = 3 data (even though equation 38 is normally reserved for cases with n ≥ 4). Figure 5 shows one of these comparisons in detail. Differences in the inferred median rate were never >10%. Approximate equation 38 tends to estimate a lower standard deviation for the combined rate in those difficult cases where all the αs are large. This effect may be due to the poor representation of the first term of equation 37 by equation 38; however, the effect is expected to become less important as n increases.

When the number of offset features is large, there are typically some strong disagreements among them. In such cases, the combined pdf computed by these methods may be multimodal. An example is shown in Figure 6.

NUMERICAL COMPUTATION

All equations in this paper have been incorporated into a Fortran 90 program. The source code is attached as file Slippery.f90; an executable version for Windows is attached as file Slippery.exe1. Each run through the program performs batch process on an entire offset-rate database; therefore, the iterative or bootstrap improvement of the prior distributions of offset rates

\(\mathit{p}_{\mathit{L}}^{prior}\)
, by employing the latest posterior distributions, requires only about three runs. Each pdf and cumulative distribution is represented internally by 4000 discrete values between nominal limits on the independent variable (xminXxmax) that are estimated (in advance) to correspond roughly to P(Xxmin)≅ 0.0001 and P(Xxmax) ≅ 0.9999. During normalization of any distribution, these bounds are treated as if they were absolute (0 ≤ P ≤ 1). When X is uniformly positive, and the ratio xmax/xmin exceeds 10, discrete values are evenly spaced on a logarithmic scale. When the range of X extends to zero (some ages and some fault offsets) or to negative numbers (ΔE), the discrete values are evenly spaced on a linear scale. Integrals are simply estimated as sums of 4000 mid-point values multiplied by their respective Δx's. All computations are done with eight-byte real numbers.

Some results will include nonzero probability densities for extremely large offset rates, which contradict our basic assumptions. This is especially true when the age of the offset feature has only an upper limit (equation 27). The same problem can arise in cases using the other types of pdf (equations 21–25) for the age of the offset feature, if laboratory uncertainties or inheritances are a large fraction of the nominal age. While unreasonable offset rates could be prevented by introducing a Bayesian prior model into all calculations, this would interfere with the planned bootstrap study of levels of agreement or disagreement between rates. In this part of the project, I take a simpler approach: I arbitrarily truncate the computation of pL at an offset rate equal to four times the relative plate velocity in California (228 mm/yr), and normalize only the part of the pdf below that limit. The very slight loss of information caused by this truncation will be handled by reporting extremely large or unbounded rates with a special code (U for unreasonable or unbounded) in the tables, instead of reporting very large numbers that might cause alarm if taken out of context.

APPLICATION TO OFFSET FEATURES IN THE WESTERN UNITED STATES

Since beginning the research that led to a preliminary paper (Bird and Rosenstock, 1984), I have kept a personal database of published offsets along faults in western North America, and their ages. Because this database also supports paleotectonic studies (e.g., Bird, 1998), it is not limited to young offset features, and it contains some faults that are considered inactive. Each item in the database has fields for fault train index (required), sense of offset (required), citation (required), an offset distance (or range, or upper or lower limit) in kilometers (if available), a standard deviation for the offset (if available), the geologic time (or range of times, or limit on the time) (in Ma) when this offset began to occur (if available), the geologic time (or range of times, or limit on the time) (in Ma) when offset had ended (if available), the inferred offset rate (or range, or limit on the rate) in km/m.y. = mm/yr, and the standard deviation of this inferred offset rate. I estimated these inferred rates and standard deviations subjectively and entered them manually. They will not be used in this project, which is intended to replace them. The fault train index identifies a particular named and digitized trace of a fault train (defined above) with a four-digit NeoKinema index number, such as “F0402 = San Andreas dextral fault, train G (Salton), CA.” The sense of offset is represented by a single-letter code: D = divergent heave across a low-angle normal fault; L = left-lateral strike slip; N = throw of high-angle normal fault; P = convergent heave across low-angle thrust plate; R = right-lateral strike slip; T = throw of thrust fault and/or associated monocline. This database structure means that there are no negative entries. This database has some serious deficiencies: (1) many offsets that were published in the past decade have not yet been entered, due to imperfect upkeep; (2) many values are qualified by equivocation characters (e.g., “∼5” or “5?”) that do not follow defined rules; (3) many optional fields are blank; (4) latitude and longitude coordinates for offset features were not recorded; (5) the type of geochronology employed was not recorded; (6) some of the literature sources cited (e.g., abstracts) do not meet the high standard discussed previously as an ideal.

Because one person's part-time effort is not sufficient to track each piece of information on an offset feature to its original source, this database also includes offsets from secondary and tertiary sources in the literature. Fortunately, these were systematically labeled. Primary sources (A) are those that give the original report of an offset by the investigators who did the mapping or trenching, or determined the age of an offset feature using laboratory geochronology. Secondary sources (B) are those that summarize the literature on fault activity in a limited region, such as around one fault, one mountain range, or one basin. Tertiary sources (C) are those that summarize the literature for an entire state or larger region (e.g., Basin and Range province). For some minor faults in unpopulated regions, a tertiary reference may be all that is available. For example, Stewart (1998) noted maximum initiation ages (varying by region) for the normal faults he compiled in the Basin and Range province; in some regions he also noted that scarp heights are generally 0.3–2 km, and that this is a lower limit on the true throw; in other cases he described the total throw as 2–5 km. Such constraints have been entered as C-class data under each individual trace that appeared in the compilation maps of Stewart (1998).

Before processing with these new algorithms, the database was lightly edited for geographic extent and for geologic time window. I kept all faults with known or potential Neogene activity in the conterminous western United States, and in those parts of the borderland underlain by continental crust. I removed most faults that are entirely in Alaska, Canada, or Mexico, and entries concerning an earlier phase of slip that probably had a sense or rate different from the Neogene rate (e.g., Howell's 1975 discussion of a proto–San Andreas fault system). Thus, the edited database corresponds roughly to the continental portion of the Gorda-California-Nevada orogen (Bird, 2003), with geographic extensions eastward to the Rio Grande rift region. It contains 849 faults, with 3577 entries. This version of this database is seen in Table 1 (files Table_1.xls or Table_1.txt; see footnote 1), which includes both input to the statistical analysis and a summary of its output. A partial version including only the first nine columns, which are read as input data by program Slippery, is the tab-delimited ASCII file f_Gorda-Cal-Nev.txt, which is attached (see footnote 1). This latter file can be used in conjunction with Slippery.exe to reproduce the analysis.

However, the algorithm presented here does not use any inferred offset rate or its associated standard deviation, but only offset distances with associated offset age where both are recorded in the same entry of the database. A human editor might combine one author's estimate of total offset with another author's estimate of the age of fault initiation, but this algorithm will not attempt to do so. With this more strict criterion, there are only 995 data concerning dated offsets on 464 distinct fault trains. The breakdown of data by offset type is: D 4%, L 3%, N 61%, P 1%, R 22%, and T 8%. The number of distinct literature citations is 221. The breakdown of data entries (not references cited) by literature type is: 24% primary (A), 13% secondary (B), and 63% tertiary (C). A complete list of the full bibliographic citations corresponding to the short citations in Table 1 is attached as file references_cited_Table_1.doc (see footnote 1).

Because this database has no field for the type of geochronology employed, I arbitrarily assumed that all ages younger than 30 ka were 14C ages, and all other ages younger than 600 ka were cosmogenic nuclide ages, and in both cases I corrected for potential inheritance with equations 23 and 24. In those cases where this correction was inappropriate or redundant, the uncertainty in the age has been artificially inflated, reducing the precision of the offset rate. This artifact could occasionally obscure real disagreements between rates.

To begin the analysis, uniform distributions were used as models of the six prior offset-rate pdfs

\(\mathit{p}_{\mathit{L}}^{prior}\)
' for offset types D, L, N, P, R, T. Their assigned ranges were:
\(\mathit{L}_{\mathit{D}}^{prior}\)
≥ mm/yr;
\(\mathit{L}_{\mathit{L}}^{prior}\)
≥ mm/yr;
\(\mathit{L}_{\mathit{N}}^{prior}\)
≥ mm/yr;
\(\mathit{L}_{\mathit{P}}^{prior}\)
≥ mm/yr;
\(\mathit{L}_{\mathit{R}}^{prior}\)
≥ mm/yr; and
\(\mathit{L}_{\mathit{T}}^{prior}\)
≥ mm/yr. As each datum-based pdf pL(') was computed, it was averaged into one of six empirical ?posterior offset-rate distributions (except that unbounded rates were not used in this step). In the second pass of this processing, these six posterior distributions were used as prior distributions, and so on. This iterative estimation of the datum-based prior distributions was well converged by the third run.

Inferring probabilities of erroneous or inapplicable data, and computing combined long-term offset rates for each fault section, are both operations that require each offset datum to be associated with a fault train. The traces of the fault trains used in this phase of the work are contained in file f_GCN_Bird2007_Table1.dig, which is attached to this paper (see footnote 1). The links between Table 1 and this file are the NeoKinema fault train numbers (e.g., F0402; described previously). These trains were mostly generalized and hand-digitized from more detailed sources, such as Jennings (1994) or the Quaternary Fault and Fold Database maintained online by the U.S. Geological Survey. A decision was made at the start of the compilation to simplify multiple ruptures and/or discontinuous traces seen in Quaternary surface units (whose strain is often small) into a single continuous trace representing the estimated fault position at seismogenic depths (where strain is often larger). Thus the entire San Jacinto–Coyote Creek dextral fault system of California is represented by a single train (F0544). This practice has the effect of minimizing the number of fault trains and maximizing the opportunities for comparing data, and for merging them into combined rates. However, I consistently divided long fault traces into trains at points where their traces intersect with those of other major faults, because it is expected that long-term offset rates change at such junctions. Thus, the San Andreas fault is represented in this compilation by seven trains.

The most interesting results concern the 323 cases in which rates from primary sources could be compared to rates from other primary sources (A/A) describing the same offset component on the same fault train. As defined above, δ ≥ 0.95 was taken as the definition of strong disagreement, and the fraction of comparisons resulting in strong disagreement was described by symbol β primary(t). I further sorted these comparisons by the difference in age between the two offset features, and divided them into seven groups of equal population along the age-discrepancy axis. This yielded seven crude estimates for β primary(t) that are presented in Figure 7. The result was that β primary(t) shows an approximately constant value of ∼13% for age differences up to ∼3 Ma, after which it abruptly increases to ∼56%.

This behavior is understandable given what is known of geologic history in the region. During the latest Cretaceous and early Tertiary (80–30 Ma) there was low-angle or horizontal subduction starting at the Pacific margin (Bird, 1998; 2002). Slab roll-back ca. 30 Ma caused elevations to increase and stress directions to shift dramatically (Bird, 2002). Tectonics changed further during the interval 28–20 Ma as new Pacific–North America transform plate boundary segments were created near the margin. A much more rapid change began ca. 20 Ma, when a subducting microplate froze onto the Pacific plate and abruptly changed its direction of motion (Nicholson et al., 1994), tearing open the Los Angeles basin and rapidly rotating the future Transverse Ranges. The last major transition occurred 10–5 Ma, when the Gulf of California began to open, the southeastern sections of the San Andreas fault were formed with a large initial left step, and the Transverse Ranges began to be elevated by thrusting. I speculate that resistance to Transverse Ranges formation may have accelerated the Sierra Nevada–Great Valley block to faster northwest motion (with respect to stable North America), thereby rotating stress fields and affecting tectonics throughout the Basin and Range province and perhaps even to the Rio Grande rift (Bird, 2002). Thus, the only two periods when we might expect relatively stability of rates are 85–30 Ma (from which only a very small number of offset features have been recorded) and 5–0 Ma (which provides the majority of the recorded offset features).

Although this history suggests that it might be safe to average rates from offset features as old as 5 Ma in neotectonic studies, my preference is to rely on an empirical model in which the risk of inapplicability begins to rise for offset features older than 3 Ma (Fig. 7). This model beta curve is also shown in Figure 7. To construct this model, I lumped all disagreements between primary offset rates since ≤3 Ma into the single value of

\({\beta}_{neotectonic}^{primary}\)
= 24/185=0.13 Then, according to equation 32,
\({\beta}_{neotectonic}^{primary}\)
≤ 0.043. That is, at least 4.3% of primary reports of neotectonic offset rates are incorrect or unrepresentative. For offset features older than 20 Ma, equation 33 implies αprimary(t) ≅ 0.522, which can be broken down as ∼5% of primary data that are incorrect or unrepresentative, and ∼47% of primary data that are correct but inapplicable to neotectonics.

There were also 397 opportunities to compare rates based on two tertiary literature sources (for the same offset component and fault train). Of these, 206 concern offset features whose ages differ by ≤3 m.y., giving

\({\beta}_{neotectonic}^{tertiary}\)
= 63/206=0.306. According to equation 32, this implies
\({\beta}_{neotectonic}^{tertiary}\)
. ≥ 0.145. This is about three times the error rate of primary literature sources. Note that clerical error is not the only cause of additional error; oversimplification of uncertainties and overgeneralization about rates of faults in some perceived (but possibly imaginary) class are probably also important.

Only 33 comparisons between offset rates based on secondary literature sources are available, and only 10 of these concern pairs of offset features with ≤3 m.y. difference in age. While the empirical

\({\beta}_{neotectonic}^{secondary}\)
= 1/10=0.1 is intermediate between those for primary sources and for tertiary sources, the number of comparisons is too small to rely on this as a good empirical estimate of
\({\alpha}_{neotectonic}^{secondary}\)
. For purposes of further computation, I assume instead that
\({\alpha}_{secondary}^{neotectonic}\)
≅(
\({\alpha}_{primary}^{neotectonic}\)
+
\({\alpha}_{tertiary}^{neotectonic}\)
)/2 = 0.094.

Finally, these empirical values of α (as a function of literature source type and age of the offset feature) were used in equations 34–39 to compute best-estimate combined pdfs for long-term offset rate of each the fault trains in my database. Table 1 presents 7 defining measures of each distribution: lower bound [of 95% confidence range, at P(L < ℓ) = 0.025], mode, median, mean, upper bound [of 95% confidence range, at P(L < ℓ) = 0.975], and standard deviation. Many distributions are skewed so as to have longer higher probability densities at low offset rates compared to high offset rates; therefore, the ordering of the three central measures is typically: (mode) < (median) < (mean). My preference is to choose the median when only a single numerical estimate of the long-term offset rate is requested, but also to note the formal standard deviation.

APPLICATION TO OFFSET FEATURES IN CALIFORNIA

I also used the same program (Slippery) to conduct a second, largely independent analysis of the PaleoSites Database subset of the U.S. Geological Survey Quaternary Fault and Fold Database. The PaleoSites Database was developed by Susan Perry with assistance from Edward Field, Chris Wills, and Vipin Gupta. It is intended to summarize all published data on offset features and/or paleoseismicity on active and potentially active sections of faults in California, supplemented by the interpretations (slip rates, numbers of paleoearthquakes) of the original authors.

Specifically, I processed the combined-events portion of the Paleo-Sites information as it appeared in spreadsheet form as of 14 February 2007. There were 471 rows referring to the same number of paleosites. However, as WG_FAULT_ID (Working Group on California Earthquake Probabilities fault section identification number) and FAULT_NAME are primary organizing fields, some paleosites are necessarily double entered under two alternate fault names (with different fault identification numbers). If entries in which FAULT_NAME includes “alt 2” are not counted, then the number of rows (and number of paleosites) is 421. Of these, 158 originated in the older Fault Activity Database of the Southern California Earthquake Center (SCEC-FAD); 87 originated in the detailed reports of the online U.S. Geological Survey Quaternary Fault and Fold Database (QFaults); and 176 originated in my personal database, having been selected for PaleoSites by Wills and/or Perry. The advantage of merging these sources is that QFaults and SCEC-FAD are more up to date than my compilation, but my compilation extends further back in geologic time.

Not all of these entries could be processed by my program. It is necessary to have entries in the SENSE_OF_MOTION field, in the MEA-SURED_COMPONENT_OF_SLIP field, in at least one of the start-time fields [PREF_START_TIME, MAX_START_TIME_(earliest), and/or MIN_START_TIME], and in at least one of the offset-distance fields [PREF_OFFSET_(m), MAX_OFFSET, and/or MIN_OFFSET]. This limited the usable rows and/or paleosites to 239 (not counting duplicate alternate-trace entries, and also only single counting very large offsets that are multiply entered under several fault sections along the same fault). Of these, 2 rows and/or paleosites were from SCEC-FAD, 64 were from QFaults, and 173 were from my personal database. Therefore, there is very substantial redundancy between the information presented in this section and that analyzed in the previous section. The principal differences are the limitation of the geographic scope to California, and the more up to date status of PaleoSites.

My method of analysis requires codes for sources in primary literature (A), secondary literature (B), or tertiary literature (C), as defined in the previous section. For the PaleoSites entries originating from my database, this code was extracted from field BIRD_FAULT_ID-REFCAT. All other entries were considered to be from primary literature (class A), as this is the practice of the online QFaults. The result was 211 class A rows, 51 class B, and 2 class C.

As in the previous section, analysis was begun with uniform distributions assumed for each of the six prior long-term offset-rate pdfs for the six types of offset (D, L, N, P, R, T). The analysis was repeated three times, with iterative replacement of the prior distributions by the posterior distributions, with good convergence.

The input and output of this analysis are both displayed in Table 2 (files Table_2.xls or Table_2.txt; see footnote 1).

As before, I compared all possible pairs of single-feature long-term offset rates for the same component of slip on the same fault section and computed the disagreement measure δ. Note that the division of California fault traces into sections in PaleoSites follows Fault Section Database v.2 (created by C.J. Wills and coworkers in the Working Group on California Earthquake Probabilities), which was based on Cao et al. (2003), which was based on Petersen et al. (1996). This is different from my division of fault traces into a minimal set of fault trains discussed in the previous section. However, any fault section should be a valid fault train.

There were not enough comparisons within literature classes B or C to draw any conclusions, but there were 447 A/A comparisons, including 304 in which the difference in ages of the two offset features was no more than 3 m.y. This population yields

\({\beta}_{neotectonic}^{primary}\)
= 41/304=0.135 (almost identical to the previous result for the western conterminous U.S. as a whole), which is not surprising as a large fraction of the offset features were also in the group previously analyzed. Then, according to equation 32,
\({\beta}_{neotectonic}^{primary}\)
≥ 0.046. Based on either database, it appears that ∼5% of reports of young (≤3 Ma) offset features in the primary literature are incorrect or unrepresentative.

The only notable difference in overall statistical results is seen in Figure 8, where subsample β(t) values are plotted against the time-dependent β(t) model that was used in the previous analysis (copied from Fig. 7). Instead of rising gradually through the interval 3–20 Ma as before, empirical β(t) from this analysis of the PaleoSites Database seems to jump suddenly upward at 8 Ma, from the same low neotectonic level of 13%–15% that was seen previously, to the same high level of 55%–60% that was seen previously. One interpretation might be that California tectonic rates have actually been stable since 8 Ma, and that offset features of 3–8 Ma age in California can be used without any additional risk of inapplicability (although we found above that they are less reliable in the conterminous western states as a whole). Another interpretation is that some degree of preselection has been applied, because of the stated intent to represent only neotectonics in SCEC-FAD, Qfaults, and PaleoSites. That is, rates obviously disagreeing with Quaternary and/or Holocene rates may have already been edited out to some extent. If so, the apparent stability of rates beyond 3 Ma could be an artifact, and not a reliable guide to the treatment of new results. In the merging of single-feature rates to produce the combined PaleoSites (California) long-term offset rates shown in Table 2, I have made the conservative choice that offset features older than 3 Ma present increasing risk of inapplicability, as in the previous section.

DISCUSSION

In this project I have attempted to define and demonstrate automated numerical procedures for quantifying long-term offset rates that are more objective and reproducible than the committee-of-experts method previously used in seismic hazard evaluation. Because these methods yield full (estimated) probability density functions for the offset rate at any datum point (and also tentative combined pdfs for the assumed-uniform rate of any fault train), it will ultimately be possible to test their validity with new data, although it may take 25–50 yr before enough new data have been collected.

The need for more primary data is highlighted by a few summary statistics showing how much remains unknown. For example, the relative precision with which a long-term fault offset rate is known could be measured by a metric defined as the ratio of the median of the combined long-term offset rate distribution to the full width of its 95% confidence interval. For a Gaussian distribution, this metric would be the mean divided by four standard deviations. In the broader Gorda-California-Nevada orogen region (Table 1), only 21 fault trains have precisions exceeding 2 by this metric. Only 48 fault trains have precisions >1, which I take as a minimum standard for a well-constrained rate. At the other end of the precision scale, 362 active (or potentially active) fault trains in my compilation for the western conterminous U.S. have no dated offsets, so their rates are estimated entirely from the bootstrap prior distributions, and they have relative precision metrics of 0.04–0.16 (depending on fault type).

To create a comparable statistic for California, consider that there are 198 active fault sections in the Working Group on California Earthquake Probabilities (WGCEP) Fault Model 2.1, which is a mutually compatible selection of fault sections from the Fault Section Database that was used by the WGCEP in its 2007 update of seismic hazard in California. Table 2 shows 97 combined long-term rates for some of these fault sections, based on the PaleoSites Database, but only 30 of these combined rates are well constrained. These figures do not count alternate traces beyond the first. Thus, ∼30/198 ≅ 15% of active fault sections in California have well-constrained rates. This is higher than the overall western conterminous U.S. fraction of 48/849 trains ≅ 6% well constrained, but still low in absolute terms. If California is excluded, the fraction of well-constrained fault trains in the other conterminous western states, based on Table 1, drops to only 15/655 ≅ 2%.

A potentially offsetting consideration is that the faults that are well-constrained tend to be those with higher offset rates. One way to quantify this is in terms of each fault's contribution to the regional seismic potency rate (integral of long-term slip rate over seismogenic fault area). For simplicity, I assume the following downdip, seismically coupled fault widths (Bird and Kagan, 2004) throughout the conterminous western U.S.: 8.6 km for strike-slip faults, 3.7 km for normal faults, and 38 km for thrust faults. Multiplying my best-estimate (median) slip rates from Table 1 by these widths and by fault lengths from digitized traces, I crudely estimate the seismic potency rate for each fault. The result is that in California, faults with well-constrained rates contribute (12.5 m3/s)/(32.7 m3/s) = 38% of the seismic potency rate. In the western conterminous U.S. as a whole, faults with well-constrained rates contribute about (12.6 m3/s)/(36.6 m3/s) = 34% of the seismic potency rate. In the western conterminous U.S. excluding California, faults with well-constrained rates contribute about (0.055 m3/s)/(3.91 m3/s) = 1.4% of the seismic potency rate. These proportions must be interpreted with caution because their denominators depend on the slip rates of many faults whose slip rates are not well constrained. Also note that none of these figures include the enormous seismic potency rate of the Cascadia subduction zone megathrust (∼130 m3/s), which is based on plate tectonic models rather than conventional offset geologic features (Bird, 2003; Bird and Kagan, 2004).

Another interesting way to look at these results is to plot the value of the metric for all fault trains versus the number of offset features used in the combined rate, as in Figure 9. Here we see a clear increase in precision of the combined rate with increasing number of data. Only when the number of offset features reaches four is there an even chance of achieving a well-constrained combined rate, and only when the number of offset features reaches seven or more is there virtual certainty. Thus, funding for neotectonic slip-rate studies should never be withheld on the basis of concerns over redundancy when the number of published offset features is less than this.

In addition to collecting new data, the community of neotectonic geologists and geophysical modelers could improve these results more rapidly in two ways. Erroneous or inappropriate assumptions that I may have made can be corrected after appropriate debate, while keeping the overall formalism. Also, increased engagement of primary field investigators with these data compilations should result in more complete and accurate tables of primary observations and their uncertainties. A process of mutual education between field geologists and mathematical geophysicists has begun, and needs to continue.

In the immediate future, these rates and uncertainties will be used as input to kinematic finite-element modeling of neotectonics with our program NeoKinema, in order to estimate the many unknown long-term offset rates by incorporation of geodetic data, stress-direction data, and kinematic consistency. Other investigators may also wish to employ them, either as input data for alternative kinematic models, or as test data to evaluate the relative success of various dynamic models.

APPENDIX: PROOF OF EQUATIONS 6 AND 7

Let X and Y be two independent random variables, with associated probability density functions (pdfs) of pX(x') and pY(y'). Let H(X,Y) be an algebraic function that is both differentiable and monotonic (hence, invertible) with respect to each of X and Y. The pdf of H, symbolized by pH(h'), is defined as:  
formula
Also, the cumulative probability of any given value of H can be expressed as the limit of sums of probabilities in mutually exclusive cases defined by an infinite series of non-overlapping, adjacent value ranges of independent variable X. In the case ∂H/∂Y > 0 where small values of Y are associated with small values of H,  
formula
Combining equations A1 and A2,  
formula
Combining the common limit operator and factor,  
formula
If the term 1/dh' is moved so that dh' only appears in the last term, the related limit can be taken locally:  
formula
Now, recall that H is assumed monotonic in Y; therefore, an infinitesimal step in H at fixed X implies an infinitesimal step in Y. This permits us to linearize the last term:  
formula
We can now recognize the left parts of this formula as equivalent to an integration of pX with respect to dx':  
formula
The alternate case that∂H/∂Y < 0 is the same except that a negative sign appears due to the reversal of inequalities along the Y axis, and the consequent use of complements of cumulative probabilities. Both cases can be combined into one formula by writing the absolute value of the partial derivative. Thus, the pdf of function H can be expressed as a weighted convolution: where all three terms inside the integral are uniformly non-negative. This is equation 6 in the paper. Equation 7 is obtained by simply switching the generic symbols X and Y.  
formula

If you are viewing the PDF of this paper, or if you are reading this offline, please visit http://dx.doi.org/10.1130/GES00127.S1 (Slippery.f90), http://dx.doi.org/10.1130/GES00127.S2 (Slippery.exe), http://dx.doi.org/10.1130/GES00127.S3 (Table_1.xls), http://dx.doi.org/10.1130/GES00127.S4 (Table_1.txt), http://dx.doi.org/10.1130/GES00127.S5 (f_Gorda-Cal-Nev.txt), http://dx.doi.org/10.1130/GES00127.S6 (references_cited_Table_1.doc), http://dx.doi.org/10.1130/GES00127.S7 (Table_2.xls), http://dx.doi.org/10.1130/GES00127.S8 (Table_2.txt), http://dx.doi.org/10.1130/GES00127.S9 (f_GCN_Bird2007_Table1.dig), or see the full-text article on www.gsajournals.org to access the table, source code, and executable files.

Susan Perry provided advance access to the combined-events portion of the PaleoSites Database. Daniel Petrizzo reviewed all California entries from the author's database prior to their submission to the PaleoSites Database, and added geographic coordinates. Chris Wills provided preprints, helpful suggestions on data sources, and clarification of the text. Research was supported in part by the U.S. Geological Survey (USGS), Department of the Interior (under USGS award number 01HQGR0021 and earlier awards), in part by the National Science Foundation (NSF grant EAR-0336950 and earlier grants), and in part by the Southern California Earthquake Center (SCEC). SCEC is funded by NSF Cooperative Agreement EAR-0106924 and USGS Cooperative Agreement 02HQAG0008. This is SCEC contribution 1077. The views and conclusions contained in this document are those of the author and should not be interpreted as necessarily representing the official policies, either expressed or implied, of the U.S. Government, the SCEC, or the Working Group on California Earthquake Probabilities.

Supplementary data