G A kinematic rupture model generator incorporating spatial interdependency of earthquake source parameters

SUMMARY Based on the results obtained by analysing a large number of dynamic strike slip and dipping fault ruptures we construct a kinematic rupture model generator that incorporates some of the key source parameters extracted from the dynamic rupture models. In this kinematic rupture model, the slip rate function includes the ﬁnal slip, the rise time, the local rupture velocity and the peak time (a measure of the impulsive part of the slip rate function). A four-dimensional correlation matrix is used to describe the spatial interdependency between the four source parameters—each of which is modelled as correlated random ﬁeld. Each source parameter has a different marginal distribution determined from the analysis of the dynamic models. The marginal distributions are allowed to change as a function of distance from the hypocentre. The autocorrelation of each parameter is modelled by a power spectrum that follows a power law. The value of the spectral decay of the power law for the different source parameters is based on the results obtained from the dynamic rupture models. Finally, the values of rise time and peak time are adjusted such that the moment rate function ﬁts a Brune spectrum for a speciﬁed corner frequency. We validate the rupture model generator using observed strong motion near ﬁeld recording for the 1994 Northridge earthquake and the 1989 Loma Prieta earthquake. We perform comparisons with the rupture model generator of Liu et al . We ﬁnd that only considering the response spectral bias of a best model is not sufﬁcient to validate a rupture model generator. Thus we test the predictive power of the two rupture model generators if multiple ruptures are computed. Overall, the new rupture model generator yields a better prediction, compared to Liu et al. , for the two validation events, especially in predicting observed PGV values and pseudospectral velocity at low frequencies ( < 1 Hz).


