Paul Spudich (U.S. Geological Survey, MS977, Menlo Park, CA 94025 USA, spudich@usgs.gov)

Mariagiovanna Guatteri (Dept of Geophysics, Stanford Univ., Stanford, CA 94305 USA, patti@pangea.stanford.edu)

Kenshiro Otsuki (Dept of Geoenvironmental Sciences, Tohoku Univ., Sendai 980, JAPAN, otsuki@dges.tohoku.ac.jp)

Jun Minagawa (Dia Consultants Co. Ltd., Ikebukuro 3-1-2, Toshima-ku, Tokyo 171, JAPAN, J.Minagawa@gen.diaconsult.co.jp)

Dislocation models of the 1995 Hyogo-ken Nanbu (Kobe) earthquake derived by Yoshida et al. (1996) show substantial changes in direction of slip with time at specific points on the Nojima and Rokko fault systems, as do striations we observed on exposures of the Nojima fault surface on Awaji Island. Spudich (1992) showed that the initial stress, i.e. the shear traction on the fault before the earthquake origin time, can be derived at points on the fault where the slip rake rotates with time if slip velocity and stress change are known at these points. From Yoshida's slip model we calculated dynamic stress changes on the ruptured fault surfaces. To estimate errors we compared the slip velocities and dynamic stress changes of several published models of the earthquake. The differences between these models had an exponential distribution, not gaussian. We developed a Bayesian method to estimate the probability density function of initial stress from the striations and from Yoshida's slip model. Striations near Toshima and Hirabayashi give initial stresses of about 13 and 7 Mpa, respectively. We obtained initial stresses of about 7 to 17 Mpa at depths of 2 to 10 km on a subset of points on the Nojima and Rokko fault systems. Our initial stresses and coseismic stress changes agree well with post-earthquake stresses measured by hydrofracturing in deep boreholes near Hirabayashi and Ogura on Awaji Island. Our results indicate that the Nojima fault slipped at very low shear stress, and fractional stress drop was complete near the surface and about 32% below depths of 2 km. Our results at depth depend on the accuracy of the rake rotations in Yoshida's model, which are probably correct on the Nojima fault but debatable on the Rokko fault. Our results imply that curved or crosscutting fault striations can be formed in a single earthquake, contradicting a common assumption of structural geology.

Most structural geologists and seismologists have an intuitive picture of earthquakes that is implicitly a high stress picture, even though a growing body of evidence from heat flow, stress orientations, and aftershock mechanism diversity suggests that many earthquakes occur at low applied shear stress, and stress drop may be complete (Lachenbruch and Sass, 1980; Zoback et al., 1987; Michael et al., 1990; Beroza and Zoback, 1993; Zoback and Beroza, 1993; Barton and Zoback, 1994). (Here we use the term "stress" to mean the component of shear traction acting tangent to the fault surface.) One of the crucial differences between the dynamics of high and low stress earthquakes is the degree to which slip rake can rotate temporally at a single point on the fault during the earthquake. Dynamic rupture models by Andrews (1994) and Guatteri and Spudich (1998) show that a fault rupturing at low stress can have slip rakes that vary both temporally during the rupture process and spatially, allowing such observable phenomena as curved or crosscutting fault striations and strong spatial variations of slip direction in large earthquakes. These low stress phenomena provide an additional possible interpretation for observed crosscutting fault striations, which are often interpreted to signify multiple earthquakes or tectonic events (e.g. Etchecopar et al. , 1981; Cashman and Ellis, 1994).

If the slip rake rotates with time at a point on a fault during dynamic rupture, then it is possible in principle to learn the absolute level of shear stress acting on that point before the earthquake initiated (Spudich, 1992). This assertion contradicts the widespread seismological lore that it is possible to learn only about stress changes from seismic data, not about the absolute level of shear stress (Aki and Richards, 1980, p. 56). This lore results from studies of slip on 1D faults, on 2D faults in which the slip rake was not allowed to rotate with time (e.g. Madariaga, 1976, who noted that such rotations were possible), and on 2D faults in which temporal rotation was permitted but which were investigated at high stress levels which minimized rotations (e.g. Das, 1981; Day, 1982).

The 1995 Hyogo-ken Nanbu (Kobe) earthquake showed two types of evidence of slip at low applied shear stress. First, curved or crosscutting fault striations were observed at several locations along the Nojima fault on Awaji Island (Otsuki et al., 1997) . Second, a number of groups using ground motion recordings and geodetic data (Ide et al., 1996; Wald, 1996; Yoshida et al., 1996) have found evidence for temporal rotations of slip rake at depth on the Nojima and Rokko fault systems.

One of the purposes of this paper is to improve the methodology of ground motion inversions by emphasizing the significance of temporal rake rotation observations. There is a problem in the inversion procedures used by most past investigators (summarized in Spudich, 1992) and by those cited above. Our current inversion procedures may be giving an overly optimistic picture of our ability to simulate ground motions for engineering purposes.

Typically in these inversion procedures the investigators allow the strike-slip and dip-slip components of slip velocity at each point on the fault to have independent time functions, meaning that slip rake can (and does) rotate temporally in their derived slip models. Although this rotation has never hitherto been regarded as significant, it doubles the number of free parameters in the slip model, and it allows the ground motions to be simulated much more accurately than would occur if no temporal rotation were allowed. If such accurate simulations of ground motion are illusory, then the ground motion community may be overestimating the effectiveness of ground motion simulations for engineering purposes.

The second purpose of this paper is to learn what the observed temporal rake rotations imply about the preseismic state of stress on the Nojima and Rokko fault systems. Otsuki et al. (1997) have previously interpreted the curved surface striations on the Nojima fault as indicating mixed mode rupture at trans-Rayleigh velocity, based on the dynamic simulations of Andrews (1994). We reanalyse this data because the subsequent work of Guatteri and Spudich (1998) indicates that temporal rake rotation is a sensitive indicator of stress level, and it can have several causes in addition to that proposed by Andrews (1994).

