Variability in synthetic earthquake ground motions caused by source variability and errors in wave propagation models

Numerical simulations of earthquake ground motions are used both to anticipate the effects of hypothetical earthquakes by forward simulation and to infer the behaviour of the real earthquake source ruptures by the inversion of recorded ground motions. In either application it is necessary to assume some Earth structure that is necessarily inaccurate and to use a computational method that is also inaccurate for simulating the waveﬁeld Green’s functions. We refer to these two sources of error as ‘propagation inaccuracies’, which might be considered to be epistemic. We show that the variance of the Fourier spectrum of the synthetic earthquake seismograms caused by propagation inaccuracies is related to the spatial covariance on the rupture surface of errors in the computed Green’s functions, which we estimate for the case of the 2009 L’Aquila, Italy, earthquake by comparing erroneous computed Green’s functions with observed L’Aquila aftershock seismograms (empirical Green’s functions). We further show that the variance of the synthetic seismograms caused by the rupture variability (aleatory uncertainty) is related to the spatial covariance on the rupture surface of aleatory variations in the rupture model, and we investigate the effect of correlated variations in Green’s function errors and variations in rupture models. Thus, we completely characterize the variability of synthetic earthquake seismograms induced by errors in propagation and variability in the rupture behaviour. We calculate the spectra of the variance of the ground motions of the L’Aquila main shock caused by propagation inaccuracies for two speciﬁc broad-band stations, the AQU and the FIAM stations. These variances are distressingly large, being comparable or in some cases exceeding the data amplitudes, suggesting that the best-ﬁtting L’Aquila rupture model signiﬁcantly overﬁts the data and might be seriously in error. If these computed variances are typical, the accuracy of many other rupture models for past earthquakes may need to be reconsidered. The results of this work might be useful in seismic hazard estimation because the variability of the computed ground motion, caused both by propagation inaccuracies and variations in the rupture model, can be computed directly, not requiring laborious consideration of multiple Earth structures.


INTRODUCTION
Numerical simulations of earthquake ground motions are used both to anticipate the effects of hypothetical earthquakes (e.g. Graves & Pitarka 2014) by forward simulation and to infer the behaviour of the real earthquake source ruptures by the inversion of recorded ground motions (e.g. Olson & Apsel 1982). In either application it is necessary to assume some Earth structure, represented as a seismic velocity, density and anelastic attenuation structure, which is necessarily inaccurate, and to use a computational method like ray theory or finite differences which is also inaccurate for simulating the wavefield Green's functions. Presently, there is no method for making quantitative estimates of the errors, caused by the errors in the assumed Earth structure and the inaccurate calculation method, in the simulated Green's functions for either type of study.
Typically in seismic hazard studies using synthetic seismograms (e.g. Wang & Jordan 2014;Villani & Abrahamson 2015)theerror in the synthetic earthquake seismograms caused by inaccuracies in the Earth structure and in the computational methods is not considered.
Recently, however, attention has turned to the inclusion of Green's function errors ('theory errors') in inversions of observed ground motions to infer the rupture process. In a very important paper Yagi & Fukahata (2011) presented a method to include the uncertainties in Green's functions into an inversion for earthquake rupture behaviour. Other significant works that look at theory errors in a Bayesian inversion context are papers by Bodin et al. (2012), Duputel et al. (2012), Dettmer et al. (2014), Minson et al. (2014) and Ragon et al. (2018). In all these papers the estimates of theory error were obtained from theoretical considerations. None of the investigators actually measured Green's function errors. Large earthquakes typically have aftershocks, which, if their rupture surfaces are physically small enough, can be considered to be point evaluations of the real Green's functions of the Earth. In our work we simulate small aftershock ground motions with (erroneous) theoretical Green's functions. The differences between the aftershock ground motions and the simulated motions are taken to be the 'theory error', and we develop a statistical model of the sources of discrepancies between the calculated and the real Green's functions. We use this with a frequency-domain version of the time-domain Yagi & Fukahata (2011) theory to derive the expected variance γ 2 caused by Green's function error in the Fourier amplitude spectra of ground motions from a larger (non-point) earthquake that we seek to model. This variance might possibly be useful both in ground motion inversions and in seismic hazard studies.
We also investigate the total variability of a synthetic seismogram for a large earthquake, which is caused by errors in the wave propagation model and by the natural variability of the earthquake source. We will show that the total variability of the ground motion is given by an equation with three terms. One term is controlled by errors in the Green's function, which are epistemic errors, one term is controlled by variations in the rupture source, which are aleatory variations, and a third term that is new to seismic hazard studies and quantifies the effects of correlated variations in the rupture models and the errors in the Green's functions.

