Seismological Research Letters
 © 2004 by the Seismological Society of America
Abstract
In order to understand the mechanics of intraplate earthquakes better, a simple 2D numerical model was developed to try to explain current seismicity in the New Madrid seismic zone, using a distinct element method. The model comprises a block geometry representing the structural framework of the New Madrid seismic zone, consisting of intersecting faults with elastic properties corresponding to the known geology. The blocks were subjected to tectonic loading for four days along the direction of the maximum horizontal stress field and the resulting patterns of stress and strain distributions were studied. The results of the modeling showed that shear stresses were higher within the Reelfoot Rift than outside it. In this 2D model the shear stresses on the horizontal plane gave a sense of rotation of the modeled blocks, and an implied sense of movement on the faults. They duplicated the rightlateral strikeslip movement along the Blytheville Fault zone and New Madrid North Fault, and leftlateral strikeslip movement along the Reelfoot Fault. Due to the twodimensional nature, however, results of modeling do not show the observed reverse motion along the Reelfoot Fault. The observed seismicity pattern was consistent with the amplitudes and signs of the maximum shear stresses along the major faults located within the Reelfoot Rift. A linear extrapolation of model results gave an annual strain rate consistent with geodetic observations. The results of modeling support the idea that in a localized volume of preexisting weak crust, fault intersections act as stress concentrators and cause anomalous stress buildup in their vicinity, resulting in observed seismicity.
INTRODUCTION
Large intraplate earthquakes are generally very devastating, although they occur less frequently than their plateboundary counterparts. Their causes are not well understood, as most of them have not been extensively studied because of their rarity, lack of surface expression, and adequate instrumental data. An exception are the earthquakes that collectively define the New Madrid seismic zone (NMSZ), a subject of several recent investigations. These include tectonic and geological studies (Van Arsdale et al., 1995, 1999; Johnston and Schweig, 1996; Van Arsdale et al., 1998; Guccione et al., 2000; Cox et al., 2001), potential field analyses (Hildenbrand, 1985; Rhea and Wheeler, 1994; Hildenbrand and Hendricks, 1995; Hildenbrand et al., 1996; Hildenbrand et al., 2001), seismicity evaluation (Chiu et al., 1992), estimates of the magnitudes of the 18111812 earthquake sequence (Johnston, 1996; Hough et al., 2000), paleoseismological investigations (Schweig and Ellis, 1994; Schweig and Tuttle, 1996; Tuttle et al., 1999; 2002), seismic refraction studies (Ginzburg et al., 1983; Mooney et al., 1983), and geodetic studies (Liu et al., 1992; Snay et al., 1994; Weber et al., 1998; Neumann et al., 1999; Gan and Prescott, 2001). These studies have not only helped define the spatial extent, seismogenic features, and structural framework of the New Madrid seismic zone but have also refined the locations of seismicity, derived the recurrence rates of large earthquakes, and discovered evidence for neotectonic activity.
Following these multidisciplinary studies resulting in extensive and highquality data, our understanding of the nature of seismicity in the New Madrid seismic zone has improved greatly and spawned some early mechanical models to explain the cause of seismicity at NMSZ (Gomberg, 1992; Gomberg and Ellis, 1994). Gomberg and Ellis (1994) employed the boundary element technique with morphological, geological, and seismicity data as constraints. They addressed the underlying driving mechanism that causes deformation in NMSZ and a method for accommodation of that deformation. More recent models by Grana and Richardson (1996), Stuart et al. (1997), Kenner and Segall (2000), Thybo et al. (2000), and Pollitz et al. (2001) address stress concentration, observed recurrence rates, surface deformation, and an apparent absence of large strain accumulation in NMSZ. They assumed that the earthquake generation cycle was initialized by perturbations in the lower crust and upper mantle, which then loaded the seismogenic structures in the upper brittle crust and triggered seismicity.
We present an alternative, simple, testable model wherein the seismicity occurs on faults with uniform frictional properties in response to stress perturbations in the strong brittle upper crust because of plate tectonic forces. Our model does not attempt to explain the initiation, history, genesis, or temporal nature of seismicity or deformation in the New Madrid seismic zone. Instead, it is limited in scope and attempts to explain the location and nature of the current seismicity, with the buildup of stresses between large earthquakes whose return periods are defined by paleoseismology. In this sitespecific model for NMSZ we have incorporated structural elements based on observations in the upper crust where most of the seismicity is located, and utilized plate tectonic stresses as driving forces. Although a 3D model is most appropriate, we tested our ideas with a preliminary 2D model based on the intersection model (Talwani, 1988).
To explain the genesis of intraplate earthquakes, based on an evaluation of five case studies, Talwani (1988) suggested that stresses can accumulate in the vicinity of intersecting tectonic features in response to an ambient stress field. His observations were supported by results from a 2D analysis of an intraplate earthquake source (Jing and Stephansson, 1990) and analytical computations by Andrews (1994). Recently, Talwani (1999) modified the model of Jing and Stephansson (1990) and illustrated the role of fault geometry with examples from three intraplate regions. In a more recent study, Dentith and Featherstone (2003) showed that essential elements of the 1968 Meckering earthquake and observed seismicity in southwestern Australia can be explained by Talwani's “intersection model.” Talwani and Gangopadhyay (2000) and Gangopadhyay and Talwani (2003) examined 39 case histories of earthquakes in twenty intraplate regions to seek common features between them. They found that within preexisting weak zones, stresses can build up in response to plate tectonic forces in the vicinity of fault intersections, plutons, and/or buried rift pillows. They inferred from these observations a possible causal mechanism for intraplate earthquakes. In this paper we test these ideas by developing a mechanical model for NMSZ using a simple 2D distinct element code called UDEC. We divided the upper crust into blocks subjected to plate tectonic stresses. We show that with appropriate fault geometry the stress buildup in the vicinity of fault intersections can explain the observed seismicity.
2D DISTINCT ELEMENT MODELING USING UDEC
Following Jing and Stephansson (1990), numerical modeling of stress accumulation in faulted rock masses was carried out by using the distinct element method, which models the rock mass as an assemblage of discrete blocks (rigid or deformable) and the joints and faults as discontinuities. The equations of motion for the blocks are solved by a central difference scheme, and mutual interactions between blocks are included. The large deformations, including block rotations, are also accounted for in this method.
The distinct element method was used to model NMSZ assuming 2D geometry and using a program called Universal Distinct Element Code (UDEC) written by Itasca Consulting Group, Inc., Minneapolis, Minnesota (Version 3.1, 1999). This 2D numerical program was first developed by Cundall (1971). The program simulates the response of the discontinuous media to either static or dynamic loading. It uses calculations in the Lagrangian scheme to model large movements and deformations of a system. Several builtin material behavior models, for both the intact blocks and the discontinuities/faults, permit the simulation of real geologic situations. Large displacements are allowed along the faults, which are treated as boundaries between blocks, allowing the blocks to move with respect to each other. The individual blocks can be made either rigid or deformable. The deformable blocks are divided into a mesh of triangular constantstrain finitedifference zones, and each zone behaves according to a prescribed stressstrain law. In the case of elastic analysis the formulation of these zones is identical to that of constantstrain finite elements. The relative motions along the discontinuities are constrained by forcedisplacement relations for movement in both the normal and shear directions.
The suitability, efficiency, and adaptability of UDEC in solving geological problems involving faulted and fractured rocks has been adequately demonstrated and established. UDEC has been successfully used to model the response of fractured rock to fluid injection (Last and Harper, 1990), effects of fault slip on fluid flow around extensional faults (Zhang and Sanderson, 1996), critical behavior of deformation and permeability of fractured rock masses (Zhang and Sanderson, 1998), effects of fracture geometry and loading direction on stability of fractured rock masses (Zhang and Sanderson, 2001), and fluid flow in deforming fractured rocks (Zhang et al., 2002). Simulations of fault slip and stress perturbations using UDEC were also performed in studies representing various tectonic settings (Dupin et al., 1993; Homberg et al., 1997). UDEC was also used to obtain a numerical model of an intraplate earthquake source and to conduct parametric studies of the effect of a fault on in situ stresses (Jing and Stephansson, 1990; Su and Stephansson, 1999).
DEVELOPMENT OF THE 2D MODEL
The first step toward modeling NMSZ was to represent the model geometry and associated parameters that mimic the geologic framework. These are described below.
Model Geometry
Figure 1 shows the structural framework of NMSZ as outlined by Hildenbrand et al. (2001) and the instrumentally located earthquakes with M 3.0 or greater adapted from the 19742002 CERI, Memphis catalog. The essential pattern of seismicity down to magnitude 1.5 (Figure 1 in Hildenbrand et al., 2001 is the same. We chose a lower cutoff at M 3.0 (Figure 1) for clarity and ease in comparing with stress calculations. In NMSZ, within the northeastsouthwesttrending, nearly 400kmlong and 100kmwide Reelfoot Rift, there are two intersecting fault zones, the ∼65kmlong Blytheville Fault zone (BFZ) oriented northeastsouthwest and the ∼60kmlong Reelfoot Fault zone (RF) oriented northwestsoutheast (Figure 1; Van Arsdale et al., 1995; Johnston and Schweig, 1996). A third fault, the New Madrid North Fault (NMNF), lies outside the edge of the floor of the Reelfoot Rift but within its edge (Figure 1; Rhea and Wheeler, 1995). This ∼30kmlong northnortheasttrending fault is considered to be the extension of the aseismic Bootheel lineament (BL) (Johnston and Schweig, 1996). The Missouri batholith, a prominent, 100kmwide, lowdensity, granitic, upper crustal structure, trends northwestsoutheast and crosscuts the Reelfoot Rift (Hildenbrand et al., 1996; 2001). The major plutons in the New Madrid seismic zone (Figure 1), including the Bloomfield and Covington plutons and Osceola Igneous Complex (OIC), are at a depth of about 210 km and are mostly mafic in nature (Hildenbrand et al., 2001).
The observed seismicity inside the Reelfoot Rift is located along the Blytheville Fault zone, Reelfoot Fault, and New Madrid North Fault (Figure 1). In particular there is clustering of seismicity at and/or near the intersections of the Bootheel lineament, Blytheville Fault zone, and Missouri batholith; Reelfoot Fault and Blytheville Fault zone; and Reelfoot Fault and New Madrid North Fault (Figure 1). A few earthquakes are also located along the southern margin of the rift (Figure 1).
The block geometry used for modeling NMSZ was modified from an outline of the structures in Hildenbrand et al. (2001). The rift boundary given by Hildenbrand et al (2001) is based on aeromagnetic data and represents the floor of the Reelfoot Rift. The western boundary of the Reelfoot Rift was taken from Rhea and Wheeler (1995). The eastern boundary of the Reelfoot Rift is not well defined, and so in our model it is based on the outline of the rift floor given by Hildenbrand et al. (2001). The block model used in our calculations is shown in Figure 2A. It is strictly twodimensional. In order to study the effect of the fault intersections alone, the plutons have been excluded. The three New Madrid faults, the Reelfoot Rift margins, the Bootheel lineament, and the margins of the Missouri batholith within the rift have been incorporated. The outermost edges of the block measure 220 × 220 km (Figure 2A). These dimensions were chosen to accommodate the most seismically active part of the Reelfoot Rift. The direction of the maximum horizontal stress, S_{Hmax}, is N80°E for NMSZ (Zoback, 1992). For computational convenience of modeling, the block boundaries including the faults, Reelfoot Rift margins, and margins of the Missouri batholith were rotated by 10° clockwise so that S_{Hmax} lies along the x axis (Figure 2A). The block corners have been rounded with a circle that is tangential to the two corresponding edges at a specified rounding distance from the corner. In practice, the rounding distance is about 1% of the typical block edge length (Itasca Consulting Group, 1999). The rounding length for the block corners in our model is 2 km. Since our model is twodimensional we have imposed the commonly used plane stress condition, wherein none of the blocks experiences stress in the vertical direction, although they can exhibit strain in that direction.
The code divides the deformable blocks into triangular finitedifference zones using a builtin automatic mesh generator that decides the size of the elements based on the block lengths, specified rounding lengths, and memory availability to perform the computations. The mesh representation of our model for NMSZ is shown in Figure 2B. All the blocks in the model are deformable and movable with respect to each other. UDEC calculates the shear stress (τ_{xy}) at each node. Figure 2C shows the sign convention. Shear stress (τ_{xy}) is positive (Figure 2C) when it tends to rotate the block in a counterclockwise manner, i.e., by leftlateral strikeslip, and negative when the block rotation is clockwise. The code calculates the amplitude and sign of the shear stress at each node to show how the block will rotate and contours to show its spatial variation.
Model Parameters
UDEC has seven builtin constitutive models for the blocks and four for the joints that can represent various geologic situations. We have used the simplest constitutive models for this preliminary twodimensional model. In our model the blocks conform to the Linearly Elastic Isotropic Model and the joints follow the Joint Area Contact Elastic/Plastic Coulomb Slip Failure Model. The Linearly Elastic Isotropic Model for the blocks describes the simplest form of material behavior, assuming homogenous and isotropic materials that exhibit linear stressstrain behavior with reversible deformation upon unloading (Itasca Consulting Group, 1999). The Joint Area Contact Elastic/Plastic Coulomb Slip Failure Model for the joints is the most commonly used Coulomb slip model that predicts failure or initiation of slip on a fault based on the accumulated shear stress (Itasca Consulting Group, 1999). It is represented by the following equations: (1) (2) where Δσ_{n}, Δμ_{n}, Δτ_{S} and Δμ_{S} are the effective normal stress, normal displacement, shear stress, and shear displacement increments respectively.
The failure criterion for the joints is given by (3) where τ_{S} is shear stress, C is cohesive strength of the joint, σ_{n} is normal stress, and φ is friction angle for the joint.
The block assembly was subjected to a horizontal force along the x axis. The value of this compressive force was derived from the plate velocity measured from geodetic observations. In geodetic studies in the New Madrid seismic zone, Liu et al. (1992) obtained a velocity of 57 mm/year across NMSZ and Weber et al. (1998) obtained about 35 mm/year. From their analysis of GPS data, however, Dixon et al. (1996) and Newman et al. (1999) argue that the differential velocity at NMSZ is not resolvable from overall North America motion within appreciable confidence limits, whereas Gan and Prescott (2001) obtained a southward velocity of 1.7 ± 0.9 mm/year for the entire Mississippi embayment. In our model we have applied an eastwest velocity of 5 mm/year to provide the driving force. We have assumed that the velocity gradient is not a function of depth and hence the whole block can be subjected to the same horizontal stress field. The behavior of the blocks is dependent on the value of the applied stress field. In our case we have chosen the larger value in order to obtain measurable response with a shorter loading time used in the model. The calculated stresses would scale down linearly if we used a velocity of ∼2 mm/year.
Model Properties
Input variables for the model calculations include elastic moduli and density for the blocks. To estimate these we utilized values of the P and Swave velocities, V_{P}, and V_{S}, the velocity ratio, V_{P}/V_{S}, and density obtained from independent geophysical studies and gravity modeling, respectively, and derived the corresponding values of the moduli using the formulae (4) (5) where μ is shear modulus, k is bulk modulus, and ρ is density.
Using nonlinear inversion and traveltime tomography, Vlahovic and Powell (2001) obtained V_{p}/V_{S} ratio values along a ∼60kmlong profile almost parallel to the Reelfoot Fault, from near its intersection with the northwestern rift boundary, southeast to near the rift axis. They found the V_{P}/V_{S} ratio value to be 1.73 for at least half the profile except for a part to the southeast of its intersection with the northeasttrending Blytheville Fault zone, where it was 1.8. In our computations of the elastic moduli to be used in model calculations, we have assumed the more “normal” value of 1.73 for the V_{P}/V_{S} ratio, for rocks both inside and outside the Reelfoot Rift.
Based on extensive seismic refraction studies in the Mississippi embayment, including the Reelfoot Rift and immediately outside it, Mooney et al. (1983) obtained V_{P} values and developed a crustal velocity model. Their V_{p}—depth profile across the Reelfoot Rift showed different layers of varying Pwave velocities and thicknesses. The top 5 km of the crust inside the rift was divided into layers of 1.5 km, 2.5 km, and 1 km thickness with Pwave velocities of 1.8 km/s, 5.95 km/s, and 4.9 km/s, respectively. The crust from a depth of 5 to 15 km inside and outside the rift showed a generally uniform Pwave velocity of 6.2 km/s. The weighted average of the Pwave velocity for the crust inside the Reelfoot Rift was found to be 5.63 km/s. In a highresolution seismic refraction/reflection survey in the region of most intense seismicity inside the rift, Williams and Catchings (1992) also obtained a Pwave velocity of 5.45 km/s below a depth of 780 m. The crustal density model of Mooney et al. (1983) comprises three layers of different thicknesses and densities within the Reelfoot Rift. The weighted average density inside the Reelfoot Rift was 2,690 kg/m^{3}. Their model suggests a density of 2,750 kg/m^{3} outside the rift. We utilize these values for blocks within and outside of the rift respectively. In calculating the elastic properties of the blocks we used the weighted averages of Pwave velocity (5.63 km/s) and density (2,690 kg/m^{3}) within the rift and 6.2 km/s and 2,750 kg/m^{3} outside the rift.
Based on gravity modeling, the Missouri batholith was located between depths of 510 km inside the Reelfoot Rift (Hildenbrand et al., 2001). The Pwave velocity obtained by Mooney et al. (1983) for this depth range is 6.2 km/s. This value was utilized for our elastic moduli computations for the Missouri batholith. The average density of the batholith was taken to be 2,705 kg/m^{3}, as inferred from gravity modeling studies of Hildenbrand et al. (2001). Thus utilizing these values of V_{P}, V_{P}/V_{S}, and density we computed the bulk and shear moduli for the various blocks (Table 1).
Other input variables required for modeling the deformation include friction angle, normal and shear stiffnesses, and cohesion of the faults (treated as joints). The values of these variables that we chose are given in Table 2. The basis for choosing these values is described next. The seismogenic parts of the Blytheville Fault zone, Reelfoot Fault, New Madrid North Fault, and rift boundary faults are assumed to lie in metamorphic, gneissic crystalline basement (Johnston and Schweig, 1996). The inference of a crystalline basement was based on borehole data reported by Hildenbrand et al. (2001), where granitic gneiss was encountered with intervals of dioritic gneiss. Based on seismic reflection data in the southeastern Reelfoot Rift margin region, Crone (1992) concluded that the rift boundary faults extend into the crystalline basement. Barton (1976) tabulated the basic friction angles for various rock types based on experimental laboratory results, estimating a friction angle for gneiss ranging between 25°29° and for granite between 31°35°. We assigned friction angles of 27° to all the faults based on the mean of the range given for gneiss by Barton (1976) (Table 2). The granitic Missouri batholith (Hildenbrand et al., 2001) has been assigned a friction angle of 33° based on the mean of the range given for granite by Barton (1976) (Table 2).
Rosso (1976) compared joint stiffness laboratory measurements and previously published results of tests on jointed samples of quartz diorite. He obtained joint normal and shear stiffnesses of 101 GPa/m and 76 GPa/m, respectively, at applied stresses of 10.5 MPa. We have used these stiffness values for the rift boundary faults and the three major faults in NMSZ (Table 2). Likewise, Rosso (1976) reported jointed triaxial test results by Brown and Swanson (1972) on samples of granite which gave joint normal and shear stiffnesses of 133 GPa/m and 100 GPa/m respectively. We have adopted these values for the granitic Missouri batholith (Table 2). Since the stiffness values vary for different magnitudes of applied stress levels, we experimented with values that were higher (two and five times), significantly higher (ten times), lower (half and onefifth), and significantly lower (onetenth) than the aforesaid values. The model results showed no qualitative change in the pattern of stresses and did not affect the overall conclusions of this paper.
McGarr and Gay (1978) concluded from different types of stress measurements made in southern Africa, North America, and Australia that the lower limits of maximum shear stress, (S_{1}  S_{3})/2, at midcrustal depths are 20 to 40 MPa. S_{1} and S_{3} are the maximum and minimum principal stresses respectively. Based on measurements made in the KTB borehole in Germany and assuming a strikeslip regime and coefficient of friction of 0.7, Zoback et al. (1993) suggested that at midcrustal depths the differential stress could reach values of ∼300 MPa, implying that the upper bound for maximum shear stress at those depths could be ∼150 MPa. Furthermore, Zoback et al. (1993) also suggested that similar conditions may exist in the seismically active part of eastern North America. Based on these studies, we assumed that the maximum shear stress in NMSZ in the brittle crust from the surface up to midcrustal depths lies between 20 MPa and 150 MPa. We assume the joint cohesion for the three major seismogenic faults in NMSZ to be zero at these stress levels (Byerlee, 1978; Homberg et al., 1997). For the aseismic rift boundary faults and the margins of the Missouri batholith we arbitrarily assign a joint friction of 0.5 MPa (Table 2).
Limitations of the Model
A notable limitation of this model is the fact that it is twodimensional. Artifact boundary effects or the edge effects also are present in the model. These manifest themselves as comparably high stress and strain values near the outermost boundaries of the blocks where the northern and southern rift margins intersect (Figures 3 and 4). The educational version of UDEC has limited memory that prohibits extensive computations and thus our ability to run the model for a geologically realistic loading time.
To study the effect of the intersections the model has been simplified by not including the plutons. The memory constraint also affects the run time of the model. The model was run for tectonic loading times corresponding to one, two, and four days. The calculated stress and strain buildup were linear with time and do not negate the conclusions of the model. In this case the longest run was for time corresponding to tectonic loading of four days.
MODEL RESULTS AND THEIR ANALYSIS
The modeling outputs were analyzed in terms of the resulting stresses, strains, and rotations in response to a tectonic loading of four days. These are described below. The NMSZ model geometry has been superimposed for convenience of comparison. The magnitudes of the stresses, strains, and rotations represent relative values of these variables at different locations. They show the pattern of stress concentration and were used to make tectonic interpretations.
Shear Stress
Shear stresses (τ_{xy}) were obtained at each node of the mesh (Figure 2B) using the sign convention shown in Figure 2C. Recall that this is a twodimensional model and that positive and negative values suggest counterclockwise and clockwise rotation (for the block) respectively. The stress values at the nodes were contoured with a contour interval of 0.25 N/m^{2} (Figure 3). Shear stresses (1 to 1.75 N/m^{2}) seen at locations near the edges of blocks are artifacts of boundary effects in the calculations and are ignored. Away from these edges the shear stress is greater in the preexisting weak zone inside the rift (0 to 2.5 N/m^{2}) compared to the stronger zone outside it (1.0 to 0 N/m^{2}). The shear stress values range from 0 ± 0.5 N/m^{2} along the rift boundaries to about 1.75 N/m^{2} and 2.5 N/m^{2} near the intersections of the Blytheville Fault zone and Bootheel lineament with the Missouri batholith, and of the Blytheville Fault zone and Reelfoot Fault respectively. These larger stress values are concentrated in very small regions near the intersections. At these locations, the shear stresses have the greatest potential for loading the faults. The area in the immediate vicinity of the intersection of the Reelfoot Fault, Bootheel lineament, Missouri batholith, and New Madrid North Fault showed stresses of ∼ 1.0 N/m^{2} (Figure 3).
Shear Strain
Figure 4 shows the corresponding shear strain with a contour interval of 5 × 10^{12}. The sign convention in this case is similar to that for shear stress. Shear strains (2 × 10^{11} to 3 × 10^{ll}) seen near the edges of blocks are artifacts of boundary effects in the calculations and are ignored. The shear strain is in general greater inside the rift (0 to 3.5 × 10^{11}) compared to that outside it (1.5 × 10^{11} to 0). The shear strain values range from 0 ± 1 × 10^{11} along the rift boundaries to about 2.5 × 10^{11} and 3.5 × 10^{11} near the intersections of the Blytheville Fault zone and Bootheel lineament with the Missouri batholith, and of the Blytheville Fault zone and Reelfoot Fault, respectively. These larger strain values are concentrated in very narrow regions within the rift, and near the intersections, and thus are locations with the greatest potential for seismicity. The area in the immediate vicinity of the intersection of the Reelfoot Fault, Bootheel lineament, and New Madrid North Fault showed strains of ∼ 1.5 × 10^{11} (Figure 4). These are lower than at the other intersections, suggesting a lesser tendency for block rotation (and strikeslip faulting ??).
Block Rotations
Different shear stresses acting on the block boundaries will tend to cause strikeslip movements. If the movement of one block is obstructed by another, the result will be rotation and uplift. Our model being twodimensional, we do not observe uplift but do note some evidence of block rotation. The direction of rotation of each block is shown in Figure 5. The maximum rotation obtained was 9.055 × 10^{10} degrees clockwise for the block immediately adjacent to the Reelfoot Fault to the east. The rotation “arc” is scaled to a maximum arc of 45°. The block rotations shown with respect to the perpendicular to the S_{Hmax} direction are intended to provide a direction of block movement only. The results indicate that the blocks outside of the rift show counterclockwise rotation, whereas the blocks inside the rift show both clockwise and counterclockwise rotation. Inside the rift the counterclockwise rotation of the block south of the Blytheville Fault zone implies rightlateral strikeslip motion along the fault. Similarly, the counterclockwise rotation of the blocks southeast of the New Madrid North Fault also imply rightlateral strikeslip motion along the fault. The directions of block movement are consistent with shear stress directions. Furthermore, the counterclockwise rotation of the blocks north of the New Madrid North Fault are consistent with the counterclockwise offset of that fault as observed in other independent geological studies (Johnston and Schweig, 1996).
Maximum Shear Stress Plots along the Faults
The concentration of shear stresses near the intersections suggested that these parameters were not uniform along the faults. To study the variation of these parameters, we plotted the maximum shear stress along the three major fault planes in NMSZ and Bootheel lineament at the end of a tectonic loading period of four days (Figure 6, curves ad). The maximum shear stress along the plane of the Blytheville Fault zone (curve a) between its intersection with the Missouri batholith (A) and Reelfoot Fault (B) is nearly constant and lies between 6.6 N/m^{2} and 7 N/m^{2}. The negative values of maximum shear stress indicate rightlateral movement along the fault. The largest absolute values on this profile are observed at the intersections of the Blytheville Fault zone, Bootheel lineament, and Missouri batholith (6.9 N/m^{2}) (A), and at the intersection of the Blytheville Fault zone and Reelfoot Fault (7 N/m^{2}) (B).
Curve b shows the variation of maximum shear stress along the plane of Reelfoot Fault between its intersection with the southern rift margin (P) and New Madrid North Fault (Q). The positive values along the Reelfoot Fault plane indicate leftlateral strikeslip movement along the fault. The maximum shear stress increases steadily from the ends to its maximum value at B (7 N/m^{2}), which is the intersection of the Reelfoot Fault and Blytheville Fault zone. It reduces to a base level of about 00.5 N/m^{2} and rises again to about 6 N/m^{2} at the intersection of the Reelfoot Fault and New Madrid North Fault.
Similarly, curve c shows the variation of maximum shear stress along the plane of the New Madrid North Fault between its intersection with the Reelfoot Fault (Q) and northern edge of the Reelfoot Rift (Y). The negative values along the New Madrid North Fault plane indicate rightlateral movement along the fault. The absolute value of maximum shear stress steadily rises from 6 N/m^{2} at its intersection with the Reelfoot Fault (Q) to 11 N/m^{2} at its intersection with the Missouri batholith (N), and then gradually reduces to about 8 N/m^{2} further away from this intersection (Y).
Curve d shows the maximum shear stress on a plane lying along the Bootheel lineament. The negative values along the Bootheel lineament indicate rightlateral movement along the fault. The largest absolute values are at the ends, near the intersection with the Blytheville Fault zone (A) (6.3 N/m^{2}) and the intersection with the New Madrid North and Reelfoot Faults (Q) (7 N/m^{2}). We discuss the tectonic interpretation of these results next.
DISCUSSION
Figures 3, 4, 5, 6 show the results of our modeling. Even though the shear moduli of the blocks inside and outside the Reelfoot Rift are different (Table 1), we notice remarkable similarities between the spatial distribution of shear stresses (Figure 3) and shear strains (Figure 4). In both cases we note that the interior of the rift has larger values. The shear stress contours show both positive and negative values, which are indicative of leftlateral and rightlateral movements, respectively. The sense of block movement (Figure 5) also confirms the inferred sense of movement along the faults. Figure 6 shows how the maximum shear stress varies along the fault planes.
Accurate seismicity data are available for NMSZ for magnitudes of 1.5 and greater. We compared the spatial pattern of seismicity for magnitude 1.5 and greater (e.g., in Figure 1 in Hildenbrand et al., 2001) with that for magnitudes 3.0 and greater (Figure 1). The seismicity pattern is similar in the two maps, continuous along the Blytheville Fault zone (BFZ) (AB) and near the ends of the Bootheel lineament where it intersects with BFZ and the Reelfoot Fault. Along the Reelfoot Fault we find intense seismicity from near its intersection with BFZ (point B), southeast 30 km to location M. There is a hiatus of intense seismicity for about 10 km northeast of B on the two maps and intense activity near its intersection with the New Madrid North Fault. In short, Figure 1 (for M 3.0 and greater) is a less cluttered but true representation of the seismicity pattern in NMSZ. We use it to compare with the maximum shear stress along the faults. We make the tacit assumption that the seismicity pattern is representative of the distribution of stresses on the faults. Comparing the stress profiles along the faults (Figure 6) with seismicity (Figure 1) we observe excellent correspondence. The absolute values of the shear stresses depend on the model parameters, and those obtained in Figure 6 are representative, but their relative values are more instructive.

The negative maximum shear stress values along the Blytheville Fault zone, New Madrid North Fault, and Bootheel lineament (Figure 6, curves a, c, and d) indicate rightlateral strikeslip movement, in agreement with the seismicity observations (Johnston and Schweig, 1996). The Bootheel lineament has the lowest absolute values, suggesting that faulting is more likely to occur on the other two in the presentday stress field. The positive maximum shear stress values along the Reelfoot Fault suggest leftlateral strikeslip motion. Faultplane solutions for three of the larger instrumentally recorded earthquakes along this fault (3 March 1963, 21 July 1967, and 4 May 1991) showed leftlateral strikeslip movements (Herrmann and Ammon, 1997). The dominant motion on the Reelfoot Fault is reverse, which this 2D model is incapable of replicating. These motions on the Blytheville Fault zone, Reelfoot Fault, and New Madrid North Fault are also indicated by the directions of block rotations (Figure 5).

The magnitudes of maximum shear stresses are not uniform along the faults. They are largest in the vicinity of intersections with other faults. Interestingly, when the seismicity (Figure 1) is compared with the stress profiles (Figure 6), we find that it occurs at locations along the faults where the absolute stress value is 6 N/m^{2} or greater. This value is a consequence of the modeling parameters, and it indicates a stress threshold for the onset of seismicity. In Figure 6 we also compare the stress changes along the fault with respect to this threshold value. For the Blytheville Fault zone the largest shear stresses are at A and B (Figure 6, curve a), the intersecting points with the Missouri batholith and Bootheel lineament; and with Reelfoot Fault respectively. Seismicity occurs all along the length of the Blytheville Fault zone, where the absolute value of the maximum shear stress is uniform at about 6.6 to 6.7 N/m^{2}, slightly greater than the threshold value. Along the Reelfoot Fault, the maximum shear stresses exceed 6 N/m^{2} for about a ∼30km (BM) section near and to the southeast of its intersection with the Blytheville Fault zone (Figure 6, curve b), the location of intense seismicity (Figure 1). Along the New Madrid North Fault the absolute value of maximum shear stress is 6 N/m^{2} at its intersection with the Reelfoot Fault (Q) and builds to the largest absolute value of about 11 N/m^{2} at its intersection with the Missouri batholith (N) (Figure 6, point N on curve c), considerably greater than the threshold value, and the location of intense seismicity (Figure 1). Along the Bootheel lineament the seismicity occurs near its ends, the regions of its intersection with BFZ and RF, where the maximum shear stress is greater than the threshold value.

To see if our calculated stresses were physically reasonable, we computed the strain rates predicted by the model and compared them with the observed geodetic results in NMSZ.
Liu et al. (1992) obtained a shear strain rate of 0.108 ± 0.045 μradian/year for the southern part of NMSZ. Snay et al. (1994) obtained a value of 0.030 ± 0.019 μradian/year for the northern part of NMSZ. Comparing results of two GPS surveys covering almost the entire NMSZ, Weber et al. (1998) divided it into two subsets, to the north and to the south. He obtained strain rates of 0.121 ± 0.083 μradian/year and 0.080 ± 0.058 μradian/year for the two subsets. Contradicting all of these measurements, however, Newman et al. (1999) argued that the NMSZ region has almost no strain as a whole when viewed along a 600kmlong line.
Figure 1 shows the locations of three points (A, B, and Q) where we computed the strain rates from our model results. The study area of Liu et al. (1992) is immediately north of A. The two other locations (B and Q) lie within the study areas of Snay et al. (1994) and Weber et al. (1998). To compare the results of our model we calculated the annual strain rates at these points from the maximum shear stress values obtained from the model at those points (Figure 6). The absolute values of the maximum shear stresses at A, B, and Q are 6.9, 7, and 6 N/m^{2} (Pa), respectively, for a tectonic loading cycle of four days. Based on shear modulus of the rocks at these points (28.48 GPa, Table 1), the annual strain rates at A, B, and Q are 0.044, 0.045, and 0.038 μradian/year respectively. These values are in general agreement with the results of the geodetic surveys.

The shear strain obtained from our model results at B is 0.045 μradian/year (engineering strain), which corresponds to a linear strain of ∼2.24 × 10^{8}/year. If extrapolated linearly for the paleoseismic recurrence time of 500 years (Tuttle et al., 2002), this generated a total strain of ∼1.12 × 10^{5}, which is a reasonable value for the generation of earthquakes (Bilham et al., 1989).
CONCLUSIONS
To explain the current seismicity in the New Madrid seismic zone we developed a simple 2D model wherein the fault geometry was represented by blocks. These were subjected to tectonic loading for four days along the direction of the maximum horizontal stress, and the resultant patterns of the shear stresses, shear strains, and block rotations were examined.
The results obtained from this simple model could explain many firstorder observations at NMSZ. The shear stresses and shear strains that developed after four days of tectonic loading were larger within the Reelfoot Rift compared to the regions outside it. The largest shear stresses and strains were obtained at and near the fault intersections, and they indicated the tendency of rotation of the blocks and the sense of movement along the faults. The model results duplicated the sense of movements along the Blytheville Fault zone, Reelfoot Fault, and New Madrid North Fault. Because our model is twodimensional, reverse motion along the Reelfoot Fault was not observed. The model did indicate leftlateral strikeslip movement, however, consistent with the faultplane solutions of three instrumentally recorded earthquakes along the Reelfoot Fault (Herrmann and Ammon, 1997). The distribution and sign of maximum shear stress along the faults were consistent with observed seismicity. The strain rate obtained from linear extrapolation of the model results was in general agreement with that observed from geodetic surveys. Overall, the model results support the idea that fault intersections within the preexisting weak crust can be the focus of stress accumulation when subjected to plate tectonic forces, and that these intersections then become locations of strain buildup, causing earthquakes.
Acknowledgments
The comments of an anonymous reviewer caused us to reflect on and clarify some parts of our paper. We also thank Drs. Lanbo Liu and Martin Chapman for their helpful comments. Westinghouse Savannah River Company is thanked for providing partial support for A.G.
Footnotes
Department of Geological Sciences, University of South Carolina
Department of Civil and Environmental Engineering, University of South Carolina
Department of Geological Sciences, University of South Carolina