The relationship between the instantaneous velocity field and the rate of moment release in the lithosphere

Summary

1 Introduction

Continental deformation generally occurs in broad and complex zones. It has often been noted (e.g. McKenzie & Jackson 1986) that the distribution of deformation in these zones is not uniform. Rather, fault populations occur on many spatial scales, and large areas of distributed deformation may be juxtaposed against similarly large areas best characterized as coherent, unfaulted blocks. One way of quantifying the distribution of deformation in continental deformation zones is through the rate of moment release, a measure of the density of the rate of production of elastic dislocations. Until the mid-1980s, estimation of moment release rate was based primarily on geological and palaeoseismological observations and, where possible, the historical record of seismicity (e.g. Wesnousky 1982). Since the advent of space-based geodesy, particularly the Global Positioning System, it has been possible to determine instantaneous surface velocity fields with increasing density and precision, to the point where lateral variations in strain rate can be mapped to a certain degree of accuracy. These strain rates are often related directly to the moment release rate and compared with other estimates of the moment release rates. However, I am not aware of a theoretical basis for this relationship. This is not surprising because the meaning of the instantaneous velocity field is intricately bound up in the mechanics of continental deformation, and opinions still diverge on appropriate conceptual models. A relationship exists between the steady-state strain rate and the moment release rate (Kostrov 1974), but a corresponding formalism connecting instantaneous strain rates with moment release rates is lacking. The purpose of this paper is to propose a new relationship between the instantaneous velocity field and the moment release rate in a manner which is, as far as possible, in harmony with prevailing ideas on the mechanics of continental deformation.

2 End-Member Models Of Continental Deformation

Active continental deformation is usually envisioned in terms of two end-member models (e.g. King 1994; Nyst 2001). In the first, the block model, it is supposed that the continental lithosphere consists of coherent, unfaulted elastic blocks moving with respect to one another (e.g. Savage & Burford 1973; Savage 1979; Hashimoto & Jackson 1993). Over very long time periods (averaged over many seismic cycles), relative motions among the blocks are accommodated by motions across the faults which divide them, and the blocks themselves behave rigidly. The block model accounts for the elastic effects associated with strain accumulation localized on the block boundaries, this strain being determined by a model of resistance to the relative motions above a locking depth (Savage 1979; Matsuúra 1986; Hashimoto & Jackson 1993). The second model is the ‘thin-sheet’ model (England & McKenzie 1982; England & Molnar 1997; Flesch 2000). It is a continuum model in which the lithosphere is defined to include elastic portions of the crust and mantle plus intervening weaker material (i.e. lower crust). This lithosphere is assumed to overlie an asthenosphere that is much weaker, i.e. it flows much more readily than any components of the lithosphere. Boundary forces acting on the thin lithosphere, such as the static force produced by lateral variations in gravitational potential energy, can be related to consequent steady-state deformation by means of an effective viscosity (e.g. Flesch 2000), which relates the induced lithospheric stresses to the lithospheric strain rate.

In the block model there is no memory from past earthquakes: instantaneous strain accumulation effects depend only on the accumulating ‘slip deficit’ on locked portions of the lithosphere at rates that bear no clear relation to the occurrence of past earthquakes. Nor in the thin-sheet model are earthquakes explicitly taken into account. The strong (elastic) portions of the lithosphere are assumed to have zero strength in the long term because the yield strength of these regions is continually being exceeded, resulting in distributed cataclastic failure. While this assumption is reasonable and implies the occurrence of earthquakes, admissible spatio-temporal effects from past earthquakes are quite limited by this simple characterization of the average strength of the elastic portions of the lithosphere.

Both the block model and the thin-sheet model are useful descriptions of continental deformation under the appropriate conditions, i.e. for the block model, the existence of large unfaulted blocks and for the thin-sheet model, steady-state continuum deformation of the lithosphere as a whole. However, they afford no insight into the properties of earthquakes that occur over time. To estimate the distribution of seismic moment release, the block model requires that the block boundaries be precisely defined and that the instantaneous velocity field be the product of specific imposed boundary motions (e.g. slip at depth) without a memory from past earthquakes. This approach may be valid in the case of an infinitely long, straight fault (see the comments below) but will lead to biased estimates in more complicated fault geometries. The effective viscosity defined in the thin-sheet model relates driving stresses to steady flow over the entire lithosphere, but it is not clear (and usually subject to great approximation) how stress in the ‘averaged’ lithosphere should be partitioned into strain rate within the elastic versus anelastic portions of the lithosphere.

There is a distinction between instantaneous velocity and steady-state (i.e. secular) velocity that cannot be overemphasized. The former is typically represented by a geodetic survey measured over a short time period, whereas the latter is best represented by the geological slip rates obtained for the faults occupying a region. Nevertheless, in a modelling context many instances could be cited in which the two types of velocity field are confused. The difference can be demonstrated by pointing out the connection between the instantaneous velocity and interseismic velocity in the context of an earthquake cycle. In a single-fault system let us consider the interseismic velocity field averaged over the entire interseismic period between successive earthquakes (Fig. 1). This generally differs from the steady-state velocity, which represents an average over a time interval that includes both the interseismic period(s) and the static offset(s) from the earthquakes themselves (Fig. 1). The distinction between ‘interseismic velocity’ and ‘steady-state velocity’ is clear in treatments of an elastic—viscoelastic coupled system. For example, in the viscoelastic coupling model of Savage & Prescott (1978), slip events occur periodically on an infinitely long strike-slip fault embedded in an elastic layer overlying a viscoelastic substrate. In this system, Savage (1999) note that the average interseismic velocity field exactly replicates the solution for strain accumulation around a strike-slip fault that is locked above a certain depth (Fig. 2a; Savage & Burford 1973), whereas the steady-state velocity is simply steady block motion of one side of the fault with respect to the other side (Fig. 2b). Clearly, the interseismic velocity field is the more appropriate description of the instantaneous velocity field.