A DISCRETIZED FREQUENCY-DOMAIN DERIVATION OF THE Yagi & Fukahata (2011)T H E O R Y , WITH ADDITIONS AND COMMENTS
In their insightful paper Yagi & Fukahata (2011) (henceforth YF) presented a time-domain theory in which they derived the effect of Green's function errors on the covariance matrix of predicted large earthquake ground motions and they included this covariance matrix in a ground motion inversion. However, their theory did not measure the errors in the Green's functions themselves. We derive a frequency-domain version of the YF theory, and we also add some terms which they have neglected, with comments.
Our notation varies somewhat from YF. Let the true traction on a fault at point x caused by a point force in the j-direction at the observer at y be g ( j) (x,ω; y) = [g ( j) 1 (x,ω; y) g ( j) 2 (x,ω; y)] T ,where 1-and 2-directions are orthogonal directions tangent to the fault at the point x , let its numerical approximation beg ( j) (x,ω; y), and the error ing ( j) be δg ( j) (x,ω; y), with expectation E[δg ( j) (x,ω; y)] = 0, so we have g ( j) (x,ω; y) =g ( j) (x,ω; y) + δg ( j) (x,ω; y) . (1) Similarly, let the relation between the true slip velocity distribution at x, the assumed slip velocity distribution and the variation in slip velocity be with expectation E[δs(x,ω)] = 0. We also have that the observed ground velocity in the j-direction d j (ω, y) equals the noise-free velocity v j (ω, y) plus the ground noise n j (ω, y), The noise-free ground velocity is given by Here, A is the rupture surface, the dot operator is the usual vector dot product of two vectors, and we have suppressed all the arguments for clarity. Eq. (5) is very important: it shows how variations in the rupture model and errors in the Green's functions contribute to the total motion. Discussing the terms individually, we definẽ Eq. (6) is just the usual finite-fault synthetic seismogram calculated using an approximate numerical method for an assumed rupture model in an erroneous assumed geologic structure. The ground velocity δv s j is caused by the variation δs in the rupture source. The variance ρ 2 j of δv s j is a measure of the aleatory variability in the ground velocity caused by the rupture source variability. The ground velocity δv g j is caused by the error δg in the assumed Green's function. The variance γ 2 j of δv g j is a measure of the epistemic variability in the ground velocity because errors in the geologic structure and the computational method can both be reduced with the addition of more information. Finally, δv sg j is the ground velocity caused by the interaction of δs and δg, and we denote its variance ξ 2 j .T h i s term has not been recognized in the seismic hazard literature and might not be negligible depending on the amplitudes of δs and δg. This term takes into account the possibility that variations in the geologic structure might cause correlated variations in the rupture behaviour, although this term will be nonzero even if δs and δg are uncorrelated. We will comment later on the meanings ofs(x,ω)and δs(x,ω).
Most of this paper will concentrate on estimating the variance γ 2 j of δv g j given real-world measurements of Green's function errors. However, we will show that the variance ρ 2 j of δv s j is given by a nearly identical result, and we will give an expression for the variance ξ 2 j of δv  Table 1). Red triangles indicate the broadband stations used in this work. We note that our γ 2 is identical to YF's τ 2 . We have changed notation to avoid confusion with the measure of interevent (betweenevent) variability τ used widely in the seismic hazard literature (Al-Atik et al. 2010).
YF do not explicitly write eq. (5). They make an important assumption. They assume that when their inversion converges, it converges to the true slip velocity model, so δs(x,ω) = 0. This enables them to drop the terms in δs above, yielding: We can de-clutter this equation by letting the index j refer to the jth channel of data, j = 1, 2, ··· , J where a channel is meant to be a single component of motion at a single observation location y and there are J total channels. Then observer location y is subsumed into index j and we have YF assemble their slip model from K constant-slip-velocity subfaults: where a qk is the total (final) slip of the kth subfault in the orthogonal q = 1andq = 2 directions on the fault, T (ω) is the Fourier spectrum of the slip-velocity time function normalized to unit total slip, assumed constant over each subfault, t k is the rupture time of the kth subfault, assumed constant over the subfault, and X k (x):= 1if x is in the kth subfault, := 0o t h e r w i s e .I ne q .( 12) we omit YF's sum over source-time functions without loss of generality. Let the integral of a single component q of the approximate Green's functions over the kth subfault bẽ with a similar equation for the integral of δg ( j) q (note -this latter integral is never evaluated). Combining these last two relations into   Table 1. the preceding one gives Further, assume that the data are filtered with some bandpass filter B(ω). Then the filtered data are Combining all these relations gives YF's eq. (7) in different notation: Here the term A qk (ω)isYF' stermP qk (t) * B(t). The first term on the right-hand side of the equal sign is just the slip-velocity model 'convolved' with the erroneous numerical Green's function, in other words it is the usual finite-fault forward synthetic in this type of inversion. Another way to write eq. (16)is where the double sum over slip-direction index q and subfault index k has been replaced by a single sum over p, with P = 2K and where Denoting the first term on the right-hand side of eq. (17)asu j (ω) and using the multidimensional delta method, the variance of u j (ω) for a single frequency ω is (Spudich et al. 2017; see the Appendix) where the dagger connotes complex-conjugate transpose, where we suppress the ω argument to avoid clutter and where We have used the transpose notation for typographical convenience to denote a column vector. Note that the right-hand side of eq. (20) contains only source terms; the data channel dependence j drops out. The pr element of the P × P complex covariance matrix C j , again with the frequency argument suppressed, is where r = 1, 2, ··· , P is a dummy index, where the asterisk denotes complex conjugation and where the expectation is taken over many realizations ofG. The covariance matrix C j has the dimension P × P, which is related to the gridding of the source on the rupture surface. In this frequency-domain derivation C j describes the covariance of Green's function errors on the rupture surface; it Downloaded from https://academic.oup.com/gji/article-abstract/219/1/346/5522610 by INGV user on 30 July 2019 Figure 5. More data-synthetic comparisons at AQU for RF structure. Observed aftershock ground velocity (red) and synthetic velocity (blue), plotted at the same scale as the observations. Numbers are peak velocities of data seismogram in cm s −1 . First 40 s of total 60 s seismograms are shown. Code AFnn is the aftershock number shown in Fig. 3(b). Data seismograms omitted from analysis owing to excessive noise or processing glitches are haloed in grey.
is not a data covariance matrix. We then have where the ω dependence is implicit. Because the covariance matrix is Hermitian, γ 2 j (ω) is real. To show that our derivation is correct, we reproduce the form of YF's result (their eq. 16, omitting the ground noise contribution) in the frequency domain by following them and letting where δ pr is the Kroeneker delta, σ 2 g (ω) is a scaling factor and S 2 pj (ω) is the maximum amplitude of the theoretical Green's function for slip component p and data channel j.Eq.(23) makes the reasonable guess that the covariance matrix is diagonal and its elements are proportional to the maximum amplitude of the theoretical Green's function. Inserting eqs (20)and(23) into eq. (19) gives where again we suppress the ω argument to avoid clutter. Note that eq. (24) lacks YF's Toeplitz matrix P qkj which contains the timeshifted source functions and is used in the time domain convolutions not needed in the frequency domain. While it is reasonable to guess that the Green's function error is proportional to the maximum value of the theoretical Green's function, it is preferable to measure the actual error, which we do later in this paper.