In this paper we apply a method introduced by Spudich (1992) to the Hyogo-ken Nanbu rake data, and we infer that the stress level on the Nojima fault immediately before the earthquake was low (7 - 17 MPa), with coseismic fractional stress drop being complete near the surface and being about 32% at depth in the crust. These numbers are roughly consistent with the post-earthquake stresses measured by Ikeda et al. (1997) in a borehole near the Nojima fault at Hirabayashi and by Tsukahara et al. (1997) in a borehole near Ogura. We will first present the observational data and our method of analysis, and we will defer to the Discussion section a detailed account of the reliability of our observations and inferences.

Striations on the exposed Nojima fault surface

The first source of evidence for temporal rake rotations during the Hyogo-ken Nanbu earthquake comes from striations engraved during the earthquake onto the uplifted fault surface of the Nojima fault. Otsuki et al. (1997) observed striations at 10 locations (Figure 1). At locations 2 and 7 sufficient lengths of striations were available to enable reconstruction of slip paths that we could use to infer initial stresses (Figure 2). At location 2 the bedrock is granite that is highly altered by weathering and hydrothermal processes. Striations I and II (Figure 2) formed first and were strong and smooth, suggesting high slip speed. Striation III formed last, and was weak and rough, suggesting low slip velocity. At location 7 the footwall is composed of fine to medium sandstone, siltstone and mudstone. The fault is associated with an open fissure several cm wide and a zone of fault clay several cm wide. The striations themselves were inscribed in a film of unconsolidated fault clay several mm thick produced by shear during the earthquake. The fragility of the clay in which the striations were inscribed guarantees that the complicated set of striations was formed in the Hyogo-ken Nanbu earthquake.

Although striations at the other locations were inadequate to allow inference of stress levels, these striations, combined with the final offset data of Awata et al. (1996) imply temporal rake rotations at several of the striation locations. Location 1 showed evidence of rotation on a small fault connecting larger fault segments, but this fault carried only a small fraction of the total slip and was not indicative of the main tectonic motion. At location 3 were found multiple striations consistent in orientation with those observed at location 2, indicating similar rake rotations. However, owing to post-earthquake slumping of the fault scarp, only short lengths (~10 cm) of striations were preserved. These striations were too short to construct a reliable slip path. Comparison of a long striation and surface offset at location 6 indicated that NW-side up motion preceded the main right-lateral motion, but no reliable slip path could be constructed because the rake of the dip-slip motion was not recorded. Location 5 striations had rakes that were consistent with the rake inferred from the total surface offset, implying straight slip paths. We did not use striation data from locations 8, 9, and 10 because they were on a branching fault and had small slip, and we did not use striations from location 4 because they were not indicative of the true slip path; faulting here was associated with crack opening and with plastic deformation of deeply weathered bedrock and surface soil.

In this paper we interpret curved striations as indicators of stress levels. Before we discuss other possible causes of curved striations, it will be helpful for the reader to understand the physical basis for our interpretation. Consequently, we defer a more detailed account of the plausibility of our interpretation until the Discussion section.

Dislocation models