I N T RO D U C T I O N
Because recordings in the near field of large earthquakes are sparse, earthquake simulations become important for seismic hazard assessment. To achieve reasonable predictions of ground motion, the earthquake process has to be understood well enough to give a reasonable approximation of the underlying physical mechanisms.
Two approaches are used to model earthquakes: dynamic models and kinematic models. A dynamic model assumes the initial stress on the fault and a friction law that describes the evolution of stress on the fault. Thus, a dynamic model of an earthquake uses stress boundary conditions coupled to a constitutive law for friction (e.g. Day 1982). The earthquake is modelled as a propagating shear crack (Kostrov 1964(Kostrov , 1966. The advantage of dynamic models is that in addition to the ground motion, the spatio-temporal evolution of slip on the fault is also computed. The disadvantage of dynamic models is the computational cost, which makes simulation of broad-band seismograms (0-10 Hz) for large magnitude earthquakes in realistic 3-D velocity structures very expensive at the current time. Furthermore, 3-D velocity and attenuation structure is typically not known for short wavelengths.
A computationally feasible description of an earthquake is a kinematic model (Haskell 1964). A kinematic model uses a slip boundary condition, that is, the complete spatio-temporal evolution of slip on the fault is prescribed. Because there is no interaction between different points on the fault (called sub-fault in the following), it is easier to implement than a dynamic model. However, a kinematic description requires the implementation of the spatiotemporal evolution of slip on the fault. Hence it is crucial to gain a better understanding of the spatio-temporal evolution of slip during an earthquake.
Several formulations have been discussed in the literature. For example, Guatteri et al. (2004) construct a pseudo-dynamic description based on their statistical analysis of a few dynamic rupture models. One of the important features incorporated in their model consists in deriving the rupture velocity from the fracture energy. The derivation is based on results obtained in 2-D (Andrews 1976b). Liu et al. (2006) construct a kinematic rupture model that assumes a spatial correlation between the final slip and the average rupture velocity, and a spatial correlation between the final slip and the rise time. Frankel (2009) develops a kinematic model based on the idea of constant stress drop; he also assumes a correlation between final slip and rupture velocity. Song & Somerville (2010) propose a model that takes into account correlations at non-zero offsets. Graves & Pitarka (2010) propose a hybrid model: (deterministic) at low-frequencies (f ≤ 1.0 Hz) a specified slip-rate function is summed over the fault; a random phase (stochastic) component is computed by assigning a stress drop to each subfault and summing their contribution in the frequency domain; after transforming to the time domain the high frequency stochastic contribution is added to the low-frequency deterministic contribution using matched filters. Mena et al. (2010) use a similar approach for the low frequencies.
Instead of enforcing a random phase they convolve an S-to-S scattering function with the discrete arrivals of a 1-D layered medium (see also Mai et al. 2010). The Fourier amplitude of the scatterogram is added to the Fourier amplitude of the low-frequency synthetic and the total is transformed to the time domain.
Kinematic inversions of the earthquake rupture process can provide some information about the rupture process, for example, about the spatial distribution of slip on the fault (see Lavallée 2008, for a summary). However, resolution for the rupture velocity is poor because changes in the rupture front radiate frequencies higher than typically utilized in the inversion process. Furthermore, there are trade-offs between the different parameters being inverted for. Hence, dynamic models of the earthquake rupture process are necessary to provide enough data that can be used to constrain the source parameters under quasi-ideal conditions. Schmedes et al. (2010a) (from now on SAL10) analysed a large (>300) number of heterogeneous dynamic strike slip ruptures; they limited their analysis to those with subshear rupture propagation. To avoid bias arising from use of a single model, they account for epistemic uncertainty by using various methods to construct initial models. In SAL10, earthquake ruptures are computed for different sets of conditions including different models for the spatial correlation of stress on the fault, homogeneous or heterogeneous modelling of the peak stress, and also homogeneous or heterogeneous modelling of the slip weakening distances. In addition, they analysed models computed by other authors (Dalguer et al. 2008;Olsen et al. 2009).
The results of the study by SAL10 allow one to infer a correlation matrix accounting for the spatial interdependency between the source parameters. It is important to note that SAL10 shows that there is no correlation between slip and rupture velocity, as assumed by other studies Frankel 2009;Graves & Pitarka 2010). SAL10 find no correlations at non-zero offset. Because this result differs from that presented in Song & Sommerville (2010), it will be worthwhile to revisit their result to explain those discrepancies. SAL10 also compute the marginal distributions of the source parameters as a function of the distance from the hypocentre. They find that as the rupture front extends from the hypocentre, the marginal distribution of the rupture velocity becomes more compact around a single rupture speed, that is, the rupture accelerates up to a limiting speed. According to SAL10, the rupture velocity depends on the fracture energy and the slope of the assumed slip weakening law.
In addition to the correlation matrix and the marginal distributions provided by SAL10, we also compute the power spectral decay (PSD; describing the autocorrelation) for a subset of models. All findings are implemented in a new kinematic rupture model generator.
While SAL10 only looked at strike slip ruptures, we have also analysed a dynamic rupture on a dipping fault and compare our results to the strike-slip ruptures. The rupture we have analysed is a dynamic model for the 1994 Northridge event (Ma & Archuleta 2006). We validate the new rupture model generator against the strong motion data from the 1994 Northridge and 1989 Loma Prieta earthquakes. We also compare our validation results with those obtained using the rupture model generator discussed by Liu et al. (2006).

B A C KG RO U N D
To construct a kinematic model of the earthquake source, the fault area is divided into N sub subfaults d i . For each subfault the evolution of slip is prescribed. In this work we describe for each subfault a slip rate function and the final slip u i .
We use the same slip rate function as SAL10; the slip rate function is illustrated in Fig. 1.
This slip rate function is computed using a convolution between a Yoffe function (Nielsen & Madariaga 2003) and a half-sine function (Madariaga, Personal Communication 2008, Eq. 1) for t ≥ 0. This is very similar to the function used by Tinti et al. (2005Tinti et al. ( , 2009. The function H(t) is the Heaviside function and the constant A is given by imposing the condition ∞ −∞ṡ (t, T ) dt = 1. The parameters T = (T 0 , T r , T p ) necessary to describe the slip rate function are rupture time T 0 (start of slipping, related to rupture Figure 1. Source time function used in SAL10 and in this study. It is computed using a convolution between a Yoffe function (left-hand side) and a half sine. The width of the half sine is called peak time in this study. The rise time describes the total duration of non-zero slip velocity in the slip rate function (right-hand side). Figure is from SAL10. velocity), rise time T r (duration of slipping) and the peak time T p (duration of impulsive part of slip rate function, see Fig. 1). The absolute slip rate functionṠ i (t) for each subfault d i is then given aṡ (2) Thus, given the functional form of the slip rate function defined in eq. (1), there are four kinematic source parameters that are distributed heterogeneously on the fault and must be specified for each subfault: (i) total slip; (ii) rupture time; (iii) rise time; (iv) peak time. The geometry of each subfault-strike, dip and rake-must also be defined.
With these parameters specified for each subfault, one can compute the ground motion at a given site. For simplicity we will only consider one component of motion.
For each subfault d i one can compute its contribution dU i (t) to the ground velocityU (t) at a given location X . This requires that the Green's functions G for the subfault and the location X are known. The contribution dU i (t) is computed by convolvinġ S i (t) with g i (t). The function g i (t) is computed using the Green's functions G and the radiation pattern of the subfault. The ground velocity is computed by summing dU i (t) over all subfaults (eq. 3). The symbol * denotes a convolution.
Each slip rate function is convolved with an appropriate (radiation pattern, component of motion) impulse response of the medium between the subfault and the observer location X. All of the filtered slip rate functions are summed to give the ground velocityU (t) resulting from the slip evolution over the whole fault .

Results from analysis of dynamic strike-slip ruptures
In this section, we provide a short summary of the results obtained by SAL10 that are important for this study. In SAL10, the parameters of the slip rate function are computed by fittingṠ i (t) (see Fig. 1) to the time dependent slip rate values generated with a dynamic model of an earthquake.
The spatial interdependency of the source parameters ( u,T 0 , T r ,T p ) is described using a covariance matrix C, which can be used to derive the correlation matrix ζ . SAL10 compute a distribution Z of correlation matrices, that is, one correlation matrix for each dynamic rupture analysed. This distribution will be useful because later we can compare the correlation matrix ζ obtained for a dipping fault rupture with Z.
For the kinematic rupture model generator discussed here, we use one covariance matrixC and one correlation matrixζ that were computed using the combination of all analysed ruptures (Fig. 2). Only regions on the fault where the local rupture velocity was subshear were considered in SAL10. Hence, at this point, the rupture model generator will be valid only for subshear rupture propagation. We will discuss the extension to supershear rupture velocities and the related challenges in a separate section.
SAL10 show that the computed correlations are robust and do not change if smaller subsets are used to compute them. As mentioned before no correlation exists between final slip and local rupture velocity. This is important because if such a correlation existed, one would expect larger peak ground motions (Schmedes & Archuleta 2008). Slip correlates with rise time, as assumed by Liu et al. (2006). The local rupture velocity correlates positively with the peak slip rate and negatively with the peak time. Peak slip rate and peak time are determined during the breakdown process, that is, during the time the instantaneous stress drops from its peak strength to sliding friction.
The marginal distributions of local rupture velocity, rise time and peak time depend on the autocorrelation of initial stress (see also Schmedes et al. 2010b). Furthermore, the marginal distributions have a dependence on distance from the hypocentre. As the rupture propagates along the fault, it accelerates, and the rise times and peak times become shorter. Consequently, the peak slip rates also increases as the distance from the hypocentre increases. This raises the question if there is also a dependency of the computed correlations with distance from the hypocentre. This question was not investigated in SAL10 but is important for this study.  . The correlation coefficient is computed for different parameter pairs and different distances from the hypocentre. Only ruptures with a power spectral decay of −2 in stress are used. The first distance bin contains points on the fault that are 0-10 km away from the hypocentre, the second points that are10-20 km away, etc. No dependency on distance is noted. The correlation matrix formed by using the mean values is used as the target correlation matrix. An iterative scheme is used to create correlated random fields with the desired target correlation matrix. Note, that there is almost no correlation between slip and rupture velocity. Rupture velocity correlates best with the peak time, which is a measure of the duration of the positive slip acceleration. The rupture velocity is controlled during the breakdown process.
To examine this possible dependence we compute the correlation matrix for 10 km wide distance bins, that is, for all subfaults 0-10 km from the hypocentre, 10-20 km from the hypocentre, etc. We use only those ruptures that have a PSD of 2 in the initial stress, that is, ruptures with a power spectrum in the initial stress that follows k −2 , where k is the radial wavenumber. The reasons for this choice will be discussed in the next section. Fig. 3 shows the correlation matrices for different distance bins; there is no obvious distance dependence. There is a slight difference for distances closest to the hypocentre; this difference might be due to the bin being close to the nucleation location.
The marginal distributions for the different distance bins are shown in Fig. 4. In the following we will use these marginal distributions, which are computed for ruptures with the same autocorrelation in initial stress.

Analysis of autocorrelation
We compute the PSD as a measure of the autocorrelation. We derive the PSD only for the ruptures with an initial stress that have a wavenumber power spectrum that follows a power law. The computation is performed in 2-D. To evaluate the exponent ν ps of the PSD we fit eq. (4) to the computed power spectrum for each parameter of interest. In our analysis, the PSD is uniquely determined by the exponent ν ps . Thus, the notations PSD and ν ps are interchangeable.
To investigate dependencies of the autocorrelation between two different parameters x and y, we plot the PSD exponent of one parameter versus the PSD exponent of the other parameter (Fig. 5). Note, that a larger value of ν ps visually translates into a stochastic field that looks smoother because the contribution of the high wavenumbers (small wavelength) is reduced.
For initial stress and final slip we find that approximately That is, final slip is more strongly correlated spatially (appears 'smoother') than stress. The Fourier spectral decay (FSD) for a parameter can be computed as ν Fs = ν ps (x) 2. Hence we find that ν Fs (stress) ≈ ν Fs (slip) − 1 as proposed by Andrews (1980).
The rise time shows a dependency with the final slip; we find approximately ν ps (rise) ≈ 0.54ν ps (slip) + 1.35.
The rupture velocity and peak time only weakly depend on the autocorrelation of initial stress (or slip). They always have a small . Probability density of rise time, rupture velocity divided by shear wave velocity, and peak time are shown. Different distance bins from the hypocentre are used. As the rupture propagates along the fault it accelerates and the rise times and peak times get shorter. This distance dependency is implemented in the kinematic rupture model generator. However, the correlations themselves do not depend on distance from the hypocentre (Fig. 3).

Figure 5.
Plotted is the slope of the 2-D power spectrum for different source parameters versus −ν ps the negative power spectrum decay for slip. Slip and initial stress show a dependency that is approximately consistent with Andrews (1980). The PSD of rupture velocity ratio (rupture velocity divided by shear wave velocity) and peak time only weakly depend on PSD of slip amplitude, while rise time shows a stronger dependency. exponent for the PSD, that is, they both appear spatially rough. We find approximately ν ps (vrup) ≈ 0.23ν ps (slip) + 0.74 ν ps (peaktime) ≈ 0.29ν ps (slip) + 0.83 As shown by SAL10 there is a dependency of the marginal distributions on the initial stress. Hence, we need to decide on one autocorrelation for slip in our kinematic rupture model generator. Analysis of subshear and supershear ruptures suggests that the autocorrelation in the Earth probably has a PSD in the range of 1 < ν ps ≤ 2 (Schmedes et al. 2010b). For this study, we chose a PSD of v ps (stress) = 2, which corresponds to a PSD in slip ν ps (slip) = 4 and FSD ν Fs (slip) = 2. This corresponds to the 'k-squared' model. The high wave number behavior of slip distributions is controlled by the discontinuities of slip, a crack like discontinuity in 3D producing a k-squared distribution (Madariaga 2011) Andrews (1980 showed that this model, which corresponds to a FSD of 1 in stress, has an average stress drop that is independent of length scale (thus magnitude). Table 1 lists the exponent of the power spectrum for the other variables in Fig. 5 corresponding to a model with ν ps (slip) = 4. These are the values, which will be used in the following sections.

2-D correlation and supershear propagation
Song & Sommerville (2010) have proposed accounting for non-zero offset correlations. For example, does large slip correlate with fast rupture velocities further along the fault? After analysing a dynamic rupture model by Dalguer et al. (2008), they conclude that there is a response distance, during which the rupture front accelerates after it passed through an area of large slip. This results in a positive correlation between slip and rupture velocity at a positive offset of about 8 km. In contrast to this result, SAL10 looked at the average 2-D correlation map of 15 dynamic ruptures and did not find any correlation between slip and rupture velocity at any offset.
The key to understanding this apparent discrepancy lies in the data considered in the analysis. First, the model analysed by Song & Sommerville (2010) has regularly spaced asperities of high stress drop and high strength excess. In between these barriers the stress drop is 0 and the strength excess is small. Between the asperities (regions of high slip) the rupture propagates at supershear speed (e.g. Burridge 1973;Andrews 1976;Freund 1979;Archuleta 1984;Bouchon 2001;Bouchon & Vallee 2003;Dunham 2003;Dunham & Archuleta 2004;Ellsworth 2004;Robinson et al. 2006;Dunham 2007;Bizarri et al. 2010). The SAL10 models include heterogeneity at all scales between the size of a subfault to the full size of the fault. SAL10 consider only points that rupture at subshear speed while the distribution of rupture velocities in the model of Dalguer et al. (2008) shows supershear propagation.
The non-zero offset correlation can be explained by the results of Liu & Lapusta (2008). When the rupture front passes from an area of high strength excess to an area of low strength excess, the stress peak in front of the crack front reaches the level of the yield strength and forms a daughter crack. The rupture front does not accelerate; rather a daughter crack is nucleated but vanishes once it hits the next area of high strength excess. This situation is physically different from a pure subshear rupture. Interpreting this situation using the erroneous concept of a single rupture front would produce a nonzero offset in the correlation between rupture velocity and slip.
Suppose the strength was a constant value everywhere on the fault. The supershear transition would occur in the area of high stress drop and high slip-areas with low strength excess. In this situation, which is discussed in Liu & Lapusta (2008), one would get a positive correlation between slip and rupture velocity at zero offset. SAL10 did not observe this correlation when only subshear propagation was considered. Given the complex patterns that can arise with supershear ruptures, we conclude that one should limit the analysis to a single rupture mode (subshear or supershear) when computing spatial coherence.
We computed the correlations matrix using 20 dynamic supershear ruptures. The correlation matrix we obtained is almost identical to the matrix obtained for subshear propagation.
There are many questions remaining about how to include supershear propagation into a kinematic model: (i) what is the correct slip rate function to use? (ii) where does the supershear transition take place? (iii) what are the correct marginal distributions? (iv) how large is the fraction of the fault that propagates at supershear speed?
(v) what is the frequency of occurrence of supershear ruptures? Schmedes (2010) found that for ruptures that have a potential to propagate with supershear speed, about 15 per cent did. For those supershear ruptures, about 40 per cent of the fault, on average, experiences a supershear speed.
The marginal distribution of rupture velocity to use for supershear propagation can be extracted from Schmedes et al. (2010b). It is bimodal, with one subshear mode and a second mode between the Eshelby speed √ 2V s and the P-wave speed. Because of the complexity related to supershear ruptures, we will exclude supershear in our rupture model generator at this time.

Dynamics of dipping fault
To check if the results from analysing strike-slip ruptures are also applicable to dipping faults, we analyse a dynamic rupture computed for the 1994 Northridge earthquake (Ma & Archuleta 2006). Of course, looking at only one event is not sufficient to validate that the correlations computed for strike slip ruptures hold for the case of dipping faults. Here we just want to be sure that we do not falsify this hypothesis.
We compare the computed correlations with the correlations of all events as discussed in SAL10. We see that the computed correlations fall within the distribution of correlations found for the strike-slip ruptures (Fig. 6). The correlation of 0.27 between slip and rupture velocity is larger than the average correlation of all events but still falls within the distribution. There is a strong negative correlation between peak time and rupture velocity, as observed for the strike slip ruptures. The correlation between slip and rise time is smaller than for most strike-slip ruptures but still within the distribution. With these observations we assume, for now, that the computed correlation matrix,ζ for strike-slip faults, can be used for dip slip faulting. Additional modelling will be needed to reach a definitive conclusion.
The marginal distribution of rupture velocity is expected to be different for a dipping fault. For a Mode III rupture, there is no forbidden range for the rupture velocity between the Rayleigh wave speed and the shear wave speed. Hence, there will be larger local subshear rupture velocities compared to a strike slip rupture (a Mode II rupture). Fig. 7 shows the marginal distribution of rupture velocity for the dynamic model of Northridge. Because the fault area is not large, we do not account for dependence on the distance to the hypocentre. The limiting rupture velocity is the shear wave velocity, as expected for a Mode III rupture. The rupture velocity also has a higher probability for slower rupture velocities than the strike-slip ruptures in SAL10. Of course one sample is not indicative of general behaviour. One possibility for this difference is how the frictional parameters are chosen, that is, is not necessarily a systematic difference between strike-slip and dip slip earthquakes. A larger likelihood for slower rupture velocities can be caused by the assumed autocorrelation (Schmedes et al. 2010b) or, for example, due to a larger number of areas with negative stress drop  (energy sinks). We will use the distribution shown in Fig. 7 for the verification with Northridge records.

Implementation of new kinematic rupture model generator
In the following, we describe how we create four stochastic fieldsslip, rise time, peak time and rupture velocity-that have the marginal distributions shown in Fig. 4 and the correlation matrix shown in Fig. 2 (without the last column and row representing peak slip rate).
In a first step, we construct a 4-D normal distribution. We use the vector of mean values η = (0,0,0,0) and the covariance matrix C, which corresponds to Fig. 2. The resulting four 1-D normal distributed vectors provide the correct amount of correlations between each pair of parameters. In the next step, we map these 1-D vectors onto the 2-D fault plane. We filter each 2-D field in the wavenumber domain to enforce the chosen PSD for each field. This is achieved by multiplying the Fourier spectrum with k −ν FS and then computing the inverse Fourier transform.
Finally, the normal distribution is mapped to the amplitude distributions shown in Fig. 4 using a normal score transform. We do this by replacing the kth largest value of the normal distribution by the kth largest value of the target distribution in each distance bin. To ensure a smooth transition between distance bins, we let the different bins overlap by 2 km and perform a linear transition in this region. For the total slip, we use the same Cauchy distribution (Lavallée et al. 2006) that is independent from the distance from the hypocentre.
Because these stochastic fields were filtered and their amplitude distribution changed, the newly created fields are unlikely to have the desired correlation matrixζ . Hence we use an iterative scheme to adjust the covariance matrix that was used to construct the initial 4-D normal distribution. Convergence is typically achieved in a few iterations. In fact, instead of iterating until we achieve a reasonable close approximation toζ , we always use five iterations. This leaves room for some variation in the correlation matrix for each event, as observed for dynamic ruptures (see histograms in Fig. 6).
To compute the rupture times from the local rupture velocity we use a 2-D finite difference code that solves the wave equation. The local rupture velocity is used as the propagation speed. The rupture times are computed by measuring the first arrival of the wave front at each point on the fault. This is essentially the same procedure used by Frankel (2009).
Finally, we adjust the vectors of rise times R 0 and peak times P 0 such that the moment rate function of the kinematic rupture has approximately the shape of a Brune-spectrum for a given corner frequency, that is, a displacement amplitude spectrum that is flat below the corner frequency with ω −2 decay at high frequencies.
The corner frequency is computed using the rupture times. First, we compute the moment rate spectrum assuming that the slip rate function is a δ-function with amplitude 1.0. Then we compute the cumulative moment rate, and we determine the time t 1 at which 7.5 per cent of the moment rate is released and the time t 2 at which 92.5 per cent of the moment rate is released. The corner frequency we use is then determined as Note, that the high frequency ground motion is very sensitive to the definition of the corner frequency. We determined 7.5 and 92.5 per cent during the validation; we do not know if these are the Downloaded from https://academic.oup.com/gji/article-abstract/192/3/1116/807260 by guest on 28 July 2018 optimum values for all cases. We compute the corner frequency to distinguish between events that have the same magnitude but might have other differing characteristics, such as unilateral versus bilateral rupture (more on this in the section about validation of Loma Prieta).
After the corner frequency is determined, we use the downhill simplex method (Nelder & Mead 1965;Press et al. 1999) to minimize a cost function based on the Brune spectrum and the moment rate spectrum (now computed using the finite slip rate function shown in Fig. 1) and estimate two scaling factors: β rise and β peak . The two scaling factor are used to rescale the amplitudes of rise time and peak time, respectively. That is, we use the vector R = β rise R 0 for rise times and P = β rise P 0 for peak times.
Before we move on to the validation, we contrast our new rupture model generator with the one by Liu et al. (2006), to which we will refer to as LAH in the following. The key differences are as follows: (1) LAH uses a different source time function and only three free parameters: slip, rise time and average rupture velocity. It is important to note, that the source time function in LAH contains a peak time, which is fixed to a correlation of 1 with the rise time.
(2) We use six different correlations (off-diagonals of correlation matrix), while LAH only use a positive correlation between slip and rise time, and slip and average rupture velocity. The latter correlation was not confirmed by SAL10.
(3) We determine the scaling of the peak and rise times by fitting the moment rate function to a Brune spectrum, while in LAH the spectral level at high frequencies is determined using the following relation (eq. 5).
with m i and τ i are respectively the seismic moment and the rise time corresponding to the subfault with index i; M 0 and f c are respectively the seismic moment and the corner frequency of the main event; N sub is the number of subfaults and β the scaling factor for which eq. (5) is solved.
(4) LAH describes the power spectrum of slip on the fault using a von Karman function (Mai & Beroza 2002), while we use a power law (Lavallée 2008).
In the following section we will illustrate how these differences translate in differences in the computed ground motion.

VA L I DAT I O N
To validate our new rupture model generator, we use near field recordings of the 1994 Northridge event, and the 1989 Loma Prieta event. We use the procedure described in LAH to compute the broad-band seismograms for a 1-D velocity model. We chose not to use a 3-D velocity model to be able to compute multiple ruptures. While neglecting 3-D effects might lead to a poor prediction for certain sites, this effect would be the same in LAH and SAL.
In the procedure outlined in LAH, Greens functions for a 1-D velocity model are computed using the FK-method (Zhu & Rivera 2002). The dip, rake and azimuth are randomised above a frequency of 3 Hz to simulate an isotropic radiation pattern (e.g. Pitarka et al. 2000). Between 1 Hz and 3 Hz there is a linear transition from a deterministic radiation pattern to a stochastic one. As mentioned earlier, for Northridge we use the rupture velocity distribution in Fig. 7 that was computed from a dynamic model of Northridge rather than the marginal distribution of rupture velocity (Fig. 4) for the validation.
To account for non-linear site response, the 1-D non-linear code NOAHW (Bonilla et al. 2005) is used to propagate the waves through the upper 300 m. The input to NOAHW depends on the site class of the station (Wills et al. 2000, LAH). Finally, because the 1-D code does not account for surface waves, we use a matchingfilter in the wavelet domain (Pengcheng Liu, Personal Communication 2008) to combine seismograms with linear response up to 1.0 Hz with the seismograms, which have a non-linear response above 1.0 Hz. This wavelet matching-filter is more accurate than a matching-filter in the frequency domain.
To compare the synthetic ground motion with the observed we initially use the response spectral bias (Abrahamson et al. 1990). This bias is computed for a 'best' model based on the average misfit of all stations. The bias with respect to the observation is also affected by the path (velocity model, scattering) and the site (nonlinear effects). The intrinsic and scattering attenuation becomes more important at high frequencies. In the method of LAH there is no frequency dependence for the seismic attenuation. We limit the validation to distances of one fault width because at larger distances such effects can dominate, and it becomes more a validation of the method to compute accurately high frequencies rather than validation of the rupture generator.
We also perform the validation using LAH and compare the results to the results based on the new kinematic description (called SAL for Schmedes, Archuleta and Lavallée in the following).
For SAL and LAH, we compute ground motion for 20 rupture models for each of the two events used in this validation. Using multiple models allow us to look at the variability of ground motion at single stations. Besides the response spectral bias, we also compare the amplitude distributions of peak ground velocity (PGV) and pseudospectral velocities in various frequency bands. Comparing the amplitude distribution is important. Even for small average bias and standard deviation of a single rupture model there is the possibility that the peak ground motion is not predicted well among all rupture models, as shown in the following section. We also test a different autocorrelation and assess its effect on ground motion. Finally, we test the stability of the computed ground motion for different subfault sizes.

Northridge
For this validation we use 44 stations that have a closest distance to the fault plane less than 25 km (Fig. 8). The site classification for each station is determined from Vs30 following Wills et al. (2000). The velocity model defined in LAH is used to compute the linear 1-D Green's functions using the FK-method (Zhu & Rivera 2002).
The fault geometry is the same as in Liu et al. (2006). Figs S1 and S2 in the electronic supplement show an example of a rupture model for Northridge (S1) created using the new rupture model generator and a comparison of ground acceleration for a few selected stations (S2).
In Fig. 9 we show the bias for the models with the smallest bias for LAH and SAL. At low frequencies, the average bias for SAL is smaller. At high frequencies, LAH and SAL have almost identical bias. Overall both yield comparable fits to the observed response spectra. Comparing the models with the minimum bias for both kinematic rupture model generators we find no significant difference.
The fact that a best rupture model exists does not guarantee that the kinematic description is useful for predicting ground motion from future earthquakes. From a statistical point of view, the important question is what difference exists between the models SAL and LAH when multiple ruptures are considered? For example, do the observed PGV at a single station fall within the probability distribution of PGV computed using all rupture models? This is a situation one could encounter when making predictions for a single fault, for example, in the context of probabilistic seismic hazard assessment.
One difference between SAL and LAH is evident if the probability densities of the 3-component PGV are compared (Fig. 10). Consider the distributions for station tpfs and grif. The observed PGV falls within the distribution from SAL but is outside the distribution determined from LAH. This can be quantified by computing the percentile that the observed PGV is within the distribution for each station. More precisely, for a given station we compute the quantity N (PGV obs > PGV syn )/N event where N event = 20 is the total number of rupture models used. Fig. 11 shows the histogram of those percentiles for the 44 stations and the two kinematic descriptions SAL and LAH. For LAH 15 of the stations have an observed PGV close to the 0th percentile (red bar in Fig. 11) meaning that all of the synthetic models over-predict the observed PGV. For eight stations the observed PGV is close to the 100th percentile (red bar in Fig. 11) meaning the observed PGV are under-predicted by all 20-rupture models. For SAL only nine stations have PGVs that are completely over-predicted (three stations) or under-predicted (six stations) compared to 23 stations for LAH. It is important to note, that overprediction of all models for a given site does not imply the computed ground motions are unusually large. It could be, for example, caused non-linear site effects, which are not reproduced in the modelling procedure or it could be that the motions are only slightly greater than the observed. We have not examined the degree to which a prediction is greater or smaller than the observed.
A similar approach can be used to compare the predictions of the two rupture models for the pseudo-velocity spectra (root mean square of three components). This approach allows us to examine the frequency dependence of the predictions. We consider the following frequency bands: (i) 0.0-1.0 Hz; (ii) 1.0-5.0 Hz; (iii) 5.0-10.0 Hz; (iv) 10.0-15.0 Hz. For each frequency band we compute the percentile of the observed spectral value to fall within the 20 synthetic pseudospectral values for a given frequency. We plot histograms that combine the results for all stations and all the frequencies in the given frequency band.
We compare results from SAL and LAH in Fig. 12. As seen in the minimum bias, SAL performs better than LAH in the lowest frequency band. For higher frequencies neither method does all that well. LAH is worse in that the total probability of either overpredicting or under-predicting is larger for LAH.
One striking difference between the two methods is that SAL shows significantly more under-prediction than over-prediction. This could be caused by a corner frequency that is too small or by a rupture velocity that is too slow. It is also possible that the path and site effects may be responsible for this underprediction. This is reasonable in that we do not observe this strong under-prediction at low frequencies for which the computation of ground motion is completely deterministic (LAH). Furthermore, the under-prediction at high frequencies is about the same for LAH and SAL. One reason could be the seismic attenuation, which especially influences the higher frequencies. The method of LAH does not account for frequency dependent seismic attenuation.
Another observation is that the under-prediction is strongest between 1.0 and 5.0 Hz for LAH and SAL (Fig. 12). One possibility is the non-linear site response. To assess the influence of the site response, we show the bias for the best model before and after the site response is included, that is, we compare the completely linear solution with the solution for which the non-linear wave   Figure 11. Histogram of the quantity N(PGV obs > PGV syn )/N event computed for all stations (N event = 20, PGV obs and PGV syn are vector sum PGV for observed and synthetic data, respectively). The red bar at 0th percentile gives the number of stations where the rupture generator overpredicts the observed PGV. The red bar close to 1 (100th percentile) gives the number of stations where the rupture generator underpredicts the observation. SAL performs significantly better than LAH in predicting the observed PGV. Figure 12. We show the probability density of N(PSV obs (f i ) > PSV syn (f i ))/(N event * N freq ) of observing a three-component vector sum of pseudospectral velocity for four different frequency bands. The left column result from 20 ground motions produced by LAH; the right column are results from SAL. The rupture generator LAH tends to overpredict (0.0 percentile) and underpredict (1.0 percentile). While this leads to a small average bias, the response spectrum at most stations is not well predicted. SAL is more likely to underpredict the spectral level.
propagation in the upper 300 m was applied. We compare the average bias and the standard deviation of the bias for SAL (Fig. 13). We find that after the site response is applied, the ground motion is smaller (positive bias) in the frequency range 1.0-5.0 Hz, whereas, the ground motion is larger (smaller bias) in the frequency range above 5 Hz. This explains the observed pattern with stronger underprediction for frequencies between 1.0 and 5.0 Hz for both, LAH and SAL.
Path and site effects are most likely responsible for this underprediction of the pseudospectral velocity in the high frequencies. We see the same pattern of under-prediction for the two different kinematic descriptions, which makes it unlikely to be caused by the source. The over-prediction for LAH, also observed in the lower frequencies, is most likely attributed to the source model, in particular, to the assumed correlation between slip and average rupture velocity.
Note that with almost the same probability LAH over-predicts and under-predicts the pseudospectral velocity for the two highest frequency bands. This implies that there is an equal number of stations that have been over predicted and under predicted. This rough equality has implications for the spectral bias. The average bias appears fine even though most stations have observed values outside the simulated probability densities. The small average bias (Fig. 9) results from the trade-off of stations that were over-predicted and those that were under-predicted stations. Judging the best model based on an average measure can be misleading.

Variation of autocorrelation, and subfault size
Using Northridge, we performed several tests to understand how the solution and the variation among models are influenced by the input variables. For SAL we change the exponent of the PSD of slip to ν ps (slip) = 3. The values for the other parameters are given in Table 1 (see also Fig. 5). With this exponent we examine a model that appears 'rougher' spatially. For LAH, we keep the Figure 13. We compute the average bias and its standard deviation for best model using linear wave propagation (solid lines) and using 1-D non-linear wave propagation in the upper 300 m (dashed lines). Linear wave propagation produces a more constant bias for all periods. Non-linear wave propagation has more variability with period. Both approaches tend to underpredict for periods greater than 0.17 s. Figure 14. Histogram of the quantity N(PGV obs > PGV syn )/N event computed for all stations (N event = 20, PGV obs and PGV syn are vector sum PGV for observed and synthetic data, respectively). The red bar at 0th percentile gives the number of stations where the rupture generator overpredicts the observed PGV. The red bar close to 1 (100th percentile) gives the number of stations where the rupture generator underpredicts the observation. For a smaller power spectral decay SAL tends to underpedict more. same PSD at high wavenumbers, but we use a power law (Lavallée 2008) rather than the von Karman function (Mai & Beroza 2002). This change primarily affects the small wavenumbers, which have a higher correlation with the power law.
With these changes we find about the same bias and standard deviation for the best model as we had before making the changes. The percentiles of the observed PGV are shown in Fig. 14. For LAH the histogram looks very similar to the case with von Karman. For SAL, there are more stations that are under-or over-predicted than we had when the PSD had an exponent of 4.0. In fact, there is a greater tendency to under-predict when the PSD has an exponent of 3.0. This under-prediction can be explained. For a smoother model ν ps = 4 the asperities become larger in spatial extent and tend to produce larger ground motion. Schmedes et al. (2010b) observed the same tendency in dynamic models, that is, larger ground motions for spatially smoother models. Smoother models lead to smoother isochrones and thus larger stopping phases.
Another important factor that could influence the outcome is the subfault size. Of course, if a coarsely sampled model is compared to a finely sampled model, one can expect differences, especially at the higher frequencies because more information is available at high wavenumbers for the finely sampled model. Overall, the spectral levels at high frequencies should not differ on average. To test this, we modified the subfault size in each dimension. For LAH, we create a finer model by using 256 × 256 subfaults instead of 128 × 128, while for SAL; we use a grid spacing of 100 m in each direction instead of 200 m. A coarse model is created using 64 × 64 subfaults for LAH, and 300 m grid spacing in each direction for SAL.
The average pseudo-velocity of 20 models for the three different subfault sizes and for two selected stations is shown in Fig. 15. For LAH, the coarse and original model show very similar amplitudes. However, the model with the smallest subfault consistently produced larger amplitudes at short periods for both stations and all the three components. For SAL, we cannot find any systematic changes. There are differences between the three sets of results, but the differences depend on frequency and the component of motion.
The reason for the systematically larger high frequency amplitudes for LAH is related to the scaling of the rise times in LAH. In Liu et al. (2006) the spectral level at high frequencies is determined using eq. (5).
The derivation of the scaling factor β using eq. (5) is not the problem. However, because the duration of positive slip acceleration scales directly with the rise times (100 per cent correlated) and with the definition of the slip rate function adopted by LAH, the duration of positive slip duration is also scaled by the factor β. This yields very short peak times for short rise times leading to large peak slip rates and consequently larger ground motion at high frequencies.
For SAL, the rise times also get shorter for smaller subfault size, but the peak times do not change significantly. Thus, high frequency motion is not larger for the model with the smallest subfault. This result illustrates the importance of introducing the additional parameter peak time in our formulation of kinematic modelling as well as the new approach of fitting the whole moment rate spectrum to a Brune-spectrum. Unfortunately, this comes at the cost of greater computation times for each kinematic rupture model.

Loma Prieta
We use 13 stations that have a closest distance to the fault of less than 20 km (Fig. 16). The Greens functions are computed using the 1-D velocity structure given in Wald et al. (1991). The Vs30 values of the stations are extracted from the PEER strong motion database. Like Northridge the site class for the non-linear wave propagation in the upper 300 m is chosen according to Wills et al. (2000).
We use a PSD of slip that has ν ps (slip) = 4; the exponent for the PSD of the other input variables is chosen according to Fig. 5. We model the fault as 41.25 km long and 20 km wide. The position of the hypocentre is 20 km along strike and 16.5 km along dip at a depth of 18 km with a dip of 70å. This rupture is bilateral with a rupture distance of about 20 km in each direction. Its rupture duration is similar to the Northridge event. For this reason we initially chose for LAH the same corner frequency 0.15 Hz that we used for Northridge. We use a seismic moment of 3.0 × 10 19 Nm.
With this corner frequency LAH completely over-predicts the high frequencies. The over-prediction by LAH is probably due to the way the rise times are scaled, as described earlier. The LAH rise times for Loma Prieta are too small for the assumed corner frequency. The effect is exacerbated because the duration of positive acceleration of slip is chosen as 0.13τ i (τ i is the ith subfault riste time in Eq. 5), yielding a too short duration of positive slip acceleration as well. In contrast to LAH, SAL fits the whole spectrum with different scaling factors for rise time and peak time. We adjusted the corner frequencies for LAH until we found a good bias. We determined a corner frequency of 0.1 Hz instead of 0.15 Hz while increasing the seismic moment to 3.7 × 10 19 Nm.
We also tested the predictive power of the rupture model generators using the observed and computed PGV values for the 13 stations and 40 rupture models, as outlined for Northridge. SAL yields a reasonable prediction of the observed PGV ( Fig. 17) with four stations being over-predicted. With a corner frequency of 0.1 Hz LAH tends to overpredict about half the stations (Fig. 17).

D I S C U S S I O N
When computing the best rupture models generated by SAL and LAH under different conditions, the best models do not show significant differences in bias and standard deviation. However when the study is extended to a larger number of realizations generated by the two models, the distribution of PGV and the distribution of pseudo-velocity spectral amplitudes at different frequencies show significant differences between the SAL and LAH. It illustrates that considering only the response spectral bias for one best rupture model is insufficient to properly validate a kinematic model. This criterion contains no information about the predictive power of the kinematic description. If prediction of ground motion for a future earthquake were the objective, one would not compute only one event. There must be multiple events to account for the aleatoric uncertainty, that is, natural variability. The objective is computing the distribution of potential ground motions. Such a distribution should include the ground motions that will be recorded in future earthquakes and assigned a probability to observe these events. A proper procedure to validate a kinematic model will insure that the observed ground motion parameter (such as PGV or pseudo-velocity spectral amplitude) falls within the distribution of computed ground motion from multiple synthetic ruptures.
Based on all the results of the validation, we conclude that the new model SAL introduced in this study performs better than LAH. The shortcomings in the prediction of the high frequencies are likely caused by path and site effects and not the source. Improving the modelling of path and site is beyond the scope of this study, which is solely concerned about the source description.

C O N C L U S I O N
A new kinematic rupture model generator was constructed based on the analysis of more than 300 dynamic rupture models (Schmedes et al. 2010a). The slip rate function is a Yoffe function (Nielsen & Madariaga 2003) convolved with a half sine (Fig. 1), similar to the slip rate function in Tinti et al. (2005Tinti et al. ( , 2009. Four source parameters are necessary to construct the slip rate function: total slip, rise time, peak time and local rupture velocity. The spatial distributions of the four parameters over the fault are heterogeneous. Their spatial Figure 17. Histogram of the quantity N(PGV obs > PGV syn )/N event computed for all stations (N event = 40, PGV obs and PGV syn are vector sum PGV for observed and synthetic data, respectively). The red bar at 0th percentile means those stations over-predict the observed PGV. The red bar close to 1 (100th percentile) means those stations under-predict. SAL performs better in predicting the observed PGV. interdependency is modelled using a 4-D correlation matrix (Figs 2  and 3).
Each parameter has a different marginal distribution (Fig. 4). The distributions are based on the analysis of the dynamic rupture models. The autocorrelation of each field is modelled using a power density spectrum that follows a power law. For each parameter, a different decay (exponent of the power law) is chosen. Those decays are also based on the interdependency found by analysing the dynamic models (Fig. 5). Our analysis shows the total slip and rise time spatial distributions are smoother than rupture velocity and peak time spatial distributions. The kinematic rupture model generator SAL introduced in this paper is very flexible. The correlation matrix, the marginal distributions and/or the autocorrelation can be modified if needed. It can be easily adjusted as new knowledge becomes available. For example, it could be tailored for a given fault, for which information about marginal distributions or correlations are available. At present the kinematic description is valid only for subshear rupture propagation. Too many open questions exist with regard to supershear propagation, for example, what is the frequency of occurrence of supershear; what is the average area of the fault that propagates at supershear; where does the transition to supershear happen; what is the correct slip rate function to use; how exactly does the transition happen (two cracks)? While some of these questions can be addressed with dynamic modelling, others need more information from observations and inversions.
We validated the rupture model generator against the 1994 Northridge and the 1989 Loma Prieta strong motion data. Only near field recordings were selected for the validation. At larger distances, path effects, especially the seismic attenuation, have a greater influence on the signal. We find that both models yield similar response spectral bias. The difference between the two models becomes evident when looking at the distribution of PGV for multiple models. If multiple models are considered, we find that the new rupture model generator performs significantly better than the method of Liu et al. (2006). Our analysis also suggests that looking at the average bias for a 'best' model can be misleading. For example, a small average bias can be achieved by trading off an overprediction with an underprediction.
There are several reasons why SAL performs better than LAH. In SAL there are six correlations (off-diagonal elements of symmetric correlation matrix) while there are only two in LAH. SAL uses a slip rate function that is more comparable to that for full dynamics compared to the slip rate function used in LAH. SAL uses an additional parameter, the peak time. In LAH the equivalent parameter is the duration of the positive acceleration of slip. This is a fixed parameter −0.13 times the rise time-implying a correlation of 1.0 between the two. This direct scaling of the duration of positive slip acceleration creates a dependency of the high frequency response on the subfault size (Fig. 15).
In the kinematic rupture model generator, we have to select a value for the parameter that governs the decay of the slip wavenumber power spectrum. With this value the autocorrelation for rise time, peak time and rupture velocity can be computed (Fig. 5). We chose a value of 4.0 for the PSD of slip for the validation using data from the Northridge earthquake. We examined a value of 3.0 as well. For Northridge, an exponent of 3.0 produced worse results than an exponent of 4.0. For the validation of the method using strong motion data from the 1989 Loma Prieta earthquake we used a PSD with an exponent of 4.0 for the slip.
We did not consider PSD values larger than 4.0. In analysis of supershear transition Schmedes et al. (2010b) concluded that values of the PSD for the initial stress of 3.0 and 4.0 (5.0 and 6.0 for slip, respectively) resulted in an exceptionally high probability for supershear rupture propagation. Schmedes et al. (2010b) show large values of the exponent often yield PGV values significantly larger than those observed, that is, larger than about 3.0 m s -1 observed during the 1999 Chi-Chi earthquake (Shin et al. 2000).
By using more than 300 dynamic simulations to determine correlations among source parameters we have derived a new rupture model generator for broad-band ground motion. We validated this generator's ability to produce realistic ground motion against data from the Northridge and Loma Prieta earthquakes. The bias and standard deviation for the 'best' model is equivalent to LAH. We show that measures of the recorded motion (PGV, spectral amplitudes) generally falls within the distribution of ground motion computed from a suite of randomly generated rupture models. This new method can be used to predict the range of ground motion one might anticipate from a future earthquake.

A C K N O W L E D G M E N T S
We thank Pengcheng Liu for the code used for LAH and the modelling of broad-band ground motion. We thank an anonymous reviewer and Martin Mai for valuable comments, which improved the quality of the paper. This work was supported by SCEC, a gift from PG&E and the National Science Foundation grant EAR0738954. SCEC is funded by NSF Cooperative Agreement EAR-0529922 and USGS Cooperative Agreement 07HQAG0008. This is SCEC contribution 1681.

S U P P O RT I N G I N F O R M AT I O N
Additional Supporting Information may be found in the online version of this article: Figure S1. Rupture model generated using the new rupture model generator and the fault geometry of the 1994 Northridge earthquake. Figure S2. Ground acceleration computed using the rupture model shown in Fig. S1. The simulated ground acceleration for each station is compared to the one observed during the 1994 Northridge earthquake. Note, that the rupture model was generated randomly, that is, no inversion was performed. Furthermore, the greens functions were computed for a 1-D earth (http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggs021/-/ DC1) .
Please note: Oxford University Press are not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.