- Split View
-
Views
-
Cite
Cite
Y. Kaneko, P. M. Shearer, Seismic source spectra and estimated stress drop derived from cohesive-zone models of circular subshear rupture, Geophysical Journal International, Volume 197, Issue 2, May 2014, Pages 1002–1015, https://doi.org/10.1093/gji/ggu030
- Share Icon Share
Abstract
Earthquake stress drops are often estimated from far-field body wave spectra using measurements of seismic moment, corner frequency and a specific theoretical model of rupture behaviour. The most widely used model is from Madariaga in 1976, who performed finite-difference calculations for a singular crack radially expanding at a constant speed and showed that |$\skew4\bar{f}_{\rm c} = k \beta /a$|, where |$\skew4\bar{f}_{\rm c}$| is spherically averaged corner frequency, β is the shear wave speed, a is the radius of the circular source and k = 0.32 and 0.21 for P and S waves, respectively, assuming the rupture speed Vr = 0.9β. Since stress in the Madariaga model is singular at the rupture front, the finite mesh size and smoothing procedures may have affected the resulting corner frequencies. Here, we investigate the behaviour of source spectra derived from dynamic models of a radially expanding rupture on a circular fault with a cohesive zone that prevents a stress singularity at the rupture front. We find that in the small-scale yielding limit where the cohesive-zone size becomes much smaller than the source dimension, P- and S-wave corner frequencies of far-field body wave spectra are systematically larger than those predicted by the Madariaga model. In particular, the model with rupture speed Vr = 0.9β shows that k = 0.38 for P waves and k = 0.26 for S waves, which are 19 and 24 per cent larger, respectively, than those of Madariaga. Thus for these ruptures, the application of the Madariaga model overestimates stress drops by a factor of 1.7. In addition, the large dependence of corner frequency on take-off angle relative to the source suggests that measurements from a small number of seismic stations are unlikely to produce unbiased estimates of spherically averaged corner frequency.
1 INTRODUCTION
Characterization of earthquake source parameters is important for understanding the physics of source processes and seismic hazard. One of the key earthquake source parameters is the static stress drop Δσ, that is, the difference between the average state of stress on the fault before and after an earthquake. Stress drop provides hints on the scaling of earthquake-source parameters (e.g. Aki 1967; Abercrombie 1995; Allmann & Shearer 2009) and insights into tectonic environments in which earthquakes occur (e.g. Kanamori & Anderson 1975; Allmann & Shearer 2007, 2009). For example, Allmann & Shearer (2009) analysed seismic spectra of 2000 mb ≥ 5.5 earthquakes worldwide and found that their stress drop estimates range from 0.3 to 50 MPa with a median value of 4 MPa. Although there are considerable variations in the estimated stress drops within earthquakes of the same magnitude, the median stress drops appear to be relatively constant with magnitude (Allmann & Shearer 2009). Stress drop is also used as a primary input parameter for ground motion simulations with stochastic techniques for quantification of seismic hazards (e.g. Graves et al.2010).
Stress drop is proportional to the ratio of total slip to rupture size. Because the physical size of the source is typically inaccessible to direct observation, it is commonly estimated from far-field body wave spectra using measurements of corner frequencies fc and a theoretical model of rupture behaviour. The source dimension and the seismic moment, obtained from the low-frequency part of the spectra, are then used to compute stress drop. Thus, although stress drop is fundamentally a static parameter (it can be obtained from geodetic observations following large earthquakes), it is often estimated in a way that relies on the validity of a specific theoretical model of rupture dynamics.
Accurate determination of stress drops from body wave observations is difficult not only because of uncertainties in observations and corrections for path effects, but also due to the fact that the methods for estimating stress drop are derived from a number of assumptions about the dynamics of the source. In particular, there is no agreement among investigators on which theoretical model should be used for estimating the source dimension, hampering comparisons among different stress-drop studies.
A key issue is that the choice of k in eq. (2) depends on which theoretical relationship is being used to relate corner frequency and source radius. Many investigators (e.g. Hanks & Thatcher 1972; Archuleta et al.1982; Ide et al.2003; Oth et al.2010; Baltay et al.2011; Cotton et al.2013) used the model of Brune (1970), which assumed a simple but somewhat ad hoc kinematic model for a circular fault and obtained k = 0.37 for S waves. Others (e.g. Prejean & Ellsworth 2002; Stork & Ito 2004; Imanishi & Ellsworth 2006) used the analytical model of Sato & Hirasawa (1973) in which the rupture nucleates at a point, spreads radially with a constant rupture speed, and stops abruptly at the source radius. In the model of Sato & Hirasawa (1973), which was constructed by assuming that the static solution of Eshelby (1957) holds at every successive instant of rupture formation under uniform stress, k depends on the rupture speed Vr with k p = 0.42 for P waves and k s = 0.29 for S waves for Vr = 0.9β. An unphysical feature of the model of Sato & Hirasawa (1973) is that particle motion ceases at the same instant, everywhere over the fault plane.
The model of Madariaga (1976) has perhaps been the most widely used (e.g. Abercrombie 1995; Prieto et al.2004; Abercrombie & Rice 2005; Shearer et al.2006; Allmann & Shearer 2007, 2009; Yang et al.2009). Madariaga (1976) performed dynamic calculations for a singular crack spreading radially at a constant rupture speed using a staggered-grid finite-difference method and found that k p = 0.32 for P waves and k s = 0.21 for S waves for Vr = 0.9β. Since the rupture front in the Madariaga model is characterized by a singular crack, the finite mesh size and smoothing procedures introduce some artificial effects, which influence the corner frequencies of the far-field spectra and hence the value of k.
As one can see from eq. (3), the difference between Brune's k = 0.37 and Madariaga's k = 0.21 leads to a factor of 5.5 difference in stress drop; thus assuming the same S-wave speed β, the Madariaga model translates to stress-drop estimates 5.5 times smaller. This means that even in an idealized scenario of circular earthquake rupture, the absolute level of stress drop is highly uncertain, motivating further investigation of this problem. (For a comprehensive summary of various k from different theoretical models, see Dong & Papageorgiou 2002.)
Here, we construct a dynamic model of expanding rupture on a circular fault, similar to that of Madariaga (1976), but with a cohesive zone that prevents a stress singularity at the rupture front. Our goal is to study the simplest class of rupture models that are physically realizable, that is, have no stress singularities, and that are generated by simulations with a proper numerical resolution. We simulate scenarios of dynamic rupture nucleating at the centre of the fault, propagating at a constant rupture speed and stopping at the edge (Section 2). Unlike the model of Sato & Hirasawa (1973), this model results in spontaneous healing of slip due to the arrivals of stopping phases over the source area. We then compute the corresponding far-field body wave displacement for a homogeneous elastic whole space using the representation theorem of Aki & Richards (2002). We analyse the behaviour of body wave spectra (Section 3), the relation between the corner frequency and the source radius and the dependence of corner frequencies on rupture speeds (Section 4). We discuss uncertainties in estimation of stress drop and the implications of our results for stress drops reported in previous observational work (Section 5).
2 COHESIVE-ZONE MODEL OF A CIRCULAR FAULT
In this model, the rupture nucleates at the centre of the fault, and the circular rupture front expands at a constant speed close to Vr. Fault growth stops instantaneously as the rupture runs into the zone of zero shear stress outside the source radius a. The level of shear stress outside of the circular source can be arbitrary as long as it is much smaller than the dynamic strength τd. In the limit where the weakening rate Aw becomes infinitely large (Fig. 1b), the model approaches a singular crack model. Note that eq. (4) is different from the commonly used slip-weakening friction law, where the resulting rupture becomes spontaneous, leading to a more complex source model with set of complicated waveforms. For simplicity, we keep Vr fixed, but the model is still dynamic in the sense that we solve for the fault motion given the rupture speed Vr and prescribed dynamic stress drop Δσd (sometimes called effective stress).
We solve the elastodynamic equation coupled with the fracture criterion (eq. 4) using a spectral element method (Komatitsch & Vilotte 1998; Kaneko et al.2008; Kaneko & Lapusta 2010). We choose a large enough computational domain such that waves reflected by the boundaries of the finite domain do not propagate back to the source area. Following the approach of Madariaga (1976), only the y-component of slip and shear traction is solved numerically (Fig. 1a); there is no rake rotation as the circular rupture front advances. To produce well-resolved numerical results, we ensure that there are at least 7–10 computational node points within a cohesive zone. For example, for a model with the smallest cohesive-zone size that we considered, there are 305 node points (or 76 spectral elements) along the radius of the circular fault. The spectral element model incorporates artificial viscosity of the Kelvin–Voight form to suppress spurious oscillations generated by the fault slip at frequencies that are too high to be resolved by the mesh (Appendix A). The dynamic rupture code has been verified through the Southern California Earthquake Center Dynamic Rupture Code Verification Exercise (Harris et al.2009) and has been used in studies of spontaneous dynamic rupture (Kaneko et al.2011; Kaneko & Fialko 2011; Konca et al.2013).
The solution of the 3-D problem formulated is similar to that of the simplified 2-D problem in Madariaga (1976) where a circular symmetry of the source is assumed. Fig. 2 shows a source model with the rate of frictional weakening |$A_{\rm w}^{\prime } = 42$| and rupture speed Vr = 0.9β. The final slip is greatest at the centre of the fault and monotonically decreases towards the edge of the source (Figs 2a and c). The final shear stress within the source is non-uniform (Figs 2b and d) and is different from the dynamic stress τd due to stress overshoot; the same feature was reported in Madariaga (1976). At the centre of the fault, the stress drop is about 2.5 times larger than the dynamic stress drop Δσd although its contribution to the spatially averaged stress drop over the source is relatively small.
Despite the similarities mentioned above, there are notable differences. Due to the finite cohesive-zone size, the resulting fracture energy is non-zero and increases with distance from the hypocentre over the source area (Fig. 2e). The stopping behaviour of the rupture is more complex than in the Madariaga model and is different along the Modes-II and III directions (Figs 2c–h). The healing waves propagate inwards from the edge of the fault at the P-wave speed α in the Mode-II direction and at ≈0.9α along the Mode-III direction, respectively. Yet the eventual slip arrest occurs sooner along the Mode-III direction (Figs 2f–h). This complex arrest leads to slightly different final slip and stress-drop distributions along the Modes-II and III directions (Figs 2a and b), which breaks the circular symmetry in the solution of the problem.
We find a similar degree of stress overshoot in models with a different cohesive-zone size and without the cohesive zone (or instantaneous strength drop) discussed in Section 4. Hence the large stress overshoot is not caused by the difference in fracture energy along the Modes-II and III directions. Indeed, the healing behaviour is qualitatively similar to that of modelled spontaneous rupture on a circular fault (Madariaga et al.1998) where the fracture energy is assumed to be the same in the Modes-II and III directions.
3 FAR-FIELD BODY WAVE RADIATIONS AND SPECTRA
While far-field synthetic seismograms can be computed within the same numerical simulation used to generate the source model, that would require significantly more computational time and memory as the domain needs to be much larger than the source dimension. Since we are only interested in the far-field response of a source model for a simple Earth structure, we use the representation theorem of Aki & Richards (2002) instead.
Fig. 3 shows the magnitude of the far-field body wave displacement and their spectra at different take-off angles θ. Both corner frequencies and spectral fall-off rate depend on θ (Fig. 3). The long-period or low-frequency region is controlled by the seismic moment, whereas the intermediate frequency region is controlled by the width of the pulse and by the size of the fault. The corner frequencies in Fig. 3 are larger or smaller depending on the width of the displacement pulse. For θ ≳ 30°, the corner frequencies are generally smaller due to the directivity effect; the difference between wave arrivals from the near-side and far-side of the fault is larger, resulting in a longer pulse duration and smaller corner frequency (Fig. 3). These features are qualitatively consistent with those of Madariaga (1976). Note that some of the spectra shown in Fig. 3 are not well fit by the spectral function (14), and we will discuss applications of a more general spectral function introduced by Boatwright (1980) in Section 5.
The corner frequencies of displacement spectra at take-off angles sampled every 5° over the focal sphere are shown in Fig. 4(a). Generally, the corner frequency of the P wave is larger than that of the S wave because of the shorter duration of the P displacement pulses. For θ ≳ 30°, the variation is gradual and both the P and S corner frequencies decrease with θ. Although the pattern of the corner-frequency variation is similar to that in the Madariaga model (Fig. 4b), the amplitude of the variation is much larger in the model considered in this study (Fig. 4a). In addition, both P and S corner frequencies for a given take-off angle show some scatter and weakly depend on the angle ϕ shown in Fig. 1(a). This is due to the lack of perfect circular symmetry in the source model described in the previous section.
Since the corner frequency is a function of the take-off angle, we compute the average |$\skew4\bar{f}_{\rm c}$| over the focal sphere. For the model shown in Fig. 4(a), the spherical averages of the P and S corner frequencies are 0.38β/a and 0.26β/a, respectively, which are larger than those in Madariaga (1976). The ratio of the P-to-S corner frequencies is 1.5, consistent with that of the Madariaga model. It is generally agreed that P corner frequency is larger than the S corner frequency for many observed earthquakes (e.g. Molnar et al.1973; Abercrombie 1995).
There is some potential bias in fitting function (14) to the computed spectra. Fig. 5(a) shows spectral fall-off rates as a function of the take-off angle θ. The spectral fall-off rate also depends on θ and ranges from 1.5 to 3.3. Both the corner frequencies and spectral fall-off rates show a similar dependence on θ (Figs 4a and 5a). Generally, a larger fall-off rate corresponds to a larger corner frequency and vice versa, and hence there is a trade-off in these parameters in the fitting procedure. To make sure that this bias is not substantial in the estimation of the corner frequencies, we fix the fall-off rate to be the average value shown in Fig. 5(a) and repeat the least-squares fit. In this case, the variation of the corner frequencies over the focal sphere becomes smaller than that obtained using variable fall-off rates, but the spherically averaged values of corner frequency differ by less than 4 per cent for both P and S waves (compare Figs 4a and 5b). Hence, we conclude that the average values of corner frequencies are well determined and the uncertainty in the estimation of the corner-frequency averages is relatively small. However, because observations are typically available from only a limited number of seismic stations, more accurate spherical averages will generally be obtained when a fixed fall-off rate is used in fitting the spectra, because the scatter in the individual fc values is reduced.
4 RELATION BETWEEN CORNER FREQUENCIES AND SOURCE RADIUS
4.1 Cohesive-zone model with rupture speed Vr = 0.9β
The source model presented so far assumes the rate of frictional weakening |$A_{\rm w}^{\prime } = 42$|. Increasing |$A_{\rm w}^{\prime }$| generally results in smaller fracture energy G′, but G′ eventually becomes independent of |$A_{\rm w}^{\prime }$| for models with a cohesive-zone size much smaller than the source dimension. We call this case the ‘small-scale yielding limit’. In this limit, the spherical average of corner frequencies does not depend on the rate of frictional weakening (Fig. 6a), and the source-time histories of the models with |$A_{\rm w}^{\prime } = 84$| and |$A_{\rm w}^{\prime } = 168$| become identical (Fig. 7a).
The apparent increase in the corner frequencies with the fracture energy shown in Fig. 6(a) comes from the way the rupture front advances in the model, where the rupture speed Vr is defined at the intersection of the weakening rate Aw and τ = τd and not at the actual rupture front (Fig. 1b). Because of this formulation, the actual rupture front is ahead of the location Vrt (Fig. 1b) and encounters the edge of the circular fault sooner than the time a/Vr. As a result, the eventual source duration becomes shorter in the model with a smaller frictional-weakening rate |$A_{\rm w}^{\prime }$|, or larger fracture energy G′ (Fig. 6a). This is why the models with larger fracture energy lead to larger corner frequencies (Fig. 6a).
To understand the difference of corner frequencies in the small-scale yielding limit and Madariaga's model, we also consider a singular crack model where |$A_{\rm w}^{\prime }$| is infinite and hence the cohesive-zone size is zero. The stress change at a point on the fault for this zero-cohesive-zone model is shown in Fig. 6(b). Obviously, the shear stress at the rupture front is finite and is not well resolved numerically. In this case, we find that the spectra are contaminated by numerical noise at high frequencies that hampers the fitting of function (14). To remedy this, we use 10-times-larger numerical damping η = Δt (Appendix A) in the source model to reduce the high-frequency noise. We then follow the same procedure described above to compute the spherical averages of corner frequencies for this case.
Interestingly, we find that the spherical averages of P and S corner frequencies in the zero-cohesive-zone case are still larger than the prediction of the Madariaga model and remain about the same as those in the small-scale yielding limit (Fig. 6a). Perhaps the corner frequencies are not much affected by the details of the weakening process, with its length scale much smaller than the source dimension. The slight difference in the corner frequencies between the zero-cohesive-zone case and the small-scale yielding limit is caused by the larger numerical damping, which affects both high and low frequencies of the source-time function (Fig. 7b).
4.2 Cohesive-zone models with different rupture speeds
We consider other source models in the small-scale yielding limit for several subshear rupture speed from 0.5β to 0.9β. The spherical averages of corner frequencies for the model with different rupture speeds are shown in Table 1. The value of k generally increases with the rupture speed because the source duration becomes shorter as the rupture speed increases, consistent with the predictions of other models (Sato & Hirasawa 1973; Madariaga 1976). Compared to the P corner frequencies, the S corner frequencies are less affected by the rupture speed (Table 1). The dependence of the P and S corner frequencies on rupture speeds is similar to that of (Sato & Hirasawa 1973). Note that non-dimensional corner frequencies k do not depend on the prescribed dynamic stress drop Δσd as long as the cohesive-zone size in the source model is much smaller than the source dimension.
Vr/β . | k p . | k s . | |$k^{p}_{\rm Mada}$| . | |$k^{s}_{\rm Mada}$| . | |$k^{s}_{\rm Brune}$| . | |$k^{p}_{\rm S\&H}$| . | |$k^{s}_{\rm S\&H}$| . |
---|---|---|---|---|---|---|---|
Infinite | 0.37 | ||||||
0.9 | 0.38 | 0.26 | 0.32 | 0.21 | 0.42 | 0.29 | |
0.8 | 0.35 | 0.26 | 0.39 | 0.28 | |||
0.7 | 0.32 | 0.26 | 0.36 | 0.27 | |||
0.6 | 0.30 | 0.25 | 0.34 | 0.27 | |||
0.5 | 0.28 | 0.22 | 0.31 | 0.24 |
Vr/β . | k p . | k s . | |$k^{p}_{\rm Mada}$| . | |$k^{s}_{\rm Mada}$| . | |$k^{s}_{\rm Brune}$| . | |$k^{p}_{\rm S\&H}$| . | |$k^{s}_{\rm S\&H}$| . |
---|---|---|---|---|---|---|---|
Infinite | 0.37 | ||||||
0.9 | 0.38 | 0.26 | 0.32 | 0.21 | 0.42 | 0.29 | |
0.8 | 0.35 | 0.26 | 0.39 | 0.28 | |||
0.7 | 0.32 | 0.26 | 0.36 | 0.27 | |||
0.6 | 0.30 | 0.25 | 0.34 | 0.27 | |||
0.5 | 0.28 | 0.22 | 0.31 | 0.24 |
Vr/β . | k p . | k s . | |$k^{p}_{\rm Mada}$| . | |$k^{s}_{\rm Mada}$| . | |$k^{s}_{\rm Brune}$| . | |$k^{p}_{\rm S\&H}$| . | |$k^{s}_{\rm S\&H}$| . |
---|---|---|---|---|---|---|---|
Infinite | 0.37 | ||||||
0.9 | 0.38 | 0.26 | 0.32 | 0.21 | 0.42 | 0.29 | |
0.8 | 0.35 | 0.26 | 0.39 | 0.28 | |||
0.7 | 0.32 | 0.26 | 0.36 | 0.27 | |||
0.6 | 0.30 | 0.25 | 0.34 | 0.27 | |||
0.5 | 0.28 | 0.22 | 0.31 | 0.24 |
Vr/β . | k p . | k s . | |$k^{p}_{\rm Mada}$| . | |$k^{s}_{\rm Mada}$| . | |$k^{s}_{\rm Brune}$| . | |$k^{p}_{\rm S\&H}$| . | |$k^{s}_{\rm S\&H}$| . |
---|---|---|---|---|---|---|---|
Infinite | 0.37 | ||||||
0.9 | 0.38 | 0.26 | 0.32 | 0.21 | 0.42 | 0.29 | |
0.8 | 0.35 | 0.26 | 0.39 | 0.28 | |||
0.7 | 0.32 | 0.26 | 0.36 | 0.27 | |||
0.6 | 0.30 | 0.25 | 0.34 | 0.27 | |||
0.5 | 0.28 | 0.22 | 0.31 | 0.24 |
5 DISCUSSION
We discuss uncertainties in estimation of stress drops and implications of the results for stress drops and other source parameters reported in observational work.
5.1 Comparison of the input and inferred source parameters
. | |$M^{\prime }_{\rm o}$| (Mo/Δσda3) . | |$E^{\prime }_{\rm r}$||$(E_{\rm r}\mu /\Delta \sigma _{\rm d}^2 a^3)$| . | |$E^{\prime }_{\rm r}/M^{\prime }_{\rm o}$| . |
---|---|---|---|
Input source model | 2.87 | 0.83 | 0.29 |
Inferred from P-wave spectra | 3.16 | 0.72 | 0.23 |
Inferred from S-wave spectra | 3.39 | 0.21 | |
Inferred from P-wave displacement | 2.89 | 0.81 | 0.28 |
Inferred from S-wave displacement | 2.88 | 0.28 |
. | |$M^{\prime }_{\rm o}$| (Mo/Δσda3) . | |$E^{\prime }_{\rm r}$||$(E_{\rm r}\mu /\Delta \sigma _{\rm d}^2 a^3)$| . | |$E^{\prime }_{\rm r}/M^{\prime }_{\rm o}$| . |
---|---|---|---|
Input source model | 2.87 | 0.83 | 0.29 |
Inferred from P-wave spectra | 3.16 | 0.72 | 0.23 |
Inferred from S-wave spectra | 3.39 | 0.21 | |
Inferred from P-wave displacement | 2.89 | 0.81 | 0.28 |
Inferred from S-wave displacement | 2.88 | 0.28 |
. | |$M^{\prime }_{\rm o}$| (Mo/Δσda3) . | |$E^{\prime }_{\rm r}$||$(E_{\rm r}\mu /\Delta \sigma _{\rm d}^2 a^3)$| . | |$E^{\prime }_{\rm r}/M^{\prime }_{\rm o}$| . |
---|---|---|---|
Input source model | 2.87 | 0.83 | 0.29 |
Inferred from P-wave spectra | 3.16 | 0.72 | 0.23 |
Inferred from S-wave spectra | 3.39 | 0.21 | |
Inferred from P-wave displacement | 2.89 | 0.81 | 0.28 |
Inferred from S-wave displacement | 2.88 | 0.28 |
. | |$M^{\prime }_{\rm o}$| (Mo/Δσda3) . | |$E^{\prime }_{\rm r}$||$(E_{\rm r}\mu /\Delta \sigma _{\rm d}^2 a^3)$| . | |$E^{\prime }_{\rm r}/M^{\prime }_{\rm o}$| . |
---|---|---|---|
Input source model | 2.87 | 0.83 | 0.29 |
Inferred from P-wave spectra | 3.16 | 0.72 | 0.23 |
Inferred from S-wave spectra | 3.39 | 0.21 | |
Inferred from P-wave displacement | 2.89 | 0.81 | 0.28 |
Inferred from S-wave displacement | 2.88 | 0.28 |
. | |$E^{s}_{\rm r}/E^{p}_{\rm r}$| . |
---|---|
Inferred from stacked spectra | 19.3 |
Inferred from the integral of squared velocity | 21.8 |
Inferred from |$f_{\rm c}^{p}/f_{\rm c}^{s}$| using the theory of Boatwright & Fletcher (1984) | 7.6 |
Observational work of Boatwright & Fletcher (1984) | 27.3 ± 3.3 |
Observational work of Abercrombie (1995) | 14.3 (from 4.4 to 46.3) |
Observational work of Prieto et al. (2004) | 9 ± 1.5 |
Observational work of Kwiatek & Ben-Zion (2013) | 4.8 (from 0.5 to 50) |
. | |$E^{s}_{\rm r}/E^{p}_{\rm r}$| . |
---|---|
Inferred from stacked spectra | 19.3 |
Inferred from the integral of squared velocity | 21.8 |
Inferred from |$f_{\rm c}^{p}/f_{\rm c}^{s}$| using the theory of Boatwright & Fletcher (1984) | 7.6 |
Observational work of Boatwright & Fletcher (1984) | 27.3 ± 3.3 |
Observational work of Abercrombie (1995) | 14.3 (from 4.4 to 46.3) |
Observational work of Prieto et al. (2004) | 9 ± 1.5 |
Observational work of Kwiatek & Ben-Zion (2013) | 4.8 (from 0.5 to 50) |
. | |$E^{s}_{\rm r}/E^{p}_{\rm r}$| . |
---|---|
Inferred from stacked spectra | 19.3 |
Inferred from the integral of squared velocity | 21.8 |
Inferred from |$f_{\rm c}^{p}/f_{\rm c}^{s}$| using the theory of Boatwright & Fletcher (1984) | 7.6 |
Observational work of Boatwright & Fletcher (1984) | 27.3 ± 3.3 |
Observational work of Abercrombie (1995) | 14.3 (from 4.4 to 46.3) |
Observational work of Prieto et al. (2004) | 9 ± 1.5 |
Observational work of Kwiatek & Ben-Zion (2013) | 4.8 (from 0.5 to 50) |
. | |$E^{s}_{\rm r}/E^{p}_{\rm r}$| . |
---|---|
Inferred from stacked spectra | 19.3 |
Inferred from the integral of squared velocity | 21.8 |
Inferred from |$f_{\rm c}^{p}/f_{\rm c}^{s}$| using the theory of Boatwright & Fletcher (1984) | 7.6 |
Observational work of Boatwright & Fletcher (1984) | 27.3 ± 3.3 |
Observational work of Abercrombie (1995) | 14.3 (from 4.4 to 46.3) |
Observational work of Prieto et al. (2004) | 9 ± 1.5 |
Observational work of Kwiatek & Ben-Zion (2013) | 4.8 (from 0.5 to 50) |
Next, we compare the static stress drop of the source model and that inferred from the spectra. The stress drop of the input source model based on eq. (1) is 1.25Δσd (Table 4). Using kp and ks obtained from the small-scale yielding limit and assuming the correct Mo, the inferred stress drops agree with that of the input model for both P and S waves, as expected. If the Madariaga model is used for the estimation, the stress drop is overestimated by a factor of 1.7 for P waves and 1.9 for S waves, respectively (Table 4). Assuming the correct Er, the corresponding radiation efficiencies ηrad = 2μER/(ΔσMo) are well determined from the spectra based on the k values from the small-scale yielding limit, whereas they are underestimated by the same factors based on the Madariaga model (Table 4). This again highlights the importance of the value of k in the proper estimation of stress drop.
. | Δσ/Δσd . | Rad. eff. ηeff . |
---|---|---|
Modelled stress drop based on Δσ = (7/16)(Mo/a3) | 1.25 | 0.46 |
Modelled stress drop based on |$\overline{\Delta \sigma }_{\rm nonuni}$| (eq. 24) | 1.23 | 0.47 |
Inferred from P spectra using kp from the small-scale yielding limit | 1.25 | 0.46 |
Inferred from S spectra using ks from the small-scale yielding limit | 1.25 | 0.46 |
Inferred from P spectra using kp from Madariaga (1976) | 2.10 | 0.28 |
Inferred from S spectra using ks from Madariaga (1976) | 2.38 | 0.24 |
. | Δσ/Δσd . | Rad. eff. ηeff . |
---|---|---|
Modelled stress drop based on Δσ = (7/16)(Mo/a3) | 1.25 | 0.46 |
Modelled stress drop based on |$\overline{\Delta \sigma }_{\rm nonuni}$| (eq. 24) | 1.23 | 0.47 |
Inferred from P spectra using kp from the small-scale yielding limit | 1.25 | 0.46 |
Inferred from S spectra using ks from the small-scale yielding limit | 1.25 | 0.46 |
Inferred from P spectra using kp from Madariaga (1976) | 2.10 | 0.28 |
Inferred from S spectra using ks from Madariaga (1976) | 2.38 | 0.24 |
. | Δσ/Δσd . | Rad. eff. ηeff . |
---|---|---|
Modelled stress drop based on Δσ = (7/16)(Mo/a3) | 1.25 | 0.46 |
Modelled stress drop based on |$\overline{\Delta \sigma }_{\rm nonuni}$| (eq. 24) | 1.23 | 0.47 |
Inferred from P spectra using kp from the small-scale yielding limit | 1.25 | 0.46 |
Inferred from S spectra using ks from the small-scale yielding limit | 1.25 | 0.46 |
Inferred from P spectra using kp from Madariaga (1976) | 2.10 | 0.28 |
Inferred from S spectra using ks from Madariaga (1976) | 2.38 | 0.24 |
. | Δσ/Δσd . | Rad. eff. ηeff . |
---|---|---|
Modelled stress drop based on Δσ = (7/16)(Mo/a3) | 1.25 | 0.46 |
Modelled stress drop based on |$\overline{\Delta \sigma }_{\rm nonuni}$| (eq. 24) | 1.23 | 0.47 |
Inferred from P spectra using kp from the small-scale yielding limit | 1.25 | 0.46 |
Inferred from S spectra using ks from the small-scale yielding limit | 1.25 | 0.46 |
Inferred from P spectra using kp from Madariaga (1976) | 2.10 | 0.28 |
Inferred from S spectra using ks from Madariaga (1976) | 2.38 | 0.24 |
5.2 Implications for stress drop estimates in observational work
The results presented in this study have important implications for estimated source parameters of earthquakes. While relative stress-drop variations found in previous observational work are not affected by the difference in the assumed value of k, the absolute level of a stress drop based on the model with the small-scale yielding limit is systematically smaller by a factor of 1.7 than that based on the Madariaga model. For example, with the revised values of k, stress drop estimates reported in Allmann & Shearer (2009) range from 0.18 to 30 MPa with the median value of 2.4 MPa. This implies that less accumulated tectonic stress is released by those earthquakes than the previous interpretation. In addition, unusually high-stress-drop earthquakes (≳100 MPa) reported in observational studies may not be actually so large and thus more comparable to the tectonic overburden pressure.
Stress drop is also used to infer other source parameters, such as seismological fracture energy. Abercrombie & Rice (2005) found that the proxy for fracture energy (denoted as G′ in their work) increases with the slip S, from 103 J m−2 at S = 1 mm to 106–107 J m−2 at S = 1 m. Since G′ = (Δσ − 2μEr/Mo)(S/2) and depends on Δσ, which was estimated based on the model of Madariaga (1976), G′ would be smaller with Δσ based on the model with the small-scale yielding limit.
5.3 Uncertainty in corner-frequency estimation due to limited station coverage
Stress drop observations typically exhibit large scatter, even in relatively compact regions (e.g. Shearer et al.2006). A key question is how much of this scatter is real (i.e. true differences in earthquake stress drops) and how much may simply reflect observational uncertainties, such as inaccurate corrections for attenuation. An underappreciated aspect of the Madariaga model, confirmed by our own work, is the large dependence of fc on take-off angle relative to the source (Fig. 4a). Thus, measurements from a small number of seismic stations are unlikely to produce unbiased estimates of spherically averaged fc, even if we assume attenuation corrections are perfect. A factor of two difference in |$\skew4\bar{f}_{\rm c}$| will produce a factor of eight difference in stress drop. As an example, synthetic tests of random take-off angles show that the standard error in stress drop is 33 and 120 per cent for single station P- and S-wave estimates, respectively, which is reduced to 19 and 53 per cent for five-station estimates, and 14 and 35 per cent for 10-station estimates (Fig. 8). These should be taken as lower bounds because it is unlikely that the true station distribution will mimic random take-off angles.
Allowing the fall-off rate to vary produces the best fit to the spectra and more accurate fc estimates for individual spectra. However, this approach yields a greater dependence of fc on take-off angle (compare Figs 4a and 5b) and will require more stations to get an unbiased spherical average. As discussed above, the undersampling errors can be reduced if a fixed fall-off rate is used in fitting for fc. Fig. 9 shows that the standard errors in stress drop for a fixed fall-off rate are systematically smaller than those shown in Fig. 8. Hence, although fixing the fall-off rate worsens the fit to individual spectra, it reduces the dependence of fc on take-off angle and the number of stations required to get an unbiased spherical average.
5.4 On fitting the spectral function: weighted versus unweighted fit
When fitting the spectral function (eq. 14) to the computed spectra, we weight the fit inversely with frequency, which makes all parts of the spectrum in a log–log plot contribute equally. Since the low-frequency part of the spectrum is defined by relatively few points, the weighting generally improves the fit to the low-frequency part. When the fit is done without the weight, the estimated corner frequencies are generally smaller than those with the weight, and the corresponding spherical average is smaller with k p = 0.33 and k s = 0.23 for the source model shown in Fig. 2. However, the unweighted fits are visually worse. Hence, we use weighted fit for all the cases considered in this study. Note that many observational studies (e.g. Shearer et al.2006) also use this inverse weighting in fitting spectral functions to data.
5.5 On averaging of corner frequencies from individual spectra versus from stacked spectra
Uncertainty in stress-drop estimation also includes how the average of corner frequencies |$\skew4\bar{f}_{\rm c}$| is calculated. We follow the approach of Madariaga (1976) where |$\skew4\bar{f}_{\rm c}$| is the average of all the fc estimated for individual spectra. However, often observational studies (e.g. Shearer et al.2006) use stacked log spectra in fitting a spectral function (e.g. eq. 14) because individual spectra tend to be irregular in shape and difficult to fit robustly. Fig. 10 shows log spectra, their average (stack) and the corresponding estimated corner frequencies k for the source model shown in Fig. 2. The values of k p and k s differ from those estimated from the spherical average by 2 and 7 per cent for P and S waves, respectively. Hence, corner frequencies estimated by stacked spectra generally agree reasonably well with those of the spherical average.
5.6 Alternate spectral fitting functions and the meaning of fc
To evaluate uncertainties in estimated stress drops for applications of different spectral functions, we fit eq. (23) to the computed spectra shown in Fig. 3. Compared to the spectral function of Brune (1970), applications of eq. (23) generally lead to better fits to the spectra, that is, the root-mean-square errors are systematically smaller (Fig. 11a). The corner frequencies still show a large dependence on take-off angles θ but their spherical averages are smaller, with k p = 0.34 and k p = 0.23 (Fig. 11b). This means that larger γ allows for the fit to the spectra with smaller corner frequencies. Hence even with the synthetic model, the estimation of corner frequencies and the spherical average depend on the assumed spectral function. Note that the spectral function of Brune (1970) has been used more often in observational studies where the spectra generally are rather noisy and difficult to fit with too many parameters (e.g. Prieto et al.2004; Shearer et al.2006). More stable results for fc are often obtained by reducing the number of parameters (e.g. fixing the fall-off rate), as discussed above.
This discussion highlights the fact that fc is not a precisely defined parameter. Unlike moment (Mo) and radiated seismic energy (Er), corner frequency is not directly related to a physical property of the source. fc is generally defined as the intersection of lines in log–log displacement spectral plots, connecting the flat low-frequency part of the spectrum and a power-law fall-off at high frequencies. However, observed source spectra generally appear more complicated than this simple model, leading to considerable ambiguity in how to measure fc. Even with the synthetic data for a simplified fault geometry presented here, the value of fc for a given spectrum will vary depending upon the selected bandwidth, the choice of fitting function, and the norm or weighting used for the fit. Given these variations, it is difficult to define the ‘best’ way to measure fc. However, it is probably more important for seismologists to agree on a standard way to compute corner frequency, at least from the point of view of facilitating comparisons among different studies. Thus, we have emphasized the Brune spectral function here, not because it necessarily provides the best fit, but because it has been the most widely used and thus provides a logical reference fc value. Of course, whatever method is used to compute fc, it is important to compute stress drop and other parameters using k values that are appropriate for that method.
5.7 On estimation of stress drop using models with stress overshoot
Strictly speaking, one cannot simply use Eshelby's analytical solution (eq. 1) for both the Madariaga and cohesive-zone models because they have ‘spatially non-uniform’ stress drop. Eq. (1) is valid only for ‘spatially uniform’ stress drop. One can set up a cohesive-zone model that does not allow for stress overshoot by enforcing the final shear stress to be equal to the dynamic shear strength (τf = τd) within the source area. In this case, the final slip and stress drop are consistent with Eshelby's analytical solution (eq. 1). However, because of the fixed shear stress on the fault, the sign of the slip reverses (i.e. the fault slips backwards) during the healing phase of the slip. This is inconsistent with the assumption of the source because the backward slip should only occur when the shear stress becomes more negative than the shear strength (i.e. τ < −τst).
6 CONCLUSIONS
Using simulations of dynamic rupture and the representation theorem of Aki & Richards (2002), we have analysed the behaviour of source spectra derived from the cohesive-zone model of a radially expanding rupture with velocity 0.9β on a circular fault. We have found similarities and differences in the behaviour of the source spectra compared to that of Madariaga (1976). The most important finding in this study is that, in the small-scale yielding limit where the cohesive-zone size is much smaller than the source dimension, P- and S-wave corner frequencies of displacement spectra are about 20 per cent larger than those predicted by Madariaga (1976). In this case, an application of the Madariaga model overestimates stress drops by a factor of 1.7. As a result, the seismic radiation efficiency is underestimated by the same factor. The ratio of the P-to-S corner frequencies is found to be 1.5, consistent with that of Madariaga (1976). However, the S-to-P energy ratio is about 22, much higher than predicted by simplified models with identical high-frequency P and S fall-off rates. The spherical average of corner frequencies over the focal sphere is larger for source models with a larger rupture speed due to a shorter source duration. The large dependence of corner frequency on take-off angle relative to the source suggests that measurements from a small number of seismic stations are unlikely to produce unbiased estimates of spherically averaged corner frequency. While the relative variations of inferred stress drops based on the Madariaga model reported in observational work are not affected by the result of this study, the absolute levels of the stress drops were likely overestimated.
We thank Ralph Archuleta, Harsha Bhat, Nadia Lapusta, Jean-Paul Ampuero, Rafael Benites, Anna Kaiser and Bill Fry for helpful discussion. The paper greatly benefited from comments by Raul Madariaga and Eric Dunham. We wish to acknowledge the contribution of NeSI high-performance computing facilities to the results of this research. NZ's national facilities are provided by the NZ eScience Infrastructure and funded jointly by NeSI's collaborator institutions and through the Ministry of Business, Innovation, and Employment Research Infrastructure programme. This research was supported by GNS Science and the Southern California Earthquake Center (SCEC). SCEC is funded by NSF Cooperative Agreement EAR-1033462 and USGS Cooperative Agreement G12AC20038. The SCEC contribution number for this paper is 1830.
REFERENCES
APPENDIX A: NUMERICAL DAMPING: KELVIN–VOIGHT VISCOSITY
The Kelvin–Voight viscosity we use is artificial in the sense that it is part of the numerical procedure for solving the perfectly elastic problem, and is not intended to represent physical damping. The purpose of the Kelvin–Voigt viscosity is to damp spurious oscillations generated by the fault slip at frequencies that are too high to be resolved by the mesh (e.g. Day et al.2005). The viscosity η depends on the size of the elements on the fault and must be a small fraction of the critical time step in an elastic medium for the average linear size of the elements on the fault plane (e.g. Ampuero 2002; Kaneko et al.2008). Unless noted otherwise, we set η = 0.1Δt, where Δt is the time step used in a simulation. Values substantially larger than this value visibly degrade the sharpness with which shear stress is resolved at the rupture front.