THE CONTINUOUS INTEGRAL CASE
YF pose their result in terms of constant-slip subfaults. The generalization to a continuous spatial variation of the slip function is straightforward. From eq. (11)wehave where the 2 × 2 covariance matrix of the Green's function error for the jth data channel is This result is obtained because the mean of the Green's function errors is assumed to be zero (the theoretical Green's functions are assumed to not have long-wavelength or frequency-domain biases), and inspection of the actual data show them to have a mean that is essentially zero, as we will see later. It should be noted that in the discrete case C j (ω) (with the lowered index) was a P × P matrix, where P was the number of grid points on the rupture and the covariance was taken between pairs of grid points, and C j (x, x ′ ,ω) (with a raised index) in the continuous case is a 2 × 2m a t r i x , where the covariance between points x and x ′ on the rupture will be given below by a continuous function. The use of a continuous function is one of the ways in which our work goes beyond YF. The raised and lowered indices are not meant to indicate covariant or contravariant components. Then, for a particular rupture model    s(x,ω), the variance of the ground motion caused by errors in the Green's functions is Physically we might expect two different functional forms for the spatial covariance. If it is dominated by finite-frequency effects, we might expect that an element of the covariance of the Green's function errors might be of the form where β is the shear wave speed of the medium and f is some decreasing function like a Gaussian centred at the origin. We expect scaling to be related to the shear wave speed because the S wave is typically much stronger than the P wave in finite-source seismograms. In such a model the errors in the Green's functions would be correlated at progressively longer distances as the wavelengths of the shear wave increased. On the other hand, if errors in the Green's functions are related to unmodelled spatial variations in the rigidity along the fault surface, this function might have no frequency dependence and would have the same form as the covariance of rigidity along the fault. Expanding eq. (27), we have where . (30) In these equations the ω dependence of the right-hand side terms is implied. Eq. (27) is the general result, but we can simplify for the common case that the rake in an earthquake rupture is primarily unidirectional. We can simplify by choosing theê 1 unit vector to be directed along the dominant slip direction and assuming that the other component of slip is negligible, which is true for the L'Aquila earthquake. We further simplify by assuming that theê 1 component of Green's function error is uncorrelated with theê 2 component, so C j becomes diagonal and we finally have We have verified numerically that the estimates of γ 2 j (ω) obtained from eq. (31) agree very well with the variance of an ensemble of ground motions calculated from an ensemble of δg having a given covariance function. In considering the differences between true Green's traction function g ( j) (x,ω) and its theoretical approximationg ( j) (x,ω), we expect that the theoretical approximation is developed in some simplified Earth model that is on the whole unbiased, and that the real Earth structure is equal to the theoretical structure plus some random heterogeneities. If a main-shock fault zone has an extensive low velocity zone, we assume that this feature is included in the theoretical traction functionsg ( j) (x,ω). The position vector x is a vector in 3-space, meaning that we are not assuming that our main-shock rupture is located on a 2-D surface (although Figure 11. Blue plus signs are the imaginary part of the covariance data formed from rigidity-and moment-scaled data for all three components of motion at station AQU using the RF velocity mode for the 0.26672-0.30006 Hz frequency band. Vertical positions of red plus signs indicate the median value of the covariance data in each distance bin, horizontal position is the centre of each distance bin. The median of the imaginary part of the covariance data at small separation is approximately zero, as expected. the aftershocks we use define a fairly planar structure). We are examining the variation in the traction wavefield in the main-shock source volume.
Following Frankel & Clayton (1986), the covariance function may be derived from observed teleseismic traveltime and amplitude variations between g ( j) (x,ω)andg ( j) (x,ω) that will depend on the wavenumber spectrum of the heterogeneities. Frankel & Clayton (1986) studied these variations in spatially stationary random media with Gaussian and exponential spatial correlation functions and in a self-similar medium with equal variation of seismic velocity over a broad range of length scales. Because their assumed random media were spatially stationary, their covariance functions were functions of spatial separation only. It should be noted that even though these random media had different spatial correlation functions, the amplitude distribution of the velocity perturbations was Gaussian. They concluded that random media with self-similar velocity fluctuations with a correlation length of a = 10 km can explain both teleseismic traveltime anomalies and the presence of seismic coda at high frequencies. Each of the three types of media has its own correlation function, shown in Fig. 1. The exponential and Gaussian covariance functions are normalized to unit amplitude at zero separation; the Von Karman correlation function is proportional to the Bessel function K 0 (r /a), where r is separation and is arbitrarily normalized to K 0 (1/10) = 1 in Fig. 1 (K 0 (0) =∞). Even for the self-similar heterogeneity spectrum, some nonzero covariance is expected at 10 km separation. Given a particular heterogeneity spectrum, we expect that there will be some spatially correlated traveltime and amplitude variations in the Green's functions. If we consider the covariance of true Green's traction function g ( j) (x,ω) and its theoretical approximationg ( j) (x,ω), we expect that there will be some nonzero covariance between different frequencies to account for traveltime errors. Some evidence of this appears in the next sections.
An important distinction to note is that Frankel & Clayton (1986) specified the covariance of their random seismic velocity structures, and their variations in wave amplitude were the result of the random structures. We, on the other hand, are using observations of aftershock seismograms to look directly at the random variation of the traction wavefield without the intervening mechanism of a random velocity structure. Our model is that the wavefield δg ( j) (x,ω) is spatially stationary (meaning that spatial covariance is a function of separation |x − x ′ | between δg ( j) (x,ω)a n dδg ( j) (x ′ ,ω)), and has zero mean, averaged over either many realizations or averaged over many correlation lengths. Consequently, the results of Frankel & Clayton (1986) are not directly applicable to our problem. What is needed now are finite-difference studies of the covariance of wavefields in three-dimensionally varying media. It might be that such covariance functions might have fundamentally different characteristics, such as having an oscillatory behaviour with the covariance function, being positive at some separations and negative at other separations.
YF's eq. (16) and our result (31) have implications for an important practical problem, namely, the problem of estimating the variability in a deterministic forward synthesis of ground motions from a hypothetical earthquake. When performing such a synthesis for a site-specific prediction of ground motion, the result has an uncertainty caused by ignorance of the exact seismic velocity structure of the region and an uncertainty caused by the variability of the rupture source behaviour. The YF result and the continuous integral form eq. (31) are important results because they quantify the errors caused by inability to calculate accurate Green's functions, both because of inadequate computational methods (e.g. ray theory versus elastic finite differences) and poor knowledge of Earth's structure.
Both of these sources of error can easily be better quantified after an earthquake by addition of more information, for example, a more accurate computational method or a better Earth structure model, and are therefore epistemic. Moreover, regardless of our state of knowledge of Earth's structure, some informed estimate of Green's function error covariance function can be made using functions. The second part of estimating the variability of a deterministic ground motion synthesis is the rupture variability. The aleatory part of the error in ground motion prediction is caused by the variability of the characteristics of the earthquake source. It could be that inclusion of a statistical representation of the source such as that by Mai & Beroza (2002) could yield the total variability of a Green's-function-based prediction of ground motion including both aleatoric and epistemic uncertainties. We comment more on this later in the paper.

ESTIMATING THE COVARIANCE MATRIX OF GREEN'S FUNCTION ERRORS
In order to measure the errors in theoretical Green's functions, we are assuming that recordings of ground motions from small aftershocks are true point-dislocation Green's functions (after some scaling to be discussed below). The error in the theoretical Green's function is related to the difference between the aftershock seismograms and theoretical point source simulations of the aftershocks.

The 2009 L'Aquila, central Italy, main shock and aftershock sequence
We have selected for study the 2009 April 6 M w 6.1 L'Aquila, Italy, earthquake and its aftershocks. The main shock and many aftershocks were well recorded by permanent strong motion instruments and broad-band sensors deployed in the epicentral area. In addition, Cirella et al. ( , 2012 have derived rupture models for the main shock using the simulated annealing algorithm of Piatanesi et al. (2007), so considerable preliminary work was already done. A map of the area is shown in Fig. 2. In this paper we concentrate on aftershocks recorded at broad-band stations AQU and FIAM that are essentially collocated with the strong-motion stations that recorded the main shock. Note that FIAM was collocated with strong motion station FMG. We chose to work with a set of 41 aftershocks recorded at station AQU (Fig. 3) that had focal mechanisms within 30 • of the main-shock mechanism and had corner frequencies above 2.68 Hz. We will comment on the possible biasing effects of this 30 • mechanism variation later.  As our highest frequency of interest for the study was 0.5 Hz, these events were point sources to a good approximation. 33 of these aftershocks were recorded at station FIAM. Theoretical traction Green's functions have been calculated for them using the frequency-wavenumber integration method (Saikia 1994;implemented  station FMG and the receiver function (RF) model of Bianchi et al. (2010), which was used for station AQU. The theoretical tractions were calculated for the fault geometry indicated by each aftershock moment tensor solution (Table 1)  The momenttensor solutions were obtained by fitting lower frequency data at more distant stations. Thus, fits shown in these figures are not the result of inversions of the observed aftershock data but rather of forward modelling of these events. We removed from the analysis those data for which the observed data had obvious ground noise or processing glitches. They were events labelled AF31, AF33, AF37 and AF39 for AQU, leaving 37 events, and AF28, AF23, AF20, AF12, AF11, AF07 and AF04 for FIAM, leaving 27 events.

Recovering tractions from ground velocity
Our observations in this section are differences between the Fourier transforms of aftershock ground velocities and synthetic ground velocities calculated using an assumed seismic velocity structure and the aftershock's seismic moment. Specifically, let index i be the aftershock index, i = 1, 2, ···n a ,wheren a is the number of aftershocks. Let index j indicate the jth realization of a random variable (the error associated with the jth erroneous velocity structure). For each frequency and component of motion we form the complex difference is the observed aftershock datum, which is a product of rigidity at the aftershock depth, the aftershock's true moment M i and its true traction Green's function g i ,ands j i is the synthetic from the jth erroneous velocity structure for the ith aftershock. We also implicitly assume that the above equations apply separately to the real and imaginary components of the complex quantities, so that when we write v i we mean Re(v i )orIm(v i ). LetP j i be our jth incorrect estimate of M i /μ i and letg j i be our jth incorrect synthetic traction. The incorrect moments might come from different sources, for example, a moment-tensor solution using broad-band data, or a long period spectral level from a 2 Hz geophone, and the incorrect traction Green's functions might come from different velocity structures. Then our synthetic aftershock seismogram is and (the real or imaginary part of) our complex difference is Normalizing by seismic moment, rigidity and a scaling factor yields a quantity with the units of the traction Green's function, namely, the scaled complex difference (or equivalently the empirical traction error) where c = 10 6 m 2 km −2 is a constant scaling factor whose only effect is to remove many decimal places from the plot annotations. We can use our empirical traction errors to determine whether the RF or the CIA model is better at AQU. Plots of the empirical traction errors differences in the complex plane for station AQU using the CIA and RF velocity models are shown in Figs 8 and 9, respectively. They show that the empirical traction errors grow in magnitude as frequency increases, justifying the appropriateness of our frequency-domain approach, and they also show that the CIA model is considerably better than the RF model (smaller complex differences) at AQU. However, since Cirella et al. ( , 2012 used the RF model for station AQU, we continue to use that model for AQU. The coloured dots showing the empirical traction errors for a particular aftershock are strung together with a black line. For many aftershocks the empirical traction errors form expanding helices, corresponding to progressive phase shifts as a function of frequency caused by time mismatches between the synthetic and real seismograms. The progressive phase shifts are evidence of nonzero covariance of Green's function errors at differing frequencies. Notably, we have not introduced 'static corrections' into the theoretical Green's functions to remove these time mismatches because time mismatches are errors in the theoretical Green's functions, the effect of which we hope to quantify.

Covariance as a function of separation on the fault
To determine the variance γ 2 j (ω n ) we need the spatial covariance matrix from (31), C j 11 (x i , x k ,ω n ), where j is the channel number (a single component of motion at a single station), x i and x k are two points on the fault (in this work, aftershock locations), and the 11-subscript indicates the covariance of Green's function error in theê 1 direction (taken to be the dominant direction of the mainshock slip) at x i with the Green's function error in theê 1 direction at x k . Including the index on ω, C j 11 (x i , x k ,ω n ) has six indices. For notational simplicity we suppress some of these indices by defining new spatial covariance matrix where we omit the slip direction indices entirely as the aftershock rakes are chosen to be within 30 • of the dominant slip direction.˜ j (x i ,ω n ) is the rigidity-and moment-scaled complex difference (or empirical traction) introduced earlier. We have previously speculated that the covariance would be a function of separation of the aftershocks, (which we will refer to as a covariance datum) as a function of r ik are shown by blue plus signs in Fig. 10 In order to take the expected values we should further average the covariance data within distance bins. However, the perturbations of the covariance data values caused by moment errors are skewed (because a moment error of a factor of 2 doubles some data while only halving others). Numerical tests of the effect of moment errors show that the covariance function inferred from the median of the covariance data is largely unbiased while the mean of the covariance data is systematically biased high (see the Appendix). The same argument applies to the effect of a 30 degree mechanism error. For some take-off angles the ratio of the true radiation pattern to the assumed radiation pattern will be greater than unity, and for others the ratio will be less than unity, with a ratio of infinity being possible in principle. While we have not checked this specifically, we believe that use of the median ameliorates the effects of mechanism errors.
For both stations AQU and FIAM we divided the distance range into 10 bins, with each bin width adjusted to hold the same number (about 70) of covariance data. To estimate the covariance as a function of separation we used the median of the covariance data in each distance bin. See the caption of Fig. 10 for more information.
The covariance function at zero separation is expected to be pure real, meaning that the imaginary part of the covariance data should be zero, which is shown in Fig. 11. The median of the imaginary part of the covariance data for all distance ranges and tested stations was not significantly different from zero, so we have assumed that it is identically zero for all subsequent computations.
The median covariance function for data from all components of motion in 10 frequency bands at AQU is shown in Fig. 12.
Surprisingly, if we normalize all these curves to unit amplitude at zero separation, we see different behaviours with frequency for the AQU/RF covariance data and the FIAM/CIA data. For AQU there is a conspicuous variation of the coherence functions with frequency, seen as the high-frequency covariance functions lying below the dashed average and the low-frequency data lying above the dashed average in Fig. 13, vaguely consistent with the behaviour predicted in eq. (28). On the other hand, a weak frequency dependence in the opposite sense of AQU is seen in Fig. 14 for FIAM. As the source zone for the aftershocks at AQU and FIAM are the same, it is difficult to identify a physical mechanism that would produce the differing frequency behaviours. The dashed average covariance functions in Figs 13 and 14 are similar in shape to Frankel & Clayton's (1986) covariance functions for exponential and self-similar media shown in Fig. 1, although our covariance functions suggest a covariance distance less than their 10 km. This suggests that the way to approach the theory-error problem in ground motion inversions is to treat the problem as one of waves in random media.

A covariance model
Our work goes beyond YF because we develop a smooth covariance function to encapsulate all the behaviour noted in the raw observations (e.g. Fig. 10) and hopefully to inject some theoretical insight into the smooth covariance function. We have noted empirically that the covariance functions for the three components of motion at AQU differ systematically, with the vertical component having the lowest covariance and the strike-perpendicular component having the highest covariance at high frequencies, whereas the component differences are not so strong at FIAM (Fig. 15). As the frequency Figure 20. The top row shows observed L'Aquila main-shock ground velocities (red) at station AQU and the synthetic ground velocities (blue) calculated from the best-fitting slip model shown in Fig. 16. The bottom row shows the Fourier amplitude spectrum of γ (ω) in black and the spectra of the data residual (observed minus synthetic) shown in green. Columns show the three components of motion. The expected error (black curve) for the strike-parallel component is much greater than the spectrum of the data residual, implying that the data for this component of motion are overfit. The expected error (black curve) for the vertical component is lower in the 0.2-0.4 Hz band than the spectrum of the data residual, implying that the data for this component of motion are underfit. dependence of the FIAM covariance was weaker and contrary to that of AQU, we arbitrarily chose to treat it as frequency-independent. An in-depth study on this is beyond the main goal of this work.
We are developing a continuous empirical covariance function to insert into eq. (31) in place of the discretized term K j ik (ω n )i n eq. (37). At this point it is necessary to separate the dependence on station s from the dependence on component of motion c.W e parametrize the covariance as a continuous function of r and ω as where N s (r ) for FIAM is the frequency-independent dashed average correlation function of Fig. 14 normalized to unit amplitude at zero separation, and sc (ω), c = 1, 2, or 3 and s = AQU or FIAM, is the covariance at minimum separation of the c component of motion at station s, as a function of frequency. To calculate sc (ω) we applied the same analysis as in Fig. 10 to individual components of motion, this time in three frequency bands with bandwidth 0.1667 Hz consisting of 10 Fourier components each (Fig. 15).
The centre frequencies of the three bands were 0.0917, 0.2584 and 0.4251 Hz. The piecewise linear curves in Fig. 15 were used to interpolate the value of sc (ω) at frequencies between the frequencies at which sc (ω) was measured (the circles, triangles and crosses). For AQU, because its covariance functions showed a systematic variation with frequency, which we wanted to preserve, the individual coloured piecewise linear functions in Fig. 13 were used for N s (r ) (making it frequency-dependent). Unfortunately, we do not have a smooth theoretical model of the expected behaviour of these curves N s (r ), and each curve is less well defined than the dashed average curve, so we expect that the ultimate effect of using these frequencydependent curves will be to make the frequency dependence of the resulting γ 2 (ω) somewhat irregular.

Epistemic error of finite-source synthetic seismograms
With this parametrization we now have everything necessary to estimate the epistemic error of our finite-source synthetic seismograms. Eq. (31) for the variance becomes where r =|x − x ′ |. As a computed example, for the rupture model s 1 (x,ω) we use the (unpublished) best fitting model found by , shown in Fig. 16. Because r =|x − x ′ |,eq.(39)isa spatial convolution over x (or x ′ ), performed using a 2D Fourier transform, followed by a simple area integral over the other integration variable. Derived values of standard deviation γ sc (ω)( c a l c ulated as the square root of real γ 2 sc (ω)) are shown in Figs 17 and 18 Downloaded from https://academic.oup.com/gji/article-abstract/219/1/346/5522610 by INGV user on 30 July 2019 for AQU using the RF velocity structure and Fig. 19 for the station FIAM (which is collocated with the strong motion station FMG) using the CIA velocity structure. The spectrum of standard deviation γ seen in Fig. 17 is serrated. The alternately depressed and elevated spectral levels are caused by use of the individual coloured piecewise linear median covariance functions for the ten different frequency bands from Fig. 13. Each of them is the result of averaging over three adjacent frequencies.
It will be possible to smooth these curves intelligently when some smooth theoretical model of the covariance function behaviour is available. For comparison, Fig. 18 shows the result of using the average covariance function (dashed curve in Fig. 13) for all frequencies. The results in Figs 17 and 19 are therefore used in the next section to discuss the goodness of fit of the main-shock synthetics and data.

Implications for inversion
One of the problems in most inversions for rupture behaviour is that no estimate of confidence in the resulting rupture models can be given because the effects of unknown errors in the Green's functions ('theory errors') have not been assessed. In inversion for rupture behaviour, these theory errors are usually bigger than the errors in data caused by ground noise and instrumental limitations. Consequently, it is probably common for such inversions to overfit the data, that is, to obtain an excessively good fit to the data considering the non-negligible theory errors. This meaning differs from the typical usage in statistics where overfit implies that a model has been overparametrized and fits noise as well as signal.
In Figs 20 and 21we compare the actual difference between the main-shock synthetic and observed ground velocities for stations AQU and FMG calculated using Cirella's best-fitting rupture model (shown in Fig. 16) with γ (ω), the Fourier amplitude spectrum of the expected standard deviation of the synthetic caused by errors in the Green's functions. In Fig. 20 the fact that the error spectrum (black) exceeds the data residual spectrum (red) for the strike-parallel component of motion may be telling us that these data at AQU are overfit by the inversion. For the strike-normal component the two spectra are in good agreement, possibly telling us that these data are being appropriately fit by the inversion. The residual spectrum exceeds the error spectrum in the 0.2 to 0.4 Hz band for the vertical component, telling us that these data are possibly under-fit by the inversion. The situation is different for FMG (Fig. 21), where all components of motion appear to be overfit.
These error estimates are very disturbing. Most investigators would consider the waveform fit of the strike-parallel component at AQU and all components at FMG as being good fits, and they would consider the strike-normal fit at AQU to be poor. If these error estimates are correct and representative of typical errors in rupture inversions using numerical Green's functions, it implies that the slip models of the many rupture models derived from past earthquakes are much less well resolved than currently thought, and a major reassessment of these models is necessary, a conclusion whichBeresnev(2003) has previously advocated.

Total variability of synthetic seismograms
We have thus far concentrated on finding γ 2 j (ω), the variance of δv g j , the Fourier amplitude spectrum of the ground velocity caused by errors in the Green's functions (eq. 8). Recall that eq. (5)shows that there are two other sources of variability in earthquake ground velocity: eq. (7) for variance ρ 2 j of δv s j , the Fourier amplitude spectrum (FAS) of the ground velocity caused by perturbations of the rupture model and eq. (9) for variance ξ 2 j of δv sg j , the FAS of the ground velocity caused by joint perturbations of the rupture model and errors in the Green's functions. Considering first ρ 2 j ,thev ariance of ground velocity caused by statistical variation in the rupture model is easily found, because the exact symmetry between eqs (7) and (8) implies that the variance is where and the ω dependence of the δs i terms on the right-hand side of eq. (41) is suppressed to reduce clutter. We can simplify eq. (41)b y assuming that δs 2 is negligible (this assumption is not necessarily true in all cases), yielding This shows that the variance ρ 2 j (ω) of the Fourier amplitude spectrum of the ground velocity, calculated for a specific Green's function model, depends on the spatial covariance of perturbations of the rupture model. At this point it is helpful to discuss the meaning of δs in greater detail. We foresee the primary use of eq. (42) will be in seismic hazard studies to calculate the variability of synthetic seismograms given an ensemble of rupture models. These rupture models might be obtained by performing many dynamic rupture simulations. It might be that in the seismic hazard study the investigator wants to include the effect of unknown hypocentre location, and the investigator's ensemble contains rupture models s that have many different hypocentres. Thens would be the average of all the individual rupture models (even though they have differing hypocentres), and S(x, x ′ ,ω) would be obtained from the deviations δs from the average.
The approach here differs from the pseudo-dynamic approach used by Schmedes et al. (2010), Song & Somerville (2010), Song et al. (2014) and Lee & Song (2017). In these four papers the spatial correlation of rupture parameters, such as rupture time and peak slip speed, is determined and then used to generate a family of rupture parameter models sharing the determined correlation statistics. This family of parameter models can then be used to generate a family of ground motions, from which the statistics (like the variance) of the ground motions can be determined. In eqs (40)and (42) rupture parameters are not used. The δs terms are variations in the Fourier spectrum of slip velocity itself. Further, the variance of the ground motion is determined directly, skipping the step of calculating the ground motions of many rupture models.
The final contributor to the variability of the ground velocity is the term given by eq. (9), the interaction between δg and δs. If slip is only in theê 1 direction, the variance ξ 2 j (ω)ofδv sg j , the interaction term, is The derivation of this result is straightforward following the steps in https://stats.stackexchange.com//questions/226657/variance-of-

strike -parallel strike -normal vertical
Comparison of data-synthetic mis ts with for FMG / CIA . The top row shows observed L'Aquila main-shock ground velocities (red) at station FIAM and the synthetic ground velocities (blue) calculated from the best-fitting slip model shown in Fig. 16. The bottom row shows the Fourier amplitude spectrum of γ (ω) in black and the spectra of the data residual (observed minus synthetic) shown in green. Columns show the three components of motion. The expected error (black curve) for all components is much greater than the spectrum of the data residual, implying that the data for all components of motion are overfit.
the-integral-of-a-stochastic-proceess.Wehaveverifiedthecorrectness of eq. (43) by direct calculation of the variance of an ensemble of δv sg j calculated from an ensemble of functions ψ(x) = δs(x) δg(x) having a given covariance as a function of |x − x ′ |.This expression (43) is new to seismic hazard research. Because |δs| or |δg| might not be small, ξ 2 j (ω) might be significant, even if the covariance between |δs| and |δg| is zero. A note of caution: it is tempting to use the rule that the variance of sum of the individual variances δv s j , δv g j and δv sg j is the variance of the noise free ground velocity, but correlations between the individual variance might exist that invalidate the rule in this case.

DISCUSSION
There are several unanswered questions left from this research. First, what is the proper functional form for the covariance of Green's function errors? We have previously speculated that it might be something like C j (x, x ′ ,ω) ∝ f (|x − x ′ |ω/β), but we have seen ambiguous evidence of its inverse scaling with shear wavelength β/ω. It might depend only weakly on spatial variations of material properties because some components of stress are continuous across material discontinuities. We also need to determine whether the two components of Green's function error are generally uncorrelated with each other, as we assumed in simplifying eq. (29) to get eq. (31). These questions could be answered by numerical studies of Green's function error in heterogeneous media. Similar questions apply to the covariance of rupture variations δs.
The possibility of a nonzero covariance between |δs| and |δg| in the error interaction term opens an interesting line of research. We can imagine that spatial variations of rigidity in the fault zone might cause spatial variations of |δg|. These spatial variations of rigidity might also cause correlated variations in the rupture process |δs|. It would be very interesting to look for such correlations in numerical simulations of spontaneous rupture in heterogeneous media.
A practical difficulty for applying eq. (31)forγ 2 j , the variance of ground motions calculated for a particular rupture model, is that the user has to specify a Green's function error covariance function. We have speculated that such a covariance function in regions having little ambient seismicity might be derived from measurements of teleseismic traveltime and amplitude anomalies and coda-Q, but this speculation must be checked in regions where there are ample small earthquakes that can be used as empirical Green's functions, Downloaded from https://academic.oup.com/gji/article-abstract/219/1/346/5522610 by INGV user on 30 July 2019 and direct observation of |δg| and construction of the covariance function are possible.
If we can determine Green's function error covariance functions for several different tectonic regimes, can they be regionalized? One might imagine that the covariance function for stable continental interiors might have more long-distance correlations than those for active tectonic regimes. If regionalization is possible, it would allow the investigator to use for peninsular India a covariance function derived for stable North America, for example.
We can offer only a few hints on how to include traveltime errors (or equivalently, phase errors) in the analysis. As our analysis gives the Fourier amplitude spectrum of the variance, the phase spectrum of the error is unspecified. Hence, specific realizations of the error could be generated from assumed phase spectra, which might be random or which might be phased to produce pulses at desired arrival times.
With the definition of covariance in eq. (26) including a complex conjugate, the resulting variance is real. If the complex conjugation is omitted, the result leads to complex valued covariances (called pseudo-covariances) and complex variances (Mandic & Goh 2009), which might be used in the future to include the effects of traveltime errors. Another alternative is to attempt to characterize the Green's function error covariance matrix in the time domain using empirical Green's functions, in other words form real covariance matrix from observed empirical Green's function errors. This is the approach of Yagi & Fukahata (2011). The resulting difficulty is that for practical calculations it is necessary to have some smooth, physically motivated interpolating function, like ours (38).
Finally, are ground motion inversions as corrupted by Green's function errors as Figs 20 and 21 seem to imply? There is one decision that we made that might have led to excessively large γ , namely our decision not to set the covariances to zero at large separations (see the Appendix). However, we did not have a good theoretical reason for doing so. One might argue that our comparison of Green's functions to aftershock waveforms (Figs 4-7) is quite bad; for example, there is no hope of fitting the aftershock codas. However, the observed main-shock seismograms contain the codas of early rupturing parts of the rupture, so it is necessary to consider the codas in the aftershock modelling. It might be possible to create a timevariable weighting scheme to reduce the weighting of aftershock late codas, but that is out of the scope of this paper. One might also argue that even after our removal of noisy aftershocks there seem to be some aftershocks left in our data set that have noise before the first arrival. However, such apparent noise might be actually a difference in the theoretical and observed first arrival times, which introduce phase errors which contribute to the covariances. This is the reason that our 'data tornadoes' (Figs 8 and 9) grow with frequency. An objective assessment of errors in source models must also include error caused by pulse traveltime errors.
It should also be noted that no measurement of Green's function error is typically done for crustal earthquakes, so it might be that rather than being catastrophically bad, our comparison is one of the best achievable. Figure A1. Test covariance data for the case σ r = 10, σ M = 1 (no moment error). Solid black line shows the target covariance function with σ r = 10, green symbols are the approximating simulated covariance data from the Cholesky method. Red crosses are the mean values of the covariance data in 10 distance bins; red circles are the median values in the same bins. Medians and means are very similar for this level of moment error (none).
Downloaded from https://academic.oup.com/gji/article-abstract/219/1/346/5522610 by INGV user on 30 July 2019 Figure A2. Test covariance data for the cases σ r = 10, σ M = 1 (no moment error) and σ r = 10, σ M = 4. Green symbols are the simulated covariance data from the Cholesky method with moment errors included. Green symbols at ordinates greater than 2 are not shown. Red and blue crosses are the mean values of the covariance data in 10 distance bins for σ M = 1andσ M = 4 respectively; red and blue circles are the median values for σ M = 1andσ M = 4 respectively in the same distance bins. Blue crosses with arrows are mean value data that exceeded the maximum plotting range. Downloaded from https://academic.oup.com/gji/article-abstract/219/1/346/5522610 by INGV user on 30 July 2019 Figure A3. Test to determine the bias of the median covariance as a function of σ M , moment error factors 2 (red), 3 (green) and 4 (blue) for the case σ r = 10 km. Mean ratio and population standard deviation shown at the centres of 10 separation bins containing the same number of covariance data. Blue and green symbols are shifted left and right slightly to avoid overplotting. Downloaded from https://academic.oup.com/gji/article-abstract/219/1/346/5522610 by INGV user on 30 July 2019 Figure A4. Test to determine the bias of the median covariance as a function of σ M , moment error factors 2 (red), 3 (green) and 4 (blue) for the case σ r = 3km, which is roughly the decay distance of the data. Mean ratio and population standard deviation shown at the centres of 10 separation bins containing the same number of covariance data. Blue and green symbols are shifted left and right slightly to avoid overplotting. Medians for all values of moment error are unbiased for separations less than σ r . Downloaded from https://academic.oup.com/gji/article-abstract/219/1/346/5522610 by INGV user on 30 July 2019