Dislocation models derived by Ide et al. (1996), Wald (1996), and Yoshida et al. (1996) from geodetic data and from local and teleseismic ground motion data provide the second source of evidence for rake rotations during the earthquake. Figure 3 shows dislocation at successive time steps in an interpolated version of the Yoshida slip model. (We had to interpolate the slip models onto a denser spatial and temporal grid in order to calculate the associated stress changes, as described later. Figure 1 shows Yoshida's original 4 km 4 km subfaults and our interpolated 1.5 km 1.5 km subfaults.) A strong temporal rake rotation is shown at subfault 145 and its neighbors in Figure 3a, with the NW side of the fault initially moving at high velocity right-laterally and downward, and subsequently moving more slowly right laterally and upward. Numerous subfaults show temporal rake rotations, but because the rotations of subfault 145 and its neighbors occur at a local maximum of slip on the Nojima fault, they are the most likely rotations to be resolved by the ground motion data, and we will concentrate on determining an initial stress for these subfaults. Little temporal rake rotation is seen on the Rokko fault system (Figure 3b), except around subfault 201, which has relatively little slip and might not be well resolved by the strong motion data. We note that the Hyogo-ken Nanbu earthquake shows strong spatial rake rotations, with the final slip direction varying by as much as 80^o from the Nojima to the Rokko faults, similar to the 1989 Loma Prieta earthquake (Beroza, 1991; Wald et al, 1991; Guatteri and Cocco, 1996).

In the remainder of this paper we concentrate on the Yoshida slip model in preference to those by Ide and Wald for three reasons. First, of the three models, the surface slips in Yoshida's model are the most consistent with the observed slips and striation curvatures observed on Awaji Island. Second, as will be demonstrated later, Yoshida's model is the most consistent with certain dynamic relations between slip velocity and stress change that any valid dislocation model should satisfy. Third, Yoshida's and Ide's models agree on the direction of curvature of slip paths for subfaults 145 and 201, for which we derive initial stresses. Wald's model is rather different for these subfaults. Ironically, Yoshida's model has the least temporal rake rotation of the three models. Ide's model has spectacular temporal rake rotations.

Spudich (1992) presented a geometric method for determining initial stress from a dislocation model of an earthquake. In this section we will review his result briefly, and we will reformulate it as a Bayesian calculation of the probability density function (pdf) of initial stress, which is more suitable for handling data with errors.

Consider a planar fault in a cartesian coordinate system, with x₁ directed along strike (roughly NE for the Nojima and Rokko fault systems), x₂ directed down dip, and x₃ being directed roughly NW. Slip or dislocation u is the displacement of the x₃>0 (NW) side of the fault minus the displacement of the x₃<0 side of the fault. The following definitions and development apply to a particular point on the fault in a discretized model of faulting. Let s be the stress tensor in the medium before the earthquake, with components s _ij. We define "initial stress" s to be the two component shear traction vector ( s ₁₃ s ₂₃)^T evaluated on the fault surface, where s ₁₃ > 0 corresponds to a right-lateral shearing force and s ₂₃ > 0 corresponds to a NW-down shearing force. s here corresponds to s ^o in Spudich (1992). (We will use bold characters to denote two-component vectors lying in the fault plane, and our notation generally follows that of Spudich, 1992.) Let t _k be the perturbation in shear traction on the fault plane at time step k caused by motions on the fault, and let v_k be the slip velocity vector at time step k, directed along unit vector. Given a dislocation model of an earthquake, and t _k are uniquely specified by the model and may be obtained from the model in a variety of ways (discussed later). The initial stress vector s can then be derived from and t _k by making the geometric construction shown in Figure 4 (adapted from Spudich (1992) equation (9) and Figure 3, replacing displacement-equivalent tractions with simple tractions). Briefly, for each time step the stress drop vector ( - t _k) is plotted and a semi-infinite line in the direction of slip velocity is extended from the end of the stress drop vector. s is found at the intersection of all the lines.

Physically, the construction of Figure 4 is simply a geometric picture of the statement that slip velocity is parallel to frictional stress. At time step k the total shear stress on the fault, which is identically the frictional stress f _k, is the sum of the initial stress and the stress change, f _k = s + t _k, or s = - t _k + f _k . This latter equation is shown vectorially in Figure 4 for three time steps. Isotropic friction requires that is parallel to f _k, yielding the construction. Rakes rotate with time because part of the contribution to t _k comes from the dynamic stresses radiated from other places on the fault before time step k. Those dynamic stresses are not necessarily parallel to the initial stress, especially when the initial stress direction itself varies over the fault plane (Guatteri and Spudich, 1998), and these deviant dynamic stresses force the slip velocity away from the initial stress direction.

Bayesian approach

For practical application we have chosen to derive a probability density function (pdf) for s using a Bayesian approach derived from the geometric construction of Figure 4. We do this because in a real slip model, and t _k have errors, so that the direction and the end point of the semi-infinite G _k line are uncertain, and consequently the pdf of s is not Gaussian. Let denote the pdf of some quantity y. From Spudich (1992, equation 9) the line G _k is parameterized by (with a slight change of notation) s _k = c_k - t _k, c_k 0. Here s _k is the initial stress consistent with and t _k, and c_k is an undetermined constant. The pdf of s _k is

where the asterisk indicates convolution. Let r and q be cylindrical coordinates in the fault plane. Slip velocity v_k is not known exactly and has a pdf , so the direction q _k of has a pdf

i.e. the pdf of the slip direction is a radial projection of the pdf of the slip velocity vector. This relation gives the intuitively satisfying result that if the error in slip velocity is assumed to be independent of the magnitude of the velocity, then the direction of slip velocity is best determined when slip velocity is large. We used a slight modifiction of the above procedure when we inferred slip velocity from striations. In this case since slip direction was measured directly, we simply assumed was Gaussian with a 10^o standard deviation.

Because we know something about the slip velocity direction but nothing about the magnitude of c_k, we write

where C(r) is an undetermined function of radial distance in the stress plane (i.e. stress magnitude). C(r) is our "Bayesian prior," which we select based on a reasonable assumption. Two possible choices of C(r) are: 1) select it so that the probability of finding s in some annulus having width Ds is independent of , i.e. we want a uniform probability of deriving any stress magnitude from the data, or 2) select C(r) so that we are equally likely to find s in any region of fixed area Ds ₁ x Ds ₂ on the stress plane, which means we are more likely to derive high stress magnitudes from the data because much more of the stress plane has high stress magnitude. Initially choice 1 sounds most reasonable, but this choice leads to a form of C(r) which is rather complicated and which must be a function of the total number of time steps in order to give a sensible in the case of no rake rotation. Consequently, we selected choice 2, and we set for r< s _m' where s _m is an a priori upper bound to stress magnitude (we chose 100 MPa), and C(r) = 0 above s _m. The value of C(r) was selected to give unit probability. Tarantola and Valette (1982, section 3) give a good discussion of this type of choice. We emphasize that since our choice of C(r) is arbitrary, we will not attempt to interpret the radial dependence of in any detail. In particular, when is an open-ended wedge (an example will be given later), we will simply state that the initial stress is undetermined by the data.

Estimation of probability density functions for slip velocity and stress change

Since we have three independent dislocation models of the Hyogo-ken Nanbu earthquake, we determined the probability density functions for slip velocity and stress change by examining the statistics of the differences between the three models. The procedure consisted of three steps.

The first step was development of a method to calculate stress changes from dislocation models. Given a dislocation model for the earthquake, we calculated the dynamic stress changes using a modified version of the boundary integral method (BIM) developed by Quin and Das (1989). Rather than allowing the code to determine slip on the fault from a rupture criterion, we inserted the slip from a dislocation model as the boundary condition and allowed the code to calculate the resulting stress changes. It was necessary to modify the BIM to include approximately the effects of a vertically heterogeneous velocity structure and a free surface reflection. We loosely followed the procedure of Quin (1990). The effect of vertical heterogeneity was included by multiplying the halfspace Green's functions by the ratio of impedences at the source and receiver locations. In addition, because in the BIM the subfault dimension is fixed to be the distance a P wave travels in 2 time steps, the actual vertical extent of each subfault varies with depth proportional to P velocity. To account for the free surface reflection we added the contribution of an image source to the halfspace Green's functions. Our implementation of the approximations performed better than that of Quin (1990), and we checked it as described later.