Figure 1

A repeating earthquake occurs at a given location with periodicity T. The expected interseismic velocity field is the average velocity during an interseismic period (i.e. time 0⁺ to T⁻). The steady-state velocity is the average over a long time interval which includes several earthquakes.

Open in new tabDownload slide

Figure 2

Fault-parallel velocity in the viscoelastic coupling model (Savage & Prescott 1978), assuming uniform periodic slip on an infinitely long, straight strike-slip fault which penetrates the entire elastic layer of thickness H. (a) Expected interseismic velocity after the system has evolved through several cycles. (b) Steady-state velocity.

Open in new tabDownload slide

I adopt the principle that instantaneous surface velocity measurements represent the interseismic velocity field, not the steady-state velocity field. In addition to the implications for mechanics of continental deformation discussed below, this means that the procedure of relating geodetically determined velocity fields directly to rates of seismic moment release, as is done in numerous studies, is not altogether correct: the former is the product of the interseismic velocity field, whereas the latter is the product of the steady-state velocity field. This may help rationalize systematic discrepancies in estimates of seismic moment release using different methods (King 1994; Ward 1998).

3 Periodic Moment Release Model of a Elastic—Viscoelastic Coupled System

3.1 Model description

To capture the properties of earthquakes occurring within the strong portion of the lithosphere, I propose a physical model intermediate between the block model and thin-sheet model. The Earth structure is assumed to consist of a strong elastic lithosphere (elastic crust and possibly uppermost mantle) underlain by a weaker ductile substrate (lower crust and/or upper mantle). This picture is supported by studies of rheology structure in the western United States (Pollitz 2001, and references therein) and Mongolia (Vergnolle 2001). Under the framework of elastic—viscoelastic coupling, the time-dependent velocity in the lithosphere can be related to one or more particular earthquake sources. I postulate that many patterns of active continental strain are the product of memory from past earthquakes, i.e. the compounded effects of stress diffusion (post-seismic relaxation) from numerous past earthquakes distributed throughout a region. The likely existence of such earthquake memory is supported by the theoretical study of Bott & Dean (1973), who demonstrated that plate-like behaviour remote from a plate boundary is the product of stress diffusion into the plate interior from continually occurring moment release along a long plate boundary. This result has been verified in special geometries, i.e. 2-D strike-slip faulting (Savage & Prescott 1978). Although the precise record of earthquake occurrence in a locality, let alone for a region, is generally unknown, it is possible to relate the interseismic velocity field to the rate of local seismic moment release by assuming that this moment release occurs periodically. The periodicity may vary from one location to the next. As demonstrated in the following section, this provides a framework for relating the instantaneous velocity field to the time-averaged distribution of moment release rate in the lithosphere.

3.2 Interseismic velocity field

Let a fault occupying a portion of lithosphere undergo a slip event with average recurrence time T. For convenience I assume that the slip events occur periodically, i.e. T is the time between all successive slip events. As discussed below, however, the principal result remains valid even for non-uniform recurrence intervals as long as the expected recurrence interval is T. The expected interseismic velocity (EIV) at a particular point is defined as the velocity averaged over interseismic times, i.e. between time 0⁺ just after the last event to time T⁻ just before the succeeding event (Fig. 1). If u(t) and v(t) are the time-dependent displacement and velocity, then

(1)

Let u_n(t) be the displacement field on the viscoelastic system for the nth previous slip event. Then

(2)

For times t > 0, u₀ contains the effects of the coseismic (purely elastic) displacement field from the event which occurred at time 0 plus accumulated post-seismic effects up to time t.

For times t > 0

(3)

Evaluating v_s with eqs (1) and (3) yields

(4)

Thus the EIV is the difference between the fully relaxed response and the purely elastic response divided by the interseismic time interval.

This can be generalized to the case of the response to moment release distributed throughout the lithosphere and a 3-D interseismic velocity field. Suppose a given volume δV of lithosphere has moment release M occurring periodically with repeat time T. If G_ij(r, r₀; t) is the displacement field observed at r at time t on the viscoelastic system from unit moment tensor excitation M_ij= 1 acting at r₀ beginning at time 0, then the EIV field at point r is

(5)

It is convenient to define the moment rate tensor as and the moment rate tensor density (per unit volume) as ⁠. To account for the additional effects of a distribution of time-dependent forces, as could arise from tectonic interactions at a collision zone or shear zone, I further define G_i(r, r₀; t) as the displacement field observed at r at time t on the viscoelastic system from unit force excitation F_i= 1 acting at r₀ beginning at time 0; I then define a forcing rate density (assumed constant with time at any given location) ⁠. Combining all of the above definitions with eq. (5), the total EIV field is

(6)