The second step was calculation of stress changes from the Hyogo-ken Nanbu dislocation models of Ide, Wald, and Yoshida. Each model was specified on rather large subfaults (2.5 km to 5 km in linear extent, depending on author), so we interpolated each model onto a grid having subfaults approximately 1.5 km in length and depth (Figure 1). The slip time function in the interpolated subfaults was kept the same as that in the larger subfault of each author, and the rupture time was lagged between interpolated subfaults in order to simulate rupture propagation within the authors' original subfaults. In addition, we adjusted slip amplitude by bilinear interpolation to reduce slip discontinuities at the edges of the original subfaults. Moment was preserved in this amplitude adjustment. From each interpolated model we calculated dynamic stress changes. Our interpolated model consisted of only a single plane, whereas Yoshida's and Wald's models used two separate planes. Thus, the stresses we calculate in interpolated subfaults at the edges of the original planes are probably not correct, and we refrain from using them in subsequent interpretation. Except for regions at the edges of the original fault planes, our stresses compared well to those calculated for the Wald model by Day et al. (1997) using a finite difference method. Our interpolation method is rather crude because, strictly speaking, an inverse approach should be used to determine a finer grid of slip functions that preserve the moment rate functions of a coarser grid. However, we will demonstrate later that the interpolation does not strongly affect the inferred stresses.

The third step was examination of the differences of the slip velocities and stress changes between the three authors' models. For each interpolated subfault in the Wald model, which was spatially the most extensive, we found the nearest interpolated subfault in the Ide and Yoshida models. For each of these three subfaults we removed the leading zeros from the slip velocity and stress change time series, found the mean slip velocity and mean stress change time series, and then calculated the absolute value of the difference of each author's slip velocity and stress change time series from the mean. The distributions of differences (both slip velocity and stress change, for strike-slip and dip slip motion) were very well fit by exponentials, not Gaussians. Fitting the strike-slip and dip-slip distributions by exponentials of the form exp(ax) and exp(bx), respectively, and multiplying the a and b coefficients by to correct for the fact that the mean velocity and stress change have errors, gives a_v=-3.3 s/m and b_v=-4.6 s/m for slip velocity and a_s=-0.19 MPa^-1 and b_s=-0.20 MPa^-1 for stress change. Note that we do not actually use these numbers unmodified in the calculation of pdfs because they assume all models are equally good. We believe Yoshida's model is the most correct of the three, so we will modify the coefficients accordingly below. We will refer to the above a and b values as the "conservative error coefficients."

The form we use for the pdf of stress change is

,

where . This pdf decays exponentially away from - t _k as a function of argument e. A similar form is used for the slip velocity pdf. This form is actually an exponentially decaying pyramid, and in retrospect it probably would have been better to use a form lacking corners, but this simple analytic form facilitated certain integrations used in the numerics. In addition, and always appear in integrals, so their exact shapes are not critical to the result. Note that because and are decaying exponentials, is also basically a decaying exponential function of angle, and convolution of with preserves the decaying exponential character. This will be relevant when we assess significance below.

Calculation of posterior probability and significance

For each subfault the posterior probability is

where is normalized to have a peak amplitude of unity, and K is the number of time steps used. Strictly speaking, both and are simply density functions, not pdfs, since they are not normalized to unit integral. However, in the interpretation of data we are interested in the value of stress that maximizes , (i.e. the maximum a posteriori probability stress, or "MAP" stress), which enables us to develop a statistic discussed later that assesses the degree of self-consistency of all our assumptions. In practice we use only the K time steps in which slip velocity exceeds max (1/a_v,1/b_v), i.e. we only include pdfs for time steps in which the slip velocity is shifted by more than the 1/e distance from the origin. Equation (5) assumes statistical independence of all time steps, which is not actually true. While the dependence of the time steps probably affects the shape of the pdf , tests we performed indicated that the location of the MAP stress is negligibly biased.

We have developed a P statistic, similar to the c ² statistic, that allows us to assess whether the MAP stress fits the data too well, or too poorly, given our estimates of expected errors (a and b coefficients). For example, suppose for a particular subfault the maximum value of is 10^-43, and for its neighbor the maximum value is 1, and further suppose that for each subfault the product is taken over 10 time steps. In the former case the MAP stress is formed from all the extreme tails of the exponential distributions, possibly indicating that the actual errors in the data far exceed the nominal errors (a and b coefficients) used in our pdfs. In the latter case, the individual pdfs agree at the MAP stress far better than they should, given our nominal error estimates. Because the c ² statistic was developed for Gaussian distributions we cannot use it directly. However, we followed the derivation of Papoulis (1965, p. 251) using our exponential distribution to derive a P statistic. Let Z be a set of n independent random variables z_i,i =1,...,n, having exponential distributions. Their joint density function is g(Z) = exp( -z₁ -z₂...-z_n). Given a scalar p, 0 < p 1, the probability of drawing at random a set Z such that g(Z) p is
P(n, - 1n(p)), where P is the incomplete gamma function (Press et al., 1986, p. 160). Since is the product of decaying exponentials (in angle q , ignoring the convolution with , which basically does not change their decaying exponential nature), we choose p to be the maximum value of , and for n we use K-2, since two degrees of freedom are removed by searching over two stress directions to maximize p. P values near 0 and 1 indicate that the MAP stress fits the data much better or much worse, respectively, than expected given our nominal errors.

We have chosen to concentrate on the Yoshida model in part because it satisfies the dynamic constraint of Figure 4 considerably better than the Wald and Ide models, as quantified by our P statistic. A dynamically valid model should satisfy that constraint in every subfault. For example, we calculated and P for each subfault for all models assuming was a delta function in (1) (i.e. assuming the stress changes could be derived exactly from the models) and using a_v=b_v=20 s/m (this choice will be justified below). The average values of P for the Yoshida, Ide, and Wald models were 0.076, 0.36, and 0.52, respectively.

The MAP stresses calculated below are based upon the assumption that the rake rotations in Yoshida's model are real. If we calculate for all subfaults of the Yoshida model using the conservative error coefficients, derived from the assumption that all three dislocation models are equally good, then no unique MAP stress is obtained for any subfault, i.e. the 95% confidence region on s is enormous and seismologically uninformative (an example will be shown later). This occurs because the error estimates are so large that none of the temporal rake rotations is significant, so the usual 1D result is obtained; the motions are consistent with any s collinear with the overall slip.

For calculation of our preferred MAP stresses from the Yoshida model, we used the error coefficients a_v=b_v=20 s/m and a_s=b_s=0.5 Mpa^-1, which correpond to plausible assumed 95% confidence limits of ±15 cm/s in slip velocity and ±6 MPa in stress change. These error estimates are more optimistic than the conservative error coefficients, consistent with our belief that the Yoshida model is more accurate. These particular values were chosen by trial and error because they were the smallest coefficients (largest error estimates) that gave well defined MAP stresses for the subfaults having large rake rotations.

We caution that the temporal rake rotations in Yoshida's model might not be well resolved. Yoshida's slip model shows two areas having temporal rake rotations, one around subfault 145 on the Nojima fault and around subfault 201 on the Rokko fault system (Figure 3). In both regions Ide's model shows the same direction of temporal rotation, although Wald's model has the opposite direction of rotation. Because the region on the Nojima fault has a lot of slip, the rotations there are more likely to be real than the rotations around subfault 201 on the Rokko fault. However, there are few strong motion stations near the Nojima fault to constrain the dislocation solution there. Consequently, we guess that there is a 70% probability that Yoshida's rake rotations on the Nojima fault are correct, and there is a 50% probability that the Rokko fault rotations are correct.

The difficulties of using the Bayesian approach to estimate stresses are illustrated by three example subfaults. Figure 5 shows geometric constructions for stress for subfaults 12, 145 (a local maximum of slip on the Nojima fault), and 253, which satisfy to different degrees the geometric constraint that the G _k lines should all intersect (Figure 4). See Figure 3 for locations of subfaults. Subfaults 12, 145, and 253 satisfy the geometric constraint very well, decently, and very poorly, respectively. Figure 6 shows the correponding density functions. While Figure 5 indicates that the motions of subfault 12 are beautifully consistent with an initial stress of 10 MPa right-lateral and 3 MPa NW-up, the Bayesian density function is extremely elongated with a poorly resolved peak at 53 MPa right-lateral and with significant density at 100 MPa, our a priori upper limit s _m. The Bayesian formalism is telling us that given our assumed errors, it is extremely unlikely that the actual right lateral stress is 10 MPa. For subfault 145 the Bayesian approach with the assumed errors gives a well-resolved stress of (12, 9) MPa. Note that the pdf has a peak even though the G _k lines do not all intersect (Figure 5). This occurs because the pdf for each time step k, is typically a blurred, truncated wedge, somewhat like the pdf for subfault 253 in Figure 6. The line G _k lies along the axis of symmetry of , which has significant amplitude well away from G_k. For this reason, even when the G _k lines do not all intersect, the corresponding pdfs can have significant overlap. Subfault 253 (and many others) has a wedge shaped density function truncated at s _m. (The wedge would extend to infinity if not artificially truncated.) The literal interpretation of this density is that the stress is probably high rather than low, but recall that the radial dependence of density is controlled by our arbitrary choice of C(r), the "uninformative" prior, and a different choice (such as option 1 above) would have caused a different, lower stress result. Consequently, we will interpret such truncated wedges as meaning simply that we cannot determine the stress level from this density, i.e. that stress is unconstrained on this subfault. In practice, we declared subfaults to have unconstrained stresses if the amplitude of their density functions at the high-stress cutoff equalled or exceeded half of their peak amplitudes. Note that the geometry of Figure 4 makes resolution of large initial stresses very difficult. If the initial stress is much larger that the magnitudes of the stress drops, then the G _k lines are nearly parallel (i.e. little rake rotation), leading to an unbounded pdf like subfault 253 (Figure 5). Consequently, while our procedure can correctly discern low stress subfaults like 145 when there is substantial rake rotation, it cannot distinguish between the low stress/small rake rotation case and the high stress/small rake rotation case.

Initial stresses derived from Yoshida's model are shown in Figure 3. To produce this figure we calculated for all subfaults having more than 15^o rake rotation (219 out of 340). Of the 219, 187 had unconstrained stresses, leaving the 32 stresses plotted in Figure 3. Of course, some of these stresses are probably not reliable because the rake rotations may not be well constrained by the ground motion data. We believe the more reliable stresses are those for subfaults around a local maximum of slip having a large rake rotation, which we have indicated in Figure 3 with diamonds. The pdfs for most of these "diamond" subfaults are similar in extent to that for subfault 145 (Figure 6), which gives a rough idea of the resolution of initial stress.

We demonstrate now that the inferred initial stresses are not highly sensitive to our interpolation method. Figure 5 shows that in the simple geometric construction the initial stress is determined by the amount of temporal rake rotation and by the stress change time series, which determines the left ends of the semi-infinite lines. We will demonstrate two points. First, any sensible interpolation of Yoshida's model would give temporal rake rotations on the Nojima fault. Yoshida's subfaults are 4 x 4 km, and his rupture velocity is 2.5 km/s. In Yoshida's subfaults 33 and 34, which correspond to our subfaults with diamonds (Figure 3a), the slip direction is constant for the first two seconds and rotates to a different rake only in the third second after rupture. The rupture front propagates across a Y subfault in 1.6 s (or diagonally in 2.2 s), so Yoshida's entire subfault is rupturing before the rake of the moment rate function starts to rotate. Thus, temporal rake rotation must occur on the Nojima fault. Rake rotates earlier after rupture initiation on the Rokko fault subfault of interest, so it is possible that spatial variations of rake on the Rokko fault could be aliased into temporal variations in the moment rate function. Second, although the rise times on our interpolated subfaults are too long and the peak velocities too low, the inferred initial stresses are not strongly affected by this. To show this, we created a test interpolated model that consisted of our original interpolated model with rise times cut in half and slip velocities doubled. Although peak dynamic stress changes are amplified somewhat in this test model, static stress change is unaffected, and in general the geometric constructions for initial stress, as in Figure 5, do not change much. Using the test model, we calculated MAP initial stresses at subfaults having diamonds in Figure 3. The average initial stress for these subfaults in the test model was 86% of the average using our original model, and in all cases the test initial stresses fell within the pdf peaks of our original model (e.g. Figure 6). Consequently, we believe that any sensible interpolation of Yoshida's slip model would yield initial stresses similar to ours, at least on the Nojima fault.