This provides a direct relationship between the EIV field and the distribution of time-averaged moment release rate and forcing rate in the lithosphere, in principle, for all six components of the moment rate tensor and all three components of the forcing rate vector. In the remainder of this paper I will focus only on the contributions to eq. (6) from ⁠. As will be demonstrated in subsequent sections, eq. (6) forms the basis for inversion of elements of the moment rate density tensor from observed surface velocity fields. In practice, given the nature of GPS velocity fields, one should expect to recover only the distribution of the vertical average of ⁠.

A noteworthy property of the derived relationship is that neither the block nor continuum extreme models of deformation need be satisfied, i.e. both discontinuous or continuous distributions of in eq. (6) are admissible sources of the deformation. Another property of the relationship is that, for a given subdivision of the Earth into elastic and anelastic regions, v_s does not depend on the viscosity distribution but only on the elastic moduli distribution. The average recurrence interval T need not be constant: it may vary from one location to the next. This allowable variation of T with r₀ is implicit in the integral (6) through its dependence on ⁠. The relationship (6) is valid even for a fault system with a more complicated distribution of interevent times. More precisely, for a single fault let the expected number of events occurring within a preceding time interval nT be n. The history of events with this property need not be periodic, only the average recurrence interval must equal T. Then one can derive the expected interseismic velocity field at time t= 0 resulting from a sequence of events occurring within −nT < t < 0. The expected contribution from one event occurring within that time interval is

The contribution from n such events reduces to eq. (4) in the limit n−→∞, leading again to eq. (6).

It will be convenient to work with vertically integrated rates of moment release

(7)

where R is Earth's radius and H is a depth assumed to be 15 km. If are assumed to be independent of r and assuming no force contributions, then eq. (6) may be rewritten as

(8)

The assumption of constant with depth in the ‘elastic’ portion of the lithosphere may be a good approximation in many situations according to modelling of the steady-state moment release rate in a depth-dependent viscoelastic medium with account for frictional rheology of faults (Rolandone & Jaupart 2002).

3.3 Steady-state velocity field

The steady-state (secular) velocity is given by the EIV plus the time-averaged contribution of static offsets from the periodic moment release events (Fig. 1) and constant forces. The time-averaged static component is given by

(9)

The sum of the velocity fields in eqs (6) and (9) is then the steady-state velocity:

(10)

If are assumed to be independent of r, and assuming no force contributions, then eq. (10) may be rewritten as

(11)

There are major differences between the EIV field (eq. 6) and the steady-state velocity field (eq. 8), not just in their mathematical form but in their interpretation. The EIV field represents the time-averaged difference between the fully relaxed deformation field and the static deformation field contributed from a distribution of elastic dislocations, whereas the steady-state velocity field represents just the fully relaxed response—the time-averaged static deformation field on the elastic Earth stripped of its viscoelastic portions.

Although the formulae for interseismic and steady-state velocity presented in this and the previous section are valid for general viscoelastic models, in all applications to be presented I will evaluate these velocity fields on a spherically stratified viscoelastic medium. The required Green functions evaluated at times 0⁺ or ∞ after a particular moment release event are calculated with the methods described in Pollitz (1996, 1997).

3.4 Context of the EIV model

A secular force model, such as that produced by the relative motions of shear zone boundaries in a strike-slip fault system or by motion of a downgoing slab at a subduction zone, corresponds to non-zero ⁠. Even in the absence of slip events (⁠⁠), eq. (6) shows that the secular force model will generate a non-trivial average instantaneous velocity field, and eq. (10) demonstrates the same for the steady-state velocity field. In a static force model such as a model of lateral variations in gravitational potential energy, ⁠, and in the absence of earthquakes there would be no instantaneous or steady-state deformation of the lithosphere: the system would reach static equilibrium. However, the topographic gradients generating lateral gradients in gravitational potential energy typically produce deviatoric stress levels much larger than the yield strength of the lithosphere; earthquakes must occur. The same conclusion may obviously be reached for the stress levels generated by a secular force model, e.g. the continual occurrence of earthquakes within the San Andreas or Anatolian strike-slip fault zones. Thus, any attempt to explain strain gradients over the short-term must retain the short-term integrity of the lithosphere (i.e. its elasticity) and thus must account for the continual release of strain through discrete slip events.

Models that describe short- or long-term deformation could each be classified as either deterministic or statistical. For long-term deformation, the first category is represented by steady block motions of tectonic plates or microplates described by rigid rotations. The second category is best represented by the thin-sheet model, in which no attempt is made to describe every fault which ruptures in a large region populated by a very large number of distributed faults. Rather, an effective viscosity (Flesch 2000) describes the expected behaviour of the lithosphere over long periods of time. Over shorter time intervals, a deterministic model such as the block model or the viscoelastic coupling model (Savage & Prescott 1978) is applicable only within a laterally restricted region where the number of important faults is small. Over short time periods but covering a large region, however, the effects of slip release over a large number of faults must be accounted for, and a statistical theory for interseismic deformation is required. The EIV model of eq. (6) provides a statistical interpretation of the interseismic deformation without demanding detailed knowledge of all of the faults occupying the given region.

3.5 Kostrov's formula

Kostrov (1974) derived a formula that is often used to relate rates of moment release to the strain rate field. As another way of emphasizing the difference between the interseismic and steady-state velocity field, I investigate under what conditions this formula is valid. Kostrov's formula states that the average strain rate in a volume is proportional to the average moment release rate in the same volume, both averaged over time and the space in the volume. Kostrov assumed shear dislocations populating a volume of constant shear modulus µ, in which case the tensor form of the formula is

(12)

where is the considered volume.

One may obtain a more general form of this formula as follows. Let an elastic volume Γ bounded by area Ω enclose a large number of seismic sources and consider the deformation rate within that volume. For a single seismic source at r₀ with moment rate tensor ⁠, the force balance equation that describes the quasi-static velocity field is

(13)

where is the stress rate tensor and δ(r−r₀) is the 3-D Dirac delta function. For an arbitrary vector function g(r) we have

(14)

Simple differential relations applied to the integrands yield

(15)

In eq. (15) I eliminated the second integral on the right-hand side using Gauss' theorem and noting the vanishing of the δ-function on Ω. Let us suppose for the source element at r₀ that ⁠, and consider the contributions from distributed over the volume Γ. Applying Gauss' theorem to the first integral on the left-hand side of eq. (15), we obtain

(16)

Eq. (16) is quite general and is valid for any volume enclosing an elastic portion of an elastic—viscoelastic coupled system. I now consider the deformation averaged over a time much longer than the seismic cycle of either Γ or the surrounding actively deforming regions. Assuming that the boundary Ω does not intersect any dislocation source, in the first term on the left-hand side of eq. (16) must be zero because, over a sufficiently long time period, all built-up stresses are relieved by the distributed faulting. Since eq. (16) is valid for arbitrary choices of g, we have

(17)

If a constitutive relation for a homogeneous isotropic elastic medium is assumed, i.e.

(18)

then this is a restatement of eq. (12), i.e. the proportionality between the volume-averaged strain rate and the moment release rate in a homogeneous isotropic elastic volume. However, eq. (17) is valid for a more general constitutive relation.

Note that eq. (17) is strictly applicable only to the steady-state strain or velocity field. During an interseismic interval or any time period in which elastic stresses are not averaged out, ⁠, in the Ω-integral of eq. (16) may be large, in which case eq. (17) is not justified.

4 Synthetic Model of Distributed Deformation

4.1 Elements of synthetic model

Here I derive an EIV field based on a random distribution of faulting within a 20 × 20 deg² area (Fig. 3). The H-shaped grey region contains a network of both tensile-opening faults (region E) and strike-slip faults (regions A—D). The relationship between the fault geometry and moment tensor for these fault types follows the conventions of Aki & Richards (1980). Assuming vertical faults, these sources are defined as follows. If the normal and horizontal tangent and vertical tangent vectors to an infinitesimal section of the fault plane of area δA are given by ⁠, and ⁠, respectively, and the slip value is u, then the non-trivial moment tensor components of that section are

(19)

where δ and μ are the Lamé parameters.

Configuration of synthetic fault network. It is composed of 300 north—south-trending tensile opening faults randomly distributed in region E and four sets of 100 strike-slip faults randomly distributed in regions A—D. Each fault ruptures periodically with a random recurrence interval in the range 500–2500 yr and random time of the most recent earthquake. Post-seismic relaxation from the compounded relaxation of these faults is observed at 100 sites (top part of figure) randomly distributed within the boxed 14 × 14 deg² area. Calculations are performed in an r-θ-φ spherical geometry.

Open in new tabDownload slide

The faulting is distributed in both space and time and is intended to represent a stationary process of earthquake occurrence within any given active fault zone (i.e. one of the regions A—E) and capture the breadth of distributed deformation that makes up a broadly deforming region such as the central Basin and Range province. On any given fault, however, earthquake occurrence is constrained to be periodic in time. In detail, I distribute 700 synthetic source faults Γ_ν according to the following scheme.

• 300 tensile opening source faults are randomly distributed in region E under a uniform distribution. Each source fault is constrained to have a strike of 0°, a length of 30 km and the time history of earthquake occurrence on a particular fault follows the probability distributions f_T and f_t given below. Tensile opening at the time of each earthquake on any given fault is 1.15 m.

• 100 strike-slip faults populate each of regions A—D under a uniform distribution. Each source fault is constrained to have a strike of 90°, a length of 30 km, a width of 15 km, and follow the f_T and f_t time distributions of occurrence. Regions A and D are populated with right-lateral faults; regions B and C with left-lateral faults. Strike-slip motion at the time of occurrence of each earthquake on any given fault is 1.725 m.

• Each fault ruptures periodically, and I associate with each fault Γ_ν a recurrence interval T_ν and elapsed time since the previous event t_ν. The T_ν are considered identically distributed random variables with probability density functions f_T(T) = 1 for 500 < T < 2500 yr and f_T(T) = 0 elsewhere. The elapsed time t_ν since the last event is uniformly distributed between the present, which is arbitrarily assigned time t= 0, and time t=T_ν, i.e. it has a probability distribution f_t(t) = 1 for 0 < t < T_ν and f_t(t) = 0 elsewhere.

The recurrence intervals of 500–2500 yr are roughly consistent with those determined on several faults in the Basin and Range province and environs (e.g. Anderson & Hawkins 1984; Wills & Borchardt 1993; McCalpin & Nishenko 1996; Sanders & Slemmons 1996; Reheis & Sawyer 1997; Bell 1999), although the range in recurrence times is generally even broader because a few faults are much less active than others. All faults are constrained to rupture over the depth interval 0–15 km with uniform slip.

The described scheme is realistic in allowing a distribution of faulting on a broad range of spatial scales, in producing earthquake occurrence which, when viewed over sufficiently long time intervals, is essentially Poissonian, and in accounting for post-seismic relaxation from multiple past events. It is unrealistic in chiefly two respects:

• the earthquake occurrence on an individual fault is assumed to be periodic;

• the constant length and slip of the faults leads to a moment release distribution that does not obey a magnitude—frequency relationship.

The assumption in the first case is largely inconsequential because (as will be demonstrated in Section 4.2) different realizations of the synthetic fault catalogue, each of which results in quite different slip histories on individual faults, yield nearly the identical instantaneous velocity field. The second case arises because the faults populating a region will generally exhibit a wide variety of slip and fault areas, ranging from microearthquakes to great earthquakes, and even a single fault segment may not necessarily slip with identical magnitude in each event. The chosen system is considered a proxy for a more complicated system that would include a realistic size distribution of earthquakes. For example, the system of 300 tensile-opening faults occupying region E of length 30 km, slip 1.15 m and average recurrence interval 1500 yr, implies the occurrence of exclusively magnitude 6.8 earthquakes approximately once every 5 yr. This could be replaced with a system of 150 faults of length 30 km and slip 1.15 m and 15 faults of length 100 km and slip 3.45 m. With the same mean recurrence interval of any given fault (1500 yr), this would produce M= 6.8 events once every 10 yr and M= 7.6 events once every 100 yr. The behaviour of this more complicated system deserves investigation, but I will henceforth restrict attention to the simpler system.

To summarize, one realization of earthquake occurrence in this model consists of the set {Γ_ν, T_ν, t_ν|ν= 1700}, which defines an infinite sequence of past earthquakes. Time-dependent deformation results from the compounded post-seismic relaxation from every earthquake in the sequence on a viscoelastic earth model, shown in Fig. 4. This model possesses a lower crust of viscosity η_c= 5 × 10¹⁹ Pa s and a mantle of viscosity η_m= 2 × 10¹⁹ Pa s. This simple viscosity structure is similar to that obtained by Bills (1994) for the central Utah region and by Pollitz (2001) for the longer-term (>2 yr post-seismic) relaxation in the Mojave Desert, CA.

Figure 4

Assumed distribution of bulk and shear moduli with depth. Viscoelastic structure beneath an elastic upper crust of thickness H= 15 km is parametrized with uniform lower crust and mantle viscosities η_c and η_m, respectively.

Open in new tabDownload slide

Time-dependent post-seismic velocity at time t after a slip event consists of a linear combination of thousands of viscoelastic normal modes (Pollitz 1997), each with exponentially decaying time dependence exp(−s_jt), where 1/s_j is the characteristic decay time of mode j. Fig. 5 shows the toroidal and spheroidal mode dispersion on the viscoelastic model of Fig. 4. There are two toroidal and seven spheroidal mode branches. One property of the dispersion curves is that decay times are larger at longer wavelength. At wavelengths ≳100 km, the spheroidal decay times reach the same order of magnitude as the recurrence intervals assigned to the faults (grey regions of Fig. 5). The toroidal modes, which are less dominant at long wavelength, have decay times of ∼100 yr at 100 km wavelength. This means that, at a given observation point sufficiently far from the source, interseismic velocity contributed by a single dislocation exhibits only moderate variations with time during an interseismic period (e.g. Savage & Prescott 1978). It might be expected that an ensemble average (e.g. Kittle 1958) of post-seismic relaxation from many ruptures distributed in space and time, as is described by the EIV model, will exhibit a fairly uniform variation with time.

Figure 5

Toroidal and spheroidal mode dispersion curves for viscoelastic normal modes derived from the structure of Fig. 4. The dispersion includes the effect of coupling with Earth's gravitational acceleration (Pollitz 1997). The range of recurrence intervals on the considered faults is shown in the grey area.

Open in new tabDownload slide

Let m_ν(r₀) be the moment tensor density of the νth fault in the synthetic catalogue. The instantaneous velocity field at selected observation sites r_k for fault Γ_ν for an event at time t before the present is evaluated in the form

(20)

For a given fault geometry and viscoelastic stratification, this deformation is calculated using the method of Pollitz (1997). The total instantaneous deformation contributed by all faults rupturing with their respective periodicities is then

(21)

The slip values listed in steps (1) and (2) above have been chosen such that the secular velocity field of the region will behave as rigid opening of the two sides of the tensile opening zone (region E) apart from one another. The rate of opening equals the product of one fault length (30 km), the number of faults (300) and the slip value (1 m) divided by the length of the rupture zone (1112 km), multiplied by the average inverse recurrence interval, which is

(22)

This yields a tensile opening rate of 7.5 mm yr⁻¹. A similar calculation for the secular motion across each strike-slip segment yields 3.75 mm yr⁻¹ in the right-lateral or left-lateral sense, i.e. one-half of the tensile opening rate, which is appropriate for the strike-slip faults that bound the northern and southern edges of the tensile opening zone. The secular motion thus consists of equal but opposite block motion of the area bounded by regions A, E and C with respect to the area bounded by regions B, E and D. The secular block motion is not perfect because of the statistical nature of the fault distribution and the termination of the strike-slip fault zones at 1000 km length (perfect block motion would require infinite length for the strike-slip fault zones).

Given the distribution of faulting and the recurrence intervals, one may construct maps of the distribution of moment release rate. Such maps were constructed on the 14 × 14 deg² area outlined in Fig. 3 by expanding the moment release rate in terms of Hermite—Gauss (HG) functions (see eq. 24) up to degree 20. Under the spherical geometry of Fig. 3, Fig. 6 shows the distribution of the two independent synthetic moment release rates per unit area:

(23)

where H= 15 km. These smooth distributions show that the employed random distributions of faulting well represent the various fault zones (regions A—E) with nearly uniform sampling.

Figure 6

Representation of the two independent moment release rates in terms of smooth Hermite—Gauss functions up to degree 20. These distributions were derived from the random distribution and geometry of synthetic faults (Fig. 3) and their corresponding recurrence intervals. Contour lines are drawn (in units of 10¹⁴ N m km⁻² yr) at values −0.1, −0.05, +0.1, +0.25 for and at ±0.05, ±0.1 for ⁠.

Open in new tabDownload slide

4.2 Instantaneous velocity field

Fig. 7(a) show the instantaneous velocity field v^obs_k at 100 randomly distributed surface observations sites r_k within the 14 × 14 deg² subregion (Fig. 3) calculated from eq. (21) for one set of random fault realizations. One remarkable property of the velocity distribution is that, sufficiently far from any of the fault zones, the velocity field exhibits roughly block motion similar to the secular motion. This may not be entirely surprising since it is known in similar geometries that the compounded effect of relaxation from past earthquakes yields such block motion (Bott & Dean 1973; Savage & Prescott 1978; Pollitz 2001). It demonstrates, however, that large instantaneous velocities are generated far from the sources of faulting even though faulting events occur infrequently. Figs 7(b) and (c) show the velocity field derived from the same distribution of faulting but on a viscoelastic structure with η_m= 8 × 10¹⁸ Pa s (Fig. 7b) or η_m = 6.4 ×10¹⁹ Pa s (Fig. 7c), assuming η_c/η_m= 2.5. The amplitude and pattern of synthetic velocity is generally the same as for the case of intermediate viscosity in Fig. 7(a). Fig. 8 shows the velocity field calculated from another realization of faulting with the structure η_m = 2 × 10¹⁹ Pa s, η_c/η_m= 2.5. It exhibits the same general pattern as that in Fig. 7(a), showing that the sampling of our random fault model space is dense enough to yield a robust velocity pattern. This similarity further suggests that the regional response to distributed faulting is statistically insensitive to the exact range of recurrence intervals assigned to the faults. For example, one fault with an expected recurrence interval of 1500 years could be replaced with two faults (in close proximity to one another) of expected recurrence interval 3000 years without appreciably altering the velocity field; both realizations are associated with the same rate of moment release.

Figure 7

Black arrows: the synthetic velocity field v^obs_k at 100 randomly distributed observation sites r_k evaluated for an infinite sequence of past earthquakes (eq. 21) using the random fault network parametrized by {Γ_ν, T_ν, t_ν|ν= 1700}. Different cases correspond to relaxation on the Earth model with viscosity: (a) η_m= 2 × 10¹⁹ Pa s, (b) η_m= 8 × 10¹⁸ Pa s, (c) η_m= 6.4 × 10¹⁹ Pa s. A constant crust-to-mantle viscosity ratio η_c/η_m= 2.5 is assumed. Grey arrows: the velocity fields derived from inversion of the synthetic velocity fields, calculated from the corresponding distributions of and in Figs 9(a)-(c) and eq. (8).

Open in new tabDownload slide

Figure 8

The velocity field v^obs_k sampled at 100 observation sites as explained in Fig. 7 with η_m= 2 × 10¹⁹ Pa s and η_c/η_m= 2.5, but with a different realization of the random fault network.

Open in new tabDownload slide

The representation of the EIV field derived in Section 3.2 is an ensemble average (e.g. Kittle 1958) over all possible realizations of faulting. If the actual system is described by rupture of numerous faults at random past times and a viscoelastic model, as in this example, then one could describe the phase space of the system in terms of a range of velocity distributions that could result from such faulting. One could, in principle, define the entropy of such a system in terms of the number of realizations of faulting that would yield a restricted range in phase space; this would depend upon the number of faults and their recurrence intervals as well as the viscoelastic structure. I will not pursue this avenue of analysis further except to note that the results to be presented in Section 5.2 suggest that entropy tends to increase with increasing viscosity, other factors being equal. In other words, the system tends towards statistical equilibrium as the relaxation time(s) of the structure become large.

The velocity field depicted in Fig. 7(a) could well represent the velocity field measured by a geodetic survey in a continental interior over a short time interval. For example, similar velocity gradients are observed in the Basin and Range province (Thatcher 1999) and Mongolia (Calais 2002), regions far from a major plate boundary. As long as this time interval did not include any earthquakes it could be interpreted as the interseismic velocity field. Given the statistical nature of the faulting distribution and times of occurrence, it is appropriate to interpret it as the expected interseismic velocity field. If the EIV is an accurate description of this velocity field, then eq. (6) should provide a good relationship between this velocity field and the various ⁠.

5 Inference of Moment Release Rates

5.1 Inversion

Using the relationship given in eq. (6) as the solution to the forward problem, one may invert a regional velocity field for the distribution of moment release rates in the elastic portion of the viscoelastic model, i.e. the upper 15 km of the model shown in Fig. 4. Vertical resolution, however, is poor, and therefore I assume that each is independent of depth. I then seek to estimate the vertically integrated moment release rates defined in eq. (18). The relationship between the EIV field v_s(r) and the moment release rates is given by eq. (8). It is convenient to parametrize this distribution over a rectangular (x, y) grid with HG functions as follows:

(24)

where l +m≤ l_max for a fixed l_max= 20, the h_m are normalized Hermite polynomials such that

(25)

and L₁ and L₂ are taken to be proportional to the dimensions of the rectangular grid. In the applications presented here L₁=L₂ are chosen such that ±7°/L₁ equals the last local maximum of the HG function of degree l_max. (There is a one-to-one correspondence between points on this rectangular grid and points on the boxed area of Fig. 3 in a spherical geometry.)

Suppose that we have a surface velocity field {v^obs_k} observed at K corresponding sites {r_k|k= 1, K}. In the inverse problem I minimize a functional of the form

(26)

The first term represents the data misfit and the second term represents the roughness of the lateral gradients in the distribution of moment release rates, weighted by α. The elements of C can account for a priori covariance structure among the synthetic observables; for simplicity I choose uncorrelated data with unit standard deviation: C_kk′=δ_kk′. With parametrized by eq. (24) and v_s by eq. (8), minimization of eq. (26) with respect to the model parameters a_ijlm leads to the normal equations

(27)

where β and γ span the set of model parameter indices {i, j, l, m}. The Frechet derivatives (terms in parentheses in eq. (27)) are determined straightforwardly from eq. (26) and integral and recurrence relations among the HG functions. In particular, the contribution of the roughness term (i.e. that proportional to α in eq. (26)) is given explicitly by eq. (19) of Friederich & Wielandt (1995).

In principle, given a data set {v^obs_k} one should solve the inversion problem for all six components of the vertically averaged moment release rates ⁠. In practice, however, one is likely to have a priori constraints based on geology and seismotectonics. I assume that, when modelling synthetic velocity fields such as those presented in Figs 7 and 8, the style of regional faulting is known to be strike-slip faulting on east—west-trending vertical faults and rifting across north—south-trending vertical faults. This allows us to restrict our model space to two independent moment release rate elements. Using the synthetic velocity fields presented in the previous section as data, I solve the inversion problem (26) for the distributions of and over the 14 × 14 deg² boxed area indicated in Fig. 3. In order to model a rifting source I impose the constraint (Aki & Richards 1980)

(28)

where 〈[]〉 denotes the value of the bracketed quantity averaged over the depth range 0-H of the model. The smoothing parameter α is uniquely chosen such that the variance reduction achieved with the inverted model with respect to the data is 97 per cent. This provides us with an objective method of comparing the results derived from different data sets.

5.2 Results

Figs 9(a)-(c) shows the patterns of inverted using the velocity fields of Figs 7(a)-(c), respectively, as the observables. The sets of grey arrows in Figs 7(a)-(c) show the corresponding inverted velocity fields. The inversion results depend, of course, upon the viscosity structure of the Earth model because of the viscosity dependence of the synthetic velocity fields. For the represented viscosity range, the inverted distributions of moment release rates recover well the input distribution (Fig. 6). I define the correlation between inverted (inverted) and synthetic input (input) as

(29)

For viscosity structure corresponding to Figs 7(a) and 9(a)(η_m = 2 × 10¹⁹ Pa s, η_c= 5 × 10¹⁹ Pa s), the correlations between inverted and synthetic moment release rate are 0.83 for and 0.74 for ⁠. Fig. 10 shows the dependence of these correlations on mantle viscosity η_m under the assumption that η_c/η_m= 2.5. For η_m≥ 0.6 × 10¹⁹ Pa s, the correlations with input and are ≥0.8 and ≥0.7, respectively. For η_m≳ 2 × 10¹⁹ Pa s, the correlations flatten as a function of viscosity and maintain relatively high values as viscosity is further increased. The reconstructive power of the method, as measured by the correlations, is poorer at smaller viscosity values. This is because in that case the velocity distribution is more sensitive to those faults which ruptured most recently and therefore contribute a larger post-seismic relaxation signal early into their next cycle. This leads to slightly more short-wavelength structure in the synthetic velocity fields. A more important effect at low viscosity is the tendency for much of the signal from past earthquakes to have relaxed prior to the observation time. This low-amplitude bias leads to somewhat lower synthetic velocities (compare Fig. 7a with Figs 7b and c) and consequently a steep fall-off in correlation as viscosity is lowered (Fig. 10).

Figure 9

Distribution of inverted moment release rates derived from inversion of the synthetic data sets of Figs 7(a)-(c) with the EIV model (eq. 8); separate distributions were derived from each velocity field. Different cases correspond to relaxation on the Earth model with viscosity: (a) η_m= 2 × 10¹⁹ Pa s, (b) η_m= 8 × 10¹⁹ Pa s, (c) η_m= 6.4 × 10¹⁹ Pa s. A constant crust-to-mantle viscosity ratio η_c/η_m= 2.5 is assumed. In each case, the roughness term of the misfit function (eq. 26) has been uniquely chosen such that the inverted velocity field yields 97 per cent variance reduction with respect to the synthetic velocity field. Contour lines are drawn (in units of 10¹⁴ N m km⁻² yr) at values −0.1, −0.05, +0.1, +0.25 for and at ±0.05, ±0.1 for ⁠.

Open in new tabDownload slide

Figure 10

Correlation between distributions inverted assuming either the EIV model (solid lines) or the secular model (dashed lines) and the corresponding HG representation of the input distributions (Fig. 6), as a function of η_m. A constant crust-to-mantle viscosity ratio η_c/η_m= 2.5 is assumed.