In order to determine stresses for surficial subfaults 8 and 11 where we have observed curved striations, we derived slip velocity and stress change time series for these subfaults in the following manner. Here we assume that the observed surface offsets and the slip paths from striations are more accurate than the slip in the surficial subfaults of Yoshida's model. Because the slip paths tell us little about slip speed, we first calculated slip speed time series for subfaults 8 and 11 from Yoshida's slip model, because that model agrees fairly well with the observed striations and surface offsets. Those slip speed time series were mapped onto the directions in the slip paths to determine slip velocity time series. For subfaults 7, 9, 10, 12, 13, and 14, we used a similar procedure with straight slip paths from the geologic offsets of Awata et al. (1996). This gave us slip velocities for subfaults 7-14 generally consistent with the observed surface geologic information and adequate for calculation of stress change on subfaults 8 and 11. We then calculated stress changes on all subfaults using a combined geological-seismological slip velocity model comprised of the above slip velocities in subfaults 7-14 and Yoshida's model in all other subfaults (we refer to this model as the hybrid model henceforth). Calculated surficial stress changes were rather sensitive to our geologic assumptions, but the consequent inferred initial stresses did not depend strongly on the assumptions.

Striations at location 2 (subfault 11) are fairly consistent with an initial stress (right-lateral, NW-down) of (5, 5) Mpa (Figures 7, 8). Striations at location 7 do not satisfy the geometric constraint very well, but they give a MAP stress of (12, 5) MPa. For both locations the MAP stresses are near the stress drop vectors for the last time step, meaning that stress drop is complete. For both locations the conclusion that stress drop is complete is probably more reliable than the specific values derived for the initial stresses. Because the initial slip velocities are strongly dip-slip (Figure 7), their associated pdfs pull the strike-slip component of MAP stress to the lowest value possible, which is related to (although not identical to) the final stress drop. We attempted to check the consistency between stresses derived from the striation slip paths with stresses derived for the same subfaults using Yoshida's dislocations. However, for these locations Yoshida's dislocations give a wedge-shaped pdf for initial stress (like subfault 253 in Figure 6), so the best we can say is that initial stresses derived from Yoshida's dislocations for subfaults 8 and 11 are not inconsistent with initial stresses from striations.

The direction of curvature of striations at locations 2 and 7 agrees with that in the theoretical dynamic models of Guatteri and Spudich (1998). They found that the largest temporal rake rotations occurred at places where the initial stress direction deviated most strongly from the average initial stress direction for the whole fault. At these places the initial slip direction agreed with the initial stress direction, but later slip was in the direction of overall slip for the fault. Initial slip at locations 2 and 7 was dip-slip, with later slip in the strike-slip direction, consistent with the dynamic models.

Comparison with borehole observations

The initial stresses we have inferred from striations and Yoshida's slip model agree surprisingly well with independent evidence about the state of stress obtained from three boreholes drilled after the Hyogo-ken Nanbu earthquake. These were the Nojima-Hirabayshi NIED (Nat'l Inst. for Earth Science and Disaster Prevention) borehole of Ikeda et al. (1996, 1997), the Nojima-Ogura borehole of Tsukahara et al. (1997), and the Ikuha borehole of Ito et al. (1997) of the Geological Survey of Japan (Figure 1). For brevity, we will refer to these three holes as the Hirabayashi, Ogura, and Ikuha holes, respectively.

Ikeda et al. (1996, 1997) made hydrofracture measurements in the 1.8 km deep Hirabayashi borehole located next to the Nojima fault very near to our striation location 2. They observed the maximum (SH) and minimum (Sh) horizontal stresses to increase until the drill penetrated a cataclastic zone at 1.04-1.14 km depth thought to be the Nojima fault, below which SH, Sh, and SV, the vertical stress, became approximately equal (Figure 9). The general orientation of SH was NW-SE, normal to the strike of the Nojima fault. Tsukahara et al. (1997) obtained two measurements of SH and Sh at depths of about 1.5 km (Figure 9) in the Ogura borehole (near our striation location 7), and they also determined from borehole breakouts that the azimuth of SH was N40±5W (H. Tsukahara, personal communication, 1997). At Ikuha, farther southwest along the fault (not on Figure 1), using breakouts Ito et al. (1997) observed SH oriented roughly N20W in the 0.4-0.6 km depth interval, and oriented roughly NW-SE in the 0.8-1.0 km depth interval.

To compare these borehole observations with our pre-earthquake stresses, we need the borehole stresses resolved into shear stress on an assumed vertical Nojima fault. The upper bound on this resolved shear stress is (SH - Sh) / 2, which is everywhere 12 MPa or less (Figure 9), and is 4 - 8 MPa below the cataclastic zone in the Hirabayashi borehole. For both the Hirabayashi and Ogura holes, the actual resolved postseismic shear stress on the Nojima fault is essentially zero (indicated by dashed lines in Figure 9) because of the fault-normal compression.

Our analysis of striations implies complete stress drop near the surface, which is consistent with the results from the Hirabayashi and Ogura boreholes. This can be seen in Figure 10, where we compare our initial stress inferred from location 2 and 7 striations, the stress change inferred from our hybrid slip model, and the borehole information on postseismic shear stresses. In plotting the borehole stresses we have assumed that the vertical component of shear stress on the Nojima fault is zero, which is equivalent to assuming that one of the principal stresses is vertical, which is generally assumed in borehole studies.