Open in new tabDownload slide

5.3 Secular velocity model

An alternative interpretation of the synthetic velocity fields is that they represent the secular velocity field. Assuming that this is the case, one may re-cast the inverse problem by using eq. (11) to represent the instantaneous velocity field. I have inverted the synthetic velocity fields as before, using eq. (11) in the intermediate step to relate the Δv_k to the HG expansion coefficients. The resulting distributions of moment release rate at η_m= 2 × 10¹⁹ Pa s are shown in Fig. 11, and the synthetic to inverted distribution correlations as a function of η_m are presented in Fig. 10. These figures demonstrate that the reconstructive power of the secular model is much poorer than the EIV model at all viscosities considered. There are large artefacts in both and far away from the locus of their true maxima (Fig. 6). The difference between the EIV model results (Fig. 9) and the secular model (Fig. 11) is a consequence of the very different characteristics of the Green functions that define the instantaneous velocity in the two approaches (eqs 8 and 11). In the secular model (eq. 11), the response to sources of deformation are largest at and attenuate substantially with distance from ⁠. In the EIV model (eq. 8), the corresponding response is, in fact, identically zero at and increases to a maximum at a certain distance from (a few elastic plate thickness according to calculations of Pollitz 1997) before falling off very gradually with increasing distance. In other words, in order to fit a given velocity field the sources of deformation in the secular model must be placed near the significant strain gradients, whereas the EIV model does not necessarily require sources in the straining regions themselves because of the long-range property of its Green functions. Thus, some of the artefacts in the moment release distributions in Fig. 11 are located at the western/eastern margins (longitude near ±7°), where the synthetic data field (Fig. 7a) has a large amplitude despite being far away from the true sources of deformation. The inaccuracies produced by placing sources in these remote areas is compensated by severe artefacts placed in other regions that happen to be devoid of data (for example, at latitude 2°N, longitude 3–5°E).

Figure 11

Distribution of inverted moment release rates derived from inversion of the synthetic data set of Fig. 7(a) (η_m= 2 × 10¹⁹ Pa s and η_c/η_m= 2.5) with the secular model (eq. 11).

Open in new tabDownload slide

6 Conclusions

The instantaneous velocity field is interpreted here as an ensemble average: the summed contributions of post-seismic relaxation from numerous dislocation sources distributed in the lithosphere, compounded over a history of quasi-periodic previous ruptures. The viscoelastic structure controls the temporal scale of the velocity field. The characteristic decay times of the viscoelastic deformation are wavelength-dependent, and at long wavelength these times are ≳10–20 times the Maxwell relaxation time of the dominant ductile component, depending on the details of the viscoelastic model. If Maxwell relaxation times of 10–40 yr apply to a region and recurrence intervals are ∼1000 yr, then a single dislocation contributes a signal to the velocity field that varies only moderately with time. This justifies the use of a statistical measure of deformation, and I define the expected interseismic velocity field to describe the response of an ensemble of faults distributed over a large area and rupturing during a large number of past slip events. Assuming continual past rupture of individual faults, each with a definite mean recurrence time, I derive a simple relationship between the EIV field and the rate of moment release in the lithosphere.

The synthetic examples presented here demonstrate the following.

• Velocity gradients similar to those observed in continental interiors can be produced on a viscoelastic model with 100s of faults distributed over ∼10⁶ km² with a mean recurrence interval of 1500 yr.

• The EIV formulation allows reconstruction of time-averaged moment release rate provided that the Maxwell relaxation time is sufficiently long. In the considered cases, η_m≳ 6 × 10¹⁸ Pa s, which implies a Maxwell relaxation time of ≳13 yr, approximately 1 per cent of the mean recurrence interval assigned to the faults.

• The instantaneous velocity field may exhibit velocity gradients even in localities where dislocation sources are absent.

• There is an important distinction between the interseismic velocity field and the secular velocity field. In particular, a model of secular velocity does not allow satisfactory reconstruction of the average rate of moment release in a region populated by numerous faults.

McKenzie & Jackson (1983) have pointed out that distributed faulting within a deformation zone will be stable only under special conditions (for example, a long, straight strike-slip fault in a shear zone). Such conditions will not generally be met by an unconstrained model of distributed ⁠, i.e. the fault geometry consistent with in a given locality will generally not be constant with time. In the EIV formulation one must therefore assume that the equivalent fault geometry at all localities is approximately constant over the rupture history for which significant relaxation signals are contributed. For the viscoelastic structure believed to apply to the western US, this would mean constant fault geometry on timescales of ∼10⁴ yr.

The relationship provided by eq. (6) is, in part, a kinematic model because it results in estimation of moment release rates. It is also a dynamic model because it is derived from a physical model of crustal deformation. However, this dynamic model is incomplete because deformation sources have not been related to tectonic driving forces. As is the case with purely kinematic models (e.g. Spakman & Nyst 2002), the present model is intended to provide the intermediate step in the process starting with observation of the instantaneous velocity field, then inference of moment release rates, and finally discerning the nature of tectonic forces that drive the long-term deformation pattern.

Acknowledgments

I am grateful to Jim Savage, Marleen Nyst and Wayne Thatcher for their comments and discussion of a preliminary draft of this paper. This paper benefited from constructive comments by two anonymous reviewers and the Associate Editor.

References