Combination of all stress measurements (Figure 11) shows some increase of shear stress with depth between 0 and 3 km. Between 3 km and 11 km there is no clear change of initial stress or postseismic stress with depth, although it must be remembered that we have not been able to infer a stress for the vast majority of subfaults, so our sampling is very sparse. In addition, most results below 6.5 km come from the more uncertain Rokko fault. Stress drop at the surface (at least for our two striation locations) is complete, but between 3 and 11 km fractional stress drop is about 32%. This result differs from the low fractional stress drop inferred for the 1992 Landers earthquake by Cotton and Campillo (1995). They interpreted the spatial uniformity of rake in the Landers event, which had a highly heterogeneous slip distribution, to imply a high stress and low stress drop. However, dynamic rupture models by Guatteri and Spudich (1998) show that spatial heterogeneity of strength (and consequent slip) by itself will not cause large spatial or temporal rake rotations if the initial stress direction is uniform, regardless of stress level. Consequently, the spatial uniformity of rake in Landers implies neither high stress nor low stress.

Plausibility of curved striations as stress indicators

We consider crosscutting straight striations, when formed in a single earthquake, and curved striations to be geologic evidence of coseismic temporal rake rotations. Of course, there may be several causes of coseismic temporal rake rotations in addition to the effect of low stress. Near surface soils can be weak and loose, experiencing nonlinear, spatially heterogeneous strains during an earthquake. Fault zones can consist of lenses of rock that jostle past each other as glide planes interact (e.g. Etchecopar, 1984, p. 30). Notwithstanding these possibilities, we believe that the average sense of curvature in an area may be related to the ratio of dynamic and initial stresses. Of course, curved striations formed in the top 2 m of the earth probably do not reflect stresses in the top 2 m; it is more likely that the near surface materials, particularly when weak, are riding passively on the deeper more competent materials. We estimate that the average surficial striations on the Nojima Fault probably reflect the motions at depths of 200-300 m, based on the assumption that the vertical distance over which rake rotations are correlated equals the horizontal distance over which they are correlated. Striation curvatures at locations 2 and 3, separated by about 200 m, were correlated. At larger separations, however, rake rotations were uncorrelated. Rotations were not correlated at locations 6 and 7, about 600 m apart. Slip at location 5 was pure strike-slip and was uncorrelated with the rotation at location 6, 1.1 km distant. Since location 7 is about 800 m from the Ogura borehole, the lack of rotation correlation over distances of 600 m raises the question of whether location 7 initial stresses can be compared to the Ogura borehole results. The comparison of locations 6 and 7 shows that the direction of initial stress can vary over 600 m. However, the fractional stress drop and the initial stress magnitude may have correlation distances longer than 600 m. Location 7 striations and the Ogura borehole are consistent with low initial stress and complete stress drop. Because the rake rotated at location 6, it is plausible that this location also had low initial stress and complete stress drop. Its direction of rotation was consistent with the dynamic models of Guatteri and Spudich (1998), having initial dip-slip motion followed by strike-slip motion. Consequently, none of the data from locations 6, 7, and the Ogura borehole are inconsistent with the hypothesis of low initial stress and complete stress drop. It would be very interesting to look at the slip surface of the Nojima fault as it appears in recovered cores from the Hirabayashi and Ogura boreholes to see whether the fault surface at depth shows curved or crosscutting striations.

Curved and crosscutting striations have been seen in other earthquakes, suggesting that there may be some repeated underlying tectonic contribution to their formation. After the 1992 Landers earthquake we observed curved and crosscutting straight striations on vertical scarps about 1.5 m high near Bessemer Mine Road (at 116^o 33.0' W, 34^o 32.2' N), where horizontal offsets reached 6 m. The Landers observation is significant because the exposed material was a weathered granite, not a loose soil, and the crosscutting striations were fresh and had obviously formed in the same event. Other examples of historic coseismic temporal rake rotations are the 1957 Gobi-Altai earthquake (Florensov and Solonenko, 1965, pp. 345-348), the 1969 Pariahuanca, Peru, earthquake (Philip and Megard, 1977, fig. 8), the 1974 Izu-Hanto-oki, Japan, earthquake (Kakimi et al., 1977), the 1980 El Asnam, Algeria, earthquake (Philip and Meghraoui, 1983; Philip, personal communication, 1997), and the 1995 Neftegorsk, Russia, earthquake, for which Shimamoto et al. (1996) reported striation rakes that disagree with the overall surface offset, implying temporal rake rotations. Etchecopar (1984) described a paleoearthquake having "brutal" curved striations, and another paleoearthquake left curved striations on the Polienas fold on the northern flank of the Vercors Massif, near Grenoble, France (J.P. Gratier, personal communication, 1997). Of course, such paleoearthquakes are especially interesting because the striations may have been formed at seismogenic depth away from near surface complications.

Implications of low shear stress

Our results for the Nojima fault support the hypothesis that major crustal faults are slipping at low shear stress levels. Between 3 and 11 km depth we see no consistent depth dependence of initial stress, and we see a very slight tendency of postseismic stress to decrease with depth (Figure 11). This lack of a strong depth dependence is in accordance with the hypotheses of Byerlee (1990, 1993) and Rice (1992) that elevated fluid pressures in a highly permeable fault zone sandwiched between impermeable walls would lead to depth-independent fault strength. Zhao et al. (1996) have found evidence for elevated fluid pressure at the hypocenter of the Hyogo-ken Nanbu earthquake. Our observations of low shear stresses agree with those of Tsukahara et al. (1996), who observed very small values of (SH-Sh) in fracture zones penetrated by a 2 km deep borehole near Ashio, 100 km north of Tokyo. However, contrary to the hypotheses of Byerlee and Rice, they did not find highly elevated water pressure in the fracture zones; rather, they appeal to the effect of clay minerals to decrease frictional strength in fault zones.

Our work and that of Bouchon and Streiff (1997) and Guatteri and Spudich (1998) emphasize the important interrelationships between the pre-seismic stress tensor in a region, the shape of the fault surface, and the resulting dynamics of rupture. If faults slip at low shear stress, then one of the principal stresses will be nearly fault normal (in the case of the San Andreas and the Nojima faults, that stress is S₁, the maximum principle stress). If the fault surface has a heterogeneous orientation in such a stress field, then the direction of the initial stress can vary strongly over the fault surface (Árnadóttir and Segall, 1994), causing the spatial variations of rake, such as those seen in the Hyogoken-Nanbu and Loma Prieta earthquakes (Beroza, 1991; Árnadóttir and Segall, 1994; Guatteri and Cocco, 1996), and consequent temporal rake rotations. The strong spatial variation of rake in the Hyogoken-Nanbu earthquake indicates that the initial stress direction was highly heterogeneous on the fault, possibly caused by variations in the fault geometry. The results of Guatteri and Spudich (1998) imply that the heterogeneous initial stress direction caused the temporal rake rotations we discuss here. This mechanism is completely different from that discussed by Andrews (1994).

It is not clear whether the existence of temporal rake rotations necessarily implies that fault friction is isotropic. Faults can have grooves of varying scales (distinct from striations) that define a preferred slip direction (e.g. Scholz, 1990, pp. 146-147). Such a preferred direction would imply an anisotropic friction law, in which slip velocity is not collinear with frictional traction acting across an average fault surface, and would impede temporal rake rotations. (We note that anisotropic friction is a macroscopic artifice to simulate the effects of isotropic friction acting on groove surfaces that are not parallel to the assumed average fault surface. In the geometric construction of Figure 4, anisotropic friction would cause the slip velocity vectors to rotate away from the frictional stress vectors.) Sleep (1997) has pointed out that macroscopic grooves will not develop on a fault surface that slides with various rakes during earthquakes, although he shows that even under this condition some anisotropic friction might be expected. Since fault jogs (Scholz, 1990, p. 148) are not necessarily insuperable barriers to slip, grooves may not be, either. Thus it becomes critical to make field geologic observations of fault geometry to understand how heterogeneity in the fault orientation at all scales affects friction over slip distances of a few meters.

Low strength on faults may help explain why ground motions from compressional environments are generally higher than those in strike-slip environments (e.g. Boore et al., 1997) while extensional regime earthquakes have lower motions (Spudich et al., 1997). McGarr (1984) suggested this correlation of ground motion with state of stress, and he postulated that earthquake stress drops (generally 3-10 MPa) were correlated with fault strength, which in his high stress model was about 250 MPa and 66 MPa at 10 km depth in compressional and extensional regimes, respectively. In a high stress environment it is not obvious why such relatively small stress drops should correlate with fault strengths 10 times larger. However, in an environment in which pore pressures reduce the fault strength to 10-20 MPa, it is clear that stress drop should correlate with any systematic difference in shear stress level in different environments.

The occurrence of temporal and spatial rake rotations can affect the relation between radiated energy and seismic moment. In the case of the Hyogo-ken Nanbu earthquake, at very long periods the NW-down motion of the Nojima fault and the NW-up motion of the Rokko fault (Figure 3) sum to nearly zero net vertical motions and contribute almost nothing to the long period moment. However, at shorter periods viewed from many azimuths (for example, in the city of Kobe) the ground motions from these two regions arrive at distinctly different times and do not destructively interfere, leading to the possibility that the radiated energy is unusually high for the moment of this event. A similar situation can be caused by temporal rake rotations such as that of subfault 145. The basic principle is that if during an earthquake a fault slips 1 m right lateral and immediately afterward slips 1 m left lateral so as to have no final offset, the earthquake has zero moment but radiates energy. The temporal rotation of subfault 145 occurs because the dip-slip component of slip velocity reverses during slip. The rotation appears in the dislocation solution because it causes additional (unidentified) energy to be observed at some seismic stations, compared to the energy that would have been radiated if the slip path had been straight. However, the seismic moment of subfault 145 depends only on the final offset, not on the sinuous slip path leading to that final offset.

Low stress and consequent temporal rake rotations may help explain some unusual polarizations observed in strong ground motion recordings. During the 1985 M 8.1 Michoacan earthquake Anderson et al. (1986) observed that the Caleta de Campos station initially moved NW, and subsequently moved southward. This station was directly above the rupture surface, and the periods were so long that the station displacement may have indicated the motion of the hanging wall of the fault. If so, these motions indicate temporal rake rotations, and the sense of the rotation is consistent with the theoretical results of Guatteri and Spudich (1998), since the later motion is consistent with the overall mechanism of the earthquake. Interestingly, Anderson et al. (1986) use the event's low stress drop and its lack of a heat flow anomaly to argue that the initial stress for this event was probably less than 10 MPa. Another candidate for temporal rake rotation is the 1994 Northridge earthquake. Spudich et al. (1996) noted that the largest acceleration pulses at several strong motions in the south, including Tarzana and Santa Monica, occured on the east component of motion, contradictory to the polarization expected from the overall mechanism of the event.

Finally, the issues raised in this paper suggest two possible ways to improve the practice of ground motion inversion. If future investigators believe that initial stresses are high or that friction is highly anisotropic, then they can remove from their inversions the free parameters that enable temporal rake rotations. Otherwise, the inversions may be improved by inclusion of initial stress constraints, like that of Figure 4, that reduce the number of free parameters in the solution. While the constraint of Figure 4 would be cumbersome to implement in present-day single-step inversion algorithms, useful information could come from iterative inversion algorithms that alternate kinematic inversions and dynamic simulations, like that of Fukuyama and Mikumo (1993). Slip rakes determined from an initial kinematic inversion could be used as the initial stress rakes. Subsequent dynamic simulations could be done for a set of differing initial stress magnitudes. Some value of initial stress will fit the data best.