On the point-source approximation of earthquake dynamics

The focus on the present study is on the point-source approximation of a seismic source. First, we compare the synthetic motions on the free surface resulting from different analytical evolutions of the seismic source (the Gabor signal (G), the Bouchon ramp (B), the Cotton and Campillo ramp (CC), the Yoffe function (Y) and the Liu and Archuleta function (LA)). Our numerical experiments indicate that the CC and the Y functions produce synthetics with larger oscillations and correspondingly they have a higher frequency content. Moreover, the CC and the Y functions tend to produce higher peaks in the ground velocity (roughly of a factor of two). We have also found that the falloff at high frequencies is quite different: it roughly follows ~−2 in the case of G and LA functions, it decays more faster than ~−2 for the B function, while it is slow than ~−1 for both the CC and the Y solutions. Then we perform a comparison of seismic waves resulting from 3-D extended ruptures (both supershear and subshear) obeying to different governing laws against those from a single pointsource having the same features. It is shown that the point-source models tend to overestimate the ground motions and that they completely miss the Mach fronts emerging from the supershear transition process. When we compare the extended fault solutions against a multiple point-sources model the agreement becomes more significant, although relevant discrepancies still persist. Our results confirm that, and more importantly quantify how, the point-source approximation is unable to adequately describe the radiation emitted during a real world earthquake, even in the most idealized case of planar fault with homogeneous properties and embedded in a homogeneous, perfectly elastic medium.


Scientific rationale
One of the major challenges of the modern days seismology consists in the simulation of the motion resulting from a natural earthquake.In turn, this requires the complete understanding of the physical and chemical behavior of the fault zone (namely, the forces acing on the discontinuity interface), the knowledge of the geometry and structure of the fault system, as well as that of the medium where the seismic waves propa-gate.Geological (field) observations can provide some clues about the seismogenic structures.On the other hand, tomographic models (and inverse models in general) can help us to constrain the properties of the medium surrounding the seismogenic source, together with GPS and InSAR inversions, which, in turn, can provide useful information regarding the accommodation of the resulting deformation on the ground.Experiments conducted at the laboratory-scale and accurate numerical forward models represent powerful tools in the attempt to clarify the traction evolution on the fault surfaces, with the pivotal and amenable goal to formulate a comprehensive and realistic governing model which accounts all the possible energy-dissipating mechanisms occurring during a faulting episode (readers can refer to the thorough discussion on this subject made by Bizzarri [2011b]).
There is no doubt that there are two main categories of numerical models to simulate the wave propagation, the first being based on the so-called point-source approximation and the second being focused on the simulation of the chemico-physical processes on extended faults.The first approach is inherently simple, in that it contains only a point-source (or a double-couple), which accounts for the geometry of the fault (in an average sense) and for the final scalar seismic moment of a given earthquake.On the other hand, within the framework of the second approach it is possible to consider all the geometrical complexities of the source, its spatial heterogeneities (in terms of rheological properties) and, more importantly, the details of the stress release (the latter being described by the fault constitutive equations), which is ultimately responsible of the seismic waves excitation.
There is no reason to overemphasize that the first methodology has a trivial implementation in numerical algorithms, it is not computationally expensive and, moreover, it does not suffers of the epistemic uncer-tainty related to the state of the stress of the seismic source, its rheological behavior and all the details of its geometry (segmentation, branching, bending, etc.).This is a very important conceptual simplification of the whole problem.It is also important to remark here that the natural extension of the single point-source approximation mentioned above consists in the superposition of multiple point-sources, applied at different times (see Section 5), as retrieved by inversion of recorded data (i.e., kinematic models).Both of these approaches (single or multiple point-sources), by construction, cannot describe the complexity of the ground motions resulting from the various processes that take place on the fault [e.g., Andrews 2002, Bizzarri and Cocco 2006, Brantut et al. 2010, Bizzarri 2011a, 2012b, 2012d, among others].Finally, we also mention here the so-called stochastic method [Hanks 1979, McGuire and Hanks 1980, Hanks and McGuire 1981, Boore 2003, Boore 2009], which are used to simulate the mean ground motion for a given earthquake at a specific station.This method is still based upon a point-source approximation and it basically assumes that the source spectra is described by a single corner-frequency model [e.g., Brune 1970].It does not include any phase effects and the near-and intermediate-fields term are lacking.Finally, another approach is to somehow combine stochastic and deterministic models of earthquake sources, where the physics of earthquakes is considered only at low frequencies and at high frequencies ground motion modeling relies on the stochastic models [e.g., Liu et al. 2006, Ameri et al. 2009, Frankel 2009, Graves and Pitarka 2010].This kind of approach is more accurate than pure stochastic models in the simulation of the directivity effects.
In the literature there is a plethora of codes dealing with the point-source approximation; far of being exhaustive we mention here Virieux [1986], Coutant [1989], Coutant et al. [1995], Yomogida and Etgen [1993], Graves [1996], Moczo et al. [2002], Festa and Nielsen [2003], Tromp et al. [2008].Some of these codes also incorporate additional complications, such as the anelastic attenuation, the topography, the structural complexity of the medium surrounding the point-source, etc.At the same time, there is a very long list of functions used to simulate the time evolution of a seismic source; the box-car, the triangular function (or multiple triangular functions), the Gaussian, the Dirac delta function, the Gabor signal, the Yoffe function, the Bouchon ramp, the power laws, etc. [Yoffe 1951, Bouchon 1981, Archuleta 1984, Cotton and Campillo 1994, Hartzell et al. 1996, Ji et al. 2002, Liu and Archuleta 2004, Emolo and Zollo 2005, among many others].Some of these functions are not ballistic, other are singular.The most important point is that the choice of such a temporal evolution is not obvious, neither in the case of wave propagation problems (where the point-source approximation is used to simulate the resulting strong motions), nor in the case of kinematic rupture models (where the time evolution of the fault slip velocity is prior-assumed and used as a fault boundary condition to retrieve the stress drop and other friction parameters; e.g., Ide and Takeo [1997], Tinti et al. [2005]).Remarkabily, the point-source approximation is routinely used in the context of the CMT solutions.Piatanesi et al. [2004] compared three different functions to explore the effects on the stress change on the fault resulting from kinematic models, by specializing to the case of the rupture times and the final slip distribution of the 2000 Western Tottori, Japan, earthquake.Cirella et al. [2006] compared three functions and found that the peak fault slip velocity and the rise time retrieved in finite-fault inversions depend on the adopted analytical formulation.Bizzarri [2012c] compares three different analytical functions in order to see whether it is possible to reproduce the fault slip velocity resulting from a spontaneous (forward) model of an earthquake with closed-form analytical formulations.
The aim of the present study is twofold.First, we systematically compare the effects on the simulated ground motions of the assumptions of different pointsources approximations, in order to highlight their prominent features and limitations.One possible criterion to select the most desirable formulation can be the high frequency spectral decay of the resulting ground motion.Indeed, the famous ~-square model [Aki 1967, Brune 1970] -introduced theoretically for the very special case of constant stress drop modelhas been extensively used in source models [e.g.Boore 2003Boore , 2009, among many others] , among many others] and it has been also observed from some real data.However, it is important to emphasize that this spectral decay in not the unique; Boore ([2003]; see his Table 2 and references cited therein) reports a significant number of source shapes used in literature.Therefore, it emerges that there is no a universally accepted, theoretical or empirical reason to a priori select one specific formulation of the point-source time evolution; in this light the quantitative results presented in this paper can provide some general guidance.
Then we compare the results from a point-source approximation with those obtained from the modeling of spontaneous ruptures (both supershear and subshear) on an extended fault, governed by different constitutive laws, either the linear slip-weakening law or a rate-and state-dependent friction law.This comparison illuminates about the reliability and the possible limitations of synthetic seismograms computed by adopting of a single point-source instead of an extended fault model.This discussion is then complemented by the comparison between an extended fault model and a multiple point-sources model.
Both of these goals are still missing in the existing literature although they are very important.Indeed, it has not been yet clarified what are the limitations of the use of a specific formulation of a point-source time evolution or, at a more fundamental level, what are the limitations of the point-source approach per se.It should also be emphasized that the visual comparison between synthetic data from point-source models and recorded signals of given earthquake, routinely presented in many papers, is not a fully satisfactory validation of the method employed and does not answer the questions above.

The point-source approximation in earthquake dynamics
In this section we briefly recall the formalism of the point-source in the framework of the elastodynamic representation of a discontinuity interface.Readers can refer to Aki and Richards [2002] for a complete treatment of the subject.
In the point-source approximation the fault surface S is regarded as a system of couples located at a given point and the moment tensor is simply expressed as where M 0 being is the scalar seismic moment (M 0 = G <u tot >A, G being is the rigidity of the medium in which S is embedded, <u tot > being the total average slip developed by the event and A being the total fractured area at the end of the rupture process (A ⊆ S)), o and n are unit vector normal to and lying on S, respectively, and In Equation (2) 〈u(t)〉 A is the average value of the seismic potency [e.g., Rice 1980].
Equation (1) -which physically represents the strength of the resulting (j,k)-couple at the given point -has been derived under the assumption of a isotropic medium and of a shear dislocation (i.e., when no opening or material interpenetration is allowed), so that it results: u • o = 0.
The time derivative of s has an unitary time integral and it can be somehow associated to the source time function (STF thereinafter) of kinematic models of faulting processes (namely, to the quantity of Equation (6) in Bizzarri [2011b]; see also Section 5).Indeed, when an extended fault is described by mean of multiple elementary sources, s will be a function depending explicitly on the rupture time t r (p) of each source (namely, t would be t − t r (p)) and on the spatial distribution of the total cumulative slip on S (u tot (p)).The numerical factor (oj nk + ok nj) in Equation (1) expresses the geometrical features of the rupture and it can be easily expressed through the dip, strike and rake angles.
The simulation of a point dislocation source in a finite difference numerical scheme can be performed in terms of body-forces in the equation of motion; in particular, we follow the approach by Frankel [1993], which is appropriate for the displacement formulation in a conventional grid scheme.Just for an example, the equivalent force acting at grid points located ±1 nodes in the x 1 -direction from the source location equals M 11 (t)/2Dx, where Dx denotes the grid spacing and M 11 (t) is given by Equation (1).We recall here that the implementation of the moment-tensor source can be done either by particle velocity [Frankel 1993, Yomogida and Etgen 1993, Graves 1996] or by stress [Virieux 1986, Coutant et al. 1995].
Finally, we also mention that the moment tensor implementation of the source allows not only the modeling of the pure shear slip, but also that of explosive sources (equal diagonal elements and vanishing nondiagonal elements of the moment tensor) and that of a compensated linear vector dipole source.

Different temporal evolutions for the point-source
In the present section we introduce the most widely employed equations for the point-source description.Some of them has been used in kinematic inversions of recorded data, but they have been also extensively used to perform forward waveform simulations.Just for example, Moczo et al. [2002] use the Gabor signal (see Section 2.1) to simulate the wave propagation at high frequencies and to test different numerical methods.One of the most famous discrete wavenumbers code, AXITRA [Coutant 1989], used to simulate the wave propagation, originally employed the Bouchon ramp (see Section 2.2), although user can provide its own source time function.The modifications of the COMP-SYN method proposed by Spudich and Xu [2003] performed by Piatanesi et al. [2004]  (2) point-source approximation is adopted to simulate ground motions at high frequencies (far of being exhaustive, Luo et al. [2013], Lee et al. [2014], Magnoni et al. [2014], and references cited therein).

Gabor signal (G)
The function takes the form (see also Moczo et al. [2002]): and it is described by 3 input parameters f G , t G and c G .f G is the predominant frequency, c G represents the source decay rate, since it controls the width of the signal, and t G is merely a time shift.If reproduced symmetrically with respect to t = t G this function can be regarded as a smoothed triangular function (see also Komatitsch and Tromp [2002]).By changing the parameters of this signal the modeler can easily change the resulting signal from delta-like (i.e., with a broad spectrum) to monochromatic (i.e., with a narrow spectrum), still having the same analytical expression.A simplified version of (3) (i.e., without the cosinus factor) was used in Festa and Nielsen [2003] to test PML absorbing boundaries.In that case the time derivative will represent a Ricker function in slip velocity.

Bouchon ramp (B)
The so-called Bouchon ramp is written as it follows: and contains only one free parameter (t B ). Equation ( 4) represents a correction made by Cotton and Campillo [1994] to the function originally proposed by Bouchon ([1981]; cf. also the time integration of Equation (2a) of Piatanesi et al. [2004]), which was not ballistic.

Cotton and Campillo ramp (CC)
It is basically due to Cotton and Campillo [1995] and it is written as: (5) and it is described by a single parameter t CC .Equation ( 5) recalls very closely the Brune's model, in which the slip rate is expressed as it follows [Brune 1970]: in which m S is the S wave of the medium, G is its rigidity, Dx b is the breakdown stress drop (Equation (5) in Bizzarri [2011b]) and t c is the slip duration.In the Brune's model t c is controlled by the propagation speed and by the dimension of the rupture.Some possible modifications of Equation ( 5) can be also found in Beeler [2006].Moreover, the analytical form of has been also used by Imperatori and May [2013] to describe the moment release function.

Yoffe function (Y)
It descends from the time integration of the Yoffe function in slip rate [Freund 1979] and it can expressed as:

Liu and Archuleta function (LA)
In their nonlinear inversion model for the 1989 Loma Prieta, CA, earthquake Liu and Archuleta [2004] propose the following source time function: where C is a normalization constant and p ∈ [1,4].The special case p = 1 makes the spectrum of (7) practically identical to the Aki-Brune ~-square model [Aki 1967, Brune 1970].In this paper we consider the case of p = 4, which makes Equation ( 7) very similar to the function analytically found by Bizzarri [2012a], which represents the solution of the elastodynamic equation for a 1-D fault obeying the slip-weakening law.Intermediate values of p give functions similar to the other formulations presented in Section 2. When p = 4, we can write: In this case it exhibits a rapid decrease after its peaks, and this could cause significant radiation of seismic waves also at the healing front.

Synoptic comparison
A compendious synopsis of the functions described in Section 2 is reported in Table 1.From Equation (1) it descends that, since s(t) must be positive for all times, the analytical expression of a point-source has to satisfy the following condition: Moreover, since its time derivative represents the L 2 -norm of the vector slip velocity on the fault surface we also have the requirement: (10) In Figure 1 we compare the various formulations, for a typical set of parameters which guarantees that both the conditions ( 9) and ( 10) are satisfied.We can clearly see from Figure 1a that, in spite of their differences in a strict analytical sense, Equations (3), ( 4) and ( 8) have a very similar behavior, especially in the early stage of the rupture and in the increasing part; they differ near the unit slip, in that the (shifted) Bouchon ramp predicts a very gentle roll off and an asymptotic trend toward s = 1.On the other hand, Equations ( 5) and ( 6) also predict an overall similar behavior, but it is clear that in these cases the increase is abrupt (we emphasize that these functions are not of class C 1 for t = 0; see Table 1) and very fast, compared to that predicted by Equations ( 3), ( 4) and ( 8).
The same conclusions arise from the comparison of the time derivative of s(t) (Figure 1b).In particular, we can see the well-known singularity of the Yoffe function in slip velocity at t = 0 and that = 2/t CC (while the functions (3), ( 4) and ( 8) equal zero at t = 0).The peaks in Figure 2b pertaining to the three functions (3), ( 4) and ( 8) are not exactly synchronized; this reflects into the fact that the change of the convexity of s(t) is not at the same instant for the three cases.

Results from different point-sources approximations
In the present section we quantitatively compare the various formulations described in Section 2, by solving numerically the equation of motion in the case of a seismic source located at 7.3 km of depth and representing a vertical strike slip fault (the geometry is the same as that used, e.g., in Bizzarri [2011b]; his Figure 14a).The finite difference code, based on a conventional grid scheme, is the described in Bizzarri and Cocco [2005] and the implementation of the pointsource is recalled in Section 2.1 above.The parameters adopted in the numerical experiments are listed in Table 2.We also assume: M 0 = 7.45 × 10 20 Nm and <u tot > = 11.63 m.These values are not casual, but they come from a numerical simulation of a spontaneous (i.e., without prior-assumed rupture speed) 3-D, bilateral rupture, which propagates on a strike slip fault and obeying the linear slip-weakening friction law [Ida 1972].The adopted parameters are the same as Model A of Bizzarri et al. [2010], tabulated in Table 2 for completeness.
The results of the numerical experiments are reported in Figures 2 and 3, where we plot the time histories of the displacement rate at a selected receiver on the free surface.From a first, synoptic comparison of Figures 2 and 3 it emerges that results pertaining to the CC and the Y functions (Equations ( 5) and ( 6), respectively; see Figure 3) exhibit larger oscillations compared to those resulting from the other tree slip velocity functions (G, B and LA functions, Equations ( 3), ( 4) and ( 8 The function is not necessarily positive, depending on the adopted parameters.Just for an example, with f G = 0.225 Hz, t G = 1.5 s and c G = 1.5 s we have negative slip and slip velocity at t = 0.

Figure 1. (a)
Comparison between the five point-source time functions assumed in the present study, normalized by <u tot > (so that namely it is reported the behavior of s(t), as defined in Equation ( 2)).The adopted parameters are: f G = 0.225 Hz, t G = 1.5 s, c G = 1 s (Equation ( 3)); t B = 0.6 s (Equation ( 4)); t CC = 0.8 s (Equation ( 5)); t Y = 1.5 s (Equation ( 6)); t LA = 1.4 s (Equation ( 8)).In the case of the Bouchon ramp we apply a time shift of −0.respectively; see Figure 2).In particular, the results for the Y function are extremely noisy, also after the first peak, in all the three components of V.Moreover, if one compares the absolute values of the history of V it can be appreciated that CC and Y functions (Figure 3) tend to maximize the values of peaks; indeed we roughly have doubled values, compared to the other three functions (G, B and LA functions; Figure 2).On the other hand, the solutions reported in Figure 2 are quite similar; the first arrivals are at the same time and also the subsequent peaks are very similar.We can note a higher level of fluctuations in the case of the G function, compared to the B and LA functions.
Since the seminal paper by Aki [1967], it is known that one essential ingredient is the spectrum of the ground motions, where the physics of the earthquake processes, and resulting wave propagation, is encapsulated.This ground motion spectrum is directly connected (through the representation theorem; e.g., Aki and Richards [2002]) to the spectrum of the source.Indeed, the spectral analysis of the different functions also provide us with some clues (see Figures 4 and 5).As expected from the previous discussion, the solutions pertaining to the simulations with the CC and Y functions globally have a higher frequency content, compared to the other three functions.We can also see that the falloff at high frequencies roughly follows ~−2 in the case of G and LA functions, while the B solutions decay more faster than ~−2 (see Figure 4).On the other hand, both the CC and the Y solutions decay more slowly than ~−1 at high frequencies (see Figure 5).The results discussed above are not peculiar of the selected free surface receiver, but they are confirmed also for other stations.In Figures 6 and 7 we compare the results pertaining to a receiver located at the same distance from the epicenter along the strike direction as in Figures 2 and 3, but it is now more close to the fault trace (now is at a distance of 10 km instead of 20 km, as for Figures 2 and 3).The same results are confirmed at larger distance from the epicenter, as reported in Figures 8 and 9 (now the station is at 40 km from the epicenter and 5 km from the fault trace).
From the inspection of Figures 4 to 9 it emerges that at low frequencies (i.e., below 1 Hz) the solutions are somehow comparable.On the other hand, at high frequencies (i.e., above 1Hz, the range of prominent interest also for engineers) the resulting ground motions are quite different.
Figure 10 depicts the spectral falloff of the analytical STF considered in the present study.Namely, we compute the Fourier amplitude spectrum (FAS henceforth) of the time derivatives of s(t) (which are plotted in Figure 1b).If a function does not fail exactly to 0 in the target time window (as in the case of the CC and Y functions), we taper it to 0 by using a cosinus-tapering [Bizzarri and Spudich 2008]

Comparison between extended fault and pointsource approximations
In the present section we compare the solutions pertaining to a finite fault to those related to a point-source.The goal is to highlight possible limitations in the latter approach for some specific examples.In principle, we could try to compare data from a given real earthquake against synthetic data from point-source and extended fault and quantify the relative differences.However, the reproduction of a real-world event requires a very accurate description of the rheological properties of the fault, of the surrounding medium (including also possible site effects and topology).Therefore we prefer to simulate realistic, hypothetical earthquakes through a finite source modeling and then compare the results against the corresponding point-source approximation.To perform the comparison the choice of the input parameters is obviously of pivotal importance.We proceed as it follows: we first perform the extended fault simulation and then we employ the values of the scalar seismic moment and of the average cumulative fault slip as input parameters for the point-source simulation.A little bit more problematic is the choice of the rise time; the adopted value adopted in the point-source numerical experiments represents an estimate of the duration obtained in the extended fault simulations.We are not aware that this parameter can influence the results, so that the exact ratio between the resulting ground motions discussed in the following should be interpreted as an estimate and it is not intended to be read so literally.Indeed, previous findings of Bizzarri [2012c] clearly demonstrated that it is very difficult to reproduce the time evolution of the fault slip velocity of a spontaneous dynamic model with a closed-form analytical expression, especially with regard to the rise time.
Although relatively less frequent, this kind of event is receiving a great interest from the recent literature (see for instance the discussions in Bizzarri et al. [2010] and Bizzarri and Das [2012]).We expect that the directivity effects due to the finiteness of the fault are intrinsically neglected in the point source approximation of such an event.The comparison of numerical results will quantify how severe are the differences.
In Figure 11 we compare the solutions from the spontaneous model (left panels) and from the pointsource where LA function is assumed.In the case of the extended fault model we use again the values mentioned in previous section M 0 = 7.45 × 10 20 Nm and <u tot > = 11.63 m, which pertain to the parameters tabulated in Table 2.These values are then imposed as constraints for the point-source numerical experiment.We chose the LA function because it is compatible with the analytical solution of the elastodynamic problem in 1-D [Bizzarri 2012a].The adopted parameters are those of Table 2. From the spontaneous rupture model of the 3-D rupture we can clearly recognize the presence of the first peak of particle velocity, associated with the dilatational motion (the P 1 arc in Figure 11a,c,e) and the shear Mach (M S ) and Rayleigh Mach (M R ) fronts, originating because the rupture is supershear (compare panels (a), (c) and (e) with panels (a), (c) and (b), respectively of Figure 3 of Bizzarri et al. [2010]; note that the scale bars are slightly different and that the extension of the computational domain is different).These fronts are completely missed in the case of the point-source (see the right panels in Figure 11), indicating that a single point-source approximation of a seismic source is unable to account for the radiation emitted by a supershear earthquake, although the geometry, the seismic moment and the total slip of the event are the same (recall that we impose these values in the point-source simulation as obtained from the extended fault simulation).
Since the wavefields are quite different in the two models it is hard to make a direct, strict comparison of the absolute amplitudes of ground motions.Nevertheless, if we focus on the x 2 -component (Figure 11c,d) we can see the contribution due to the subshear front (the two spots close to the fault trace and located roughly at x 1 = 70.4km and x 1 = 67.5 km, respectively; see Figure 11c).In these regions, the ground motions pertaining to the finite fault are smaller than those pertaining to the point-source.This differences are even more pronounced if we consider other stations, especially in the direction x 2 perpendicular to the fault trace, where the point-source exhibits strong positive and negative pulses (see Figure 11b), contrarily to the finite fault (see Figure 11a).These differences can be related to different things.First, we must consider that the finite fault produces, due to the assumed governing model and homogeneous properties, a crack-like solution, while the LA function (as also the other STFs considered in the present study) has a pulse-like behavior.In this case, as mentioned above, the selection of the rise time is somehow problematic.
Another important issue is the role of the directivity.It is well known that all sites everywhere are af-fected by directivity; some stations are affected by ground motion amplification, other by de-amplification and other are neutral.An example is reported in Figure 1.5 of Spudich et al. [2013].Just for an example, by looking at the Shahai and Baker model (Figure 1.5b of Spudich et al. [2013]), which explicitly considers the high amplitude directivity pulses, it emerges that the receiver of Figures 2 and 3 of the present paper fails in the deamplification region, that of Figures 6 and 7 e) and (f ) Vertical component (V 3 ).In the left panels P 1 , M R , M S denote the first peak of particle velocity, the Rayleigh Mach and the shear Mach fronts, respectively.
gion.We will consider all these stations in the comparison reported in Section 4.2.

Case 2: Subshear earthquake
In this section we consider a rather different configuration.The fault is now governed by a rate-and state-dependent friction law, namely, the Ruina-Dieterich model (Ruina [1983]; see Equation ( 35) in Bizzarri [2011b]).This model is conceptually different with respect to the slip-weakening model employed in Section 4.1; further details are given in Bizzarri [2011b].We also consider a buried fault, so that the interactions of the rupture front with the free surface is inhibited by the presence of a region 1 km wide, where the rupture can not penetrate.Indeed, it is known that the free surface can induce supershear propagation also for configuration for which the subshear propagation is theoretically expected from the values of the assumed constitutive parameters [e.g., Olsen et al. 1997, Bizzarri 2010].This would happen also for the present configuration in absence of the shallow, unbreakable region we include in the present configuration.As expected, the resulting synthetic earthquake does not exhibit supershear speeds and, moreover, it can be considered more close to a pulse-like rupture.Indeed, we have recalled above that that the linear slip-weakening model with homogeneous properties, as that considered in Section 4.1, fails in the a so-called crack-like type of propagation.On the other hand, the STFs adopted here are characterized by a finite rise time (see Section 2.6 and Table 1); the present extended fault model can be therefore better compared against an analytical STF.
The adopted parameters are still listed in Table 2.The resulting scalar seismic moment, again computed from Equation (20) of Bizzarri [2011b] is: M 0 = 2.65 × 10 19 Nm and the average fault slip is: <u tot > = 0.818 m.The comparison between the resulting ground motions, at t = 7 s from the nucleation, is reported in Figure 12.Now that the Mach fronts are obviously absent, we can even better appreciate that, in general, the pointsource solution (right panels in Figure 12) predicts larger motions compared to the finite fault (left panels in Figure 12).This is especially true for the stations aligned along the direction perpendicular to the fault, as already observed in the supershear case (see Figure 11).
A more detailed comparison between the two models is reported in Figure 13, where the time histories of V are compared in three different free surface stations.To make the comparison more understandable, we normalize the displacement rate of each time serie by its peak-to-peak amplitude.It is clear that the point-source model (thin lines) predicts more oscillations with respect to the finite fault model (thick lines).
Moreover, the wave packets pertaining to the finite fault solution appears to be more spread out.The ratio between the peak-to-peak amplitude from the pointsource simulation and that from the extended fault simulation is reported in each panels.By focusing on the main pulses in the time series it is clear that the point-source overestimate the ground velocity.This degree of overestimate depends on the location of the receivers and on the selected component.It is safe to affirm that the overestimate roughly is of a factor 5 to 10 for most stations (in terms of absolute amplitudes of peaks).

Multiple point-sources
In the comparisons presented in previous Sections 4.1 and 4.2 all the seismic waves come, by definition, from a single point in the point-source approximations.On the contrary, in the cases of extended faults they come from the different points which fail and release stress.This in turn causes that the radiations pattern effects, as well as the wave attenuation, vary from one fault node to another.This geometrical effect can somehow be responsible of the differences observed and discussed above.
In the framework of the stochastic methods, Boore [2009] suggest the introduction of a correction, the effective distance measure, to tackle this effect.In the framework of the point-source forward simulation, the natural extension is to consider multiple point-sources instead of a single one.We reiterate that, from a physical point of view, in this approach all the energy available for a given event is concentrated into a single point (namely the hypocenter), where the total seismic moment (M 0 ) and the average slip (〈u tot 〉) are assumed to describe the static part of the event and where the (prescribed) temporal evolution of fault slip (〈u(t)〉) is assumed to be a proxy of the dynamic evolution of the whole fault (see Equation ( 1)).The next step in the modeling is therefore to assume that the energy is not concentrated only in a single point, but it is distributed through multiple seismic sources (or sub-faults, using a common, although misleading, nomenclature).It is postulated that in all sub-faults (namely, in all fault nodes) the components of the slip velocity are simply expressed as it follows (cf.also Spudich and Frazer [1984]): (11) in which the subscript i = 1 and 3 denotes the components (along the strike direction and down depth, respectively) and the function is the time derivative of the function s of Equation (2).In Equation (11) t r (x 1 ,x 3 ) is the rupture time distribution and u tot (x 1 ,x 3 ) is the , , , , BIZZARRI magnitude of the total cumulated fault slip.Both these two quantities are known from the spontaneous 3-D model.The quantity t r is defined as the time when o first exceeds the threshold value of o l = 1 cm/s (accordingly to previous papers; e.g., Bizzarri and Das [2012]).Note that the argument of is the retarded time t − t r , which states that the assumed fault slip velocity exactly is zero in a fault node until that point fails in the forward model.(On the contrary, in the extended fault model, the slip velocity can have values smaller than o l and non zero.)Namely, Equation ( 11) can be regarded as a kinematic model, where the fault slip velocity (and thus the rupture speed) is prior-imposed (and not retrieved as a part of the solution as in the spontaneous dynamic models).
Although we can expect to have some rake rotation, especially when rupture speed is supershear [Bizzarri andCocco 2005, Bizzarri andDas 2012], we conservatively assume here that all the slip is aligned in the direction of the initial shear stress (i.e., o 3 = 0 in Equation ( 11)); this constraint has been used by Hartzell and Heaton [1983], Kikuchi and Fukao [1985], Ward and Barrientos [1986], and Ekström [1989], among many others.
Figure 14   the kinematic model exhibits a slight delay with respect to the spontaneous model; this is because in the latter t r represents the time origin of the imposed fault slip velocity history in the former model (see Equation ( 11)).Compared to the single point-source (Figure 11b,d,f ) the actual wavefields are more similar to those of the dynamic model (Figure 11a,c,e).Indeed, in this model we recover all the features of the extended fault model that were missed in the single point-source simulation (Figure 11,right column).In particular, we can envisage the arc representing the first arrival (P 1 ), which corresponds to the supershear rupture front (located at x 1 = 75 km on the fault trace and at the selected time).Then we recognize the subshear front, located roughly at x 1 = 67 km (this front is clearly visible as the negative peak in the off-fault (x 2 ) component; Figure 14b).In between these two fronts we have, again as in the dynamic model, a positive contribution close to the fault trace (compare Figures 11c and 14b).We also remark that the distinction between the Rayleigh Mach and the shear Mach fronts is not so clear as in the case of the extended rupture.
In Figure 15 we compare the time evolutions of the displacement rate in the selected receiver R (indicated in Figure 14).From these plots we can see that the general behavior is similar in terms of peaks.We can also clearly appreciate the slight delay of the kinematic model with respect to the spontaneous dynamic model, discussed above.In terms of absolute values, it is also evident that the kinematic model tends to underestimate (of a factor of 2) the amplitudes of the signal; see the first arrival (P 1 in Figures 11 and 14) and the Mach pulse.

Discussion and conclusions
In this paper we have considered the propagation of seismic waves due to a point-source embedded in a perfectly elastic, homogeneous and isotropic medium.The wave propagation from the source in the medium surrounding the fault is computed through a Finite Different approach [Bizzarri and Cocco 2005].Contrarily to extended fault models, such as 2-D or 3-D rupture problems, the present model physically assumes that all the dissipated energy and the seismic wave excitation come from a single point, in which is concentrated all the energy available for the simulated earthquake event (namely, its seismic moment).The geometrical characterization of the fault (assumed to be uniform and planar) is provided by mean of the components of the moment tensor, i.e., through the orientation of the double-couple (see Equation ( 1)).
From the wave propagation point-of-view it is possible generate propagating waves by using many differ-ent signals, even signals which are not physical from the earthquake point-of-view.To describe the complex histories observed during real-world earthquakes the geometrical disorder and inhomogeneities should be taken into account [e.g., Xu and Knopoff 1994] and there is no doubt that, from a purely theoretical point of view, the point-source approximation should be certainly be regarded only as a first-order approximation of a fully determinist ground motion generation based upon a comprehensive model of an extended fault structure, which accounts for the underlying physics of the seismic source.Indeed, Boore [2009] suggests that there are some indications that the point-source is adequate only POINT SOURCES AND EXTENDED FAULTS  for small earthquakes and at substantial distances from the source.
Analytical expressions for the time evolution of a seismic source are currently incorporated in the forward modeling of the wave propagation [e.g., Moczo et al. 2002, Tromp et al. 2008, among many others].Moreover, they are often imposed a priori in ground motion inversion algorithms; indeed, although some authors obtain the time evolution of the slip velocity on the fault from inversions of the far-field or strong-motion data [Das and Kostrov 1990, Das and Suhadolc 1996, Das et al. 1996, Sarao et al. 1998], others have to impose it as an a priori assumption [e.g., Hartzell et al. 1996, Cirella et al. 2008, among many others].Finally, the time evolution of a seismic source is imposed, again as an a priori (and arbitrary) assumption, in kinematic rupture models, where a source time function (Equation ( 11)) is imposed to infer the stress-change fields [e.g., Bouchon 1997, Ide and Takeo 1997, Day et al. 1998, Dalguer et al. 2002, Mikumo et al. 2003, Tinti et al. 2005, among others].
Recently, also forward dynamic models can be also used to infer source properties; at a first stage, thousands of dynamic spontaneous models (in which the source time function is a part of the solution) are performed to produce seismograms, then they are compared with the observations and finally a best fit, based on some Monte Carlo inversion method, extracts those models which fit data best [Olsen et al. 1997, Ruiz andMadariaga 2011].However, the single-window kinematic (inversion) models mentioned above remain a current practice.It is therefore extremely important to clarify what kind of source function (or whether) is adequate or appropriate to describe the seismic radiation of an earthquake event and what are its main features.Indeed, since Piatanesi et al. [2004] it is known that different source time functions lead to different estimates of the stress drop in kinematic rupture models and, more importantly, we know that kinematic models may not be consistent with a well posed-dynamic fracture model.
In addition to fully deterministic, extended fault modeling [e.g., Aagaard et al. 2008, Olsen et al. 2009, Bielak et al. 2010, among many others] and to the stochastic methods [e.g, Boore 2003Boore , 2009, and references cited therein], the point-source modeling of wave propagation based upon the analytical functions described in the present paper (see Section 1.2) is a very common practice in the Geophysical community (far of being exhaustive, we only mention here a few recent examples: Tromp et al. [2008], Moczo et al. [2002], Magnoni et al. [2014]).Therefore, we believe that a thorough look into the details of the above-mentioned analytical functions, their features and limitation is timely and important.
Through a comparison of different analytical expressions proposed in the literature and adopted in many different papers we found that the ground motions pertaining to the Cotton and Campillo ramp (CC; Equation ( 5)) and to the Yoffe function (Y; Equation ( 6)) exhibit larger oscillations compared to models where the other velocity functions (see Figures 2, 3, 6a, 7a, 8a and 9a) and correspondingly they exhibit a higher frequency content (see Figures 4,5,6b,7b,8b and 9b).This is related to the abrupt increase of the slip pre- dicted by these two STFs at t = 0.Moreover, the CC and the Y functions tend to produce higher peaks in particle velocity (V), roughly of a factor of two, compared with the other expressions.We have also found that the falloff of V at high frequencies is quite different in the various cases: it roughly follows ~−2 in the case of the Gabor signal (G; Equation ( 7)) and of the Liu and Archuleta function (LA; Equation ( 8)), it decays more faster than ~−2 for the Bouchon ramp (B; Equation ( 4)), while it decays more slowly than ~−1 for both the CC and the Y solutions.The differences in the spectral signature of the solutions on the free surface stations reflect those of the adopted STFs; from Figure 10 it is clear that CC and LA decay as ~−1 and ~−2 , respectively, Y has a more flat spectrum, while G and B decays faster than ~−2 .Since there are no theoretical or empirical reasons to a priori select a specific STF these results provide a quantitative guidance in the point source modeling.
Another interesting outcome of the present study is that we confirm that the single point-source approximation is totally unable to account for the radiation emitted by a supershear earthquake (see Figure 11; we have shown this for the LA function, but the results is the same for the others expressions).This is expected since the point-source does not predict the Mach cones originated by a supershear event [e.g., Dunham andBhat 2008, Bizzarri et al. 2010].This is a very relevant limitation if we consider that there is an increasing evidence that some of major faults have the tendency to generate supershear events (see for instance Das [2010] and references cited therein).
Moreover, we found that for most of the free surface receivers, the single point-source models tend to overestimate the ground motions; this propertywidely discussed in the engineering seismology community; e.g., Boore [2009] -is explaining why the observed near-field saturation of ground motions is magnitude-dependent.This results holds for both super-and subshear events and irrespective of the choice of the adopted governing model (see Figures 11  and 12).Remarkably, these differences persist also at relatively large distances from the fault (up to several tens of chilometers; see Figure 13).
The insertion of multiple point-sources (see Section 5) is able to recover some of the features of the supershear ruptures (see Figure 14) and conceptually solves the problem that in the single point-source approximation the wave attenuation and the radiation pattern effects do not vary from one fault point to another.However, if we compare the results against those pertaining to an extended fault we can see that there is no a clear distinction between the Rayleigh Mach and the shear Mach wave fronts, (compare Figure 11, left panels and Figure 14) and that the multiple point-sources model tends to underestimate the ground motions.
It is important to remark here that the differences envisaged between the two approaches (extended fault and point-source models) can be partially influenced by the fact that it is very difficult to fit exactly the shape of the fault slip velocity resulting from a spontaneous, 3-D, fully deterministic model with an analytical function, especially with respect to the choice of the rise time (this issue has been also discussed by Bizzarri [2012c]).
On the other hand, we also emphasize that the differences in the radiated wavefields resulting from a spontaneous dynamic model and from a single pointsource would become even more pronounced if we incorporate in the forward dynamic modelling of the extended fault additional energy-dissipating mechanisms which can occur during the stress release process (e.g., thermal pressurization of fluids, melting, decomposition weakening, microcracks healing etc.; see Bizzarri [2014]) and which are inherently missed in the simplistic point-source approximation.The same is true also if we consider that earthquakes (especially large earthquakes) occur on quite complex fault systems, which include bends, step-overs, branching, jogs [e.g., Sibson 1986, Sieh et al. 1993].These complexities are neglected by definition in the point-source approximation.
As an overall conclusion, the present work provides some caveats in the adoption of a STF.Indeed, we stress that the choice of the analytical evolution of the STF does really matter, in that it can significantly affect the resulting ground motions and their spectral signature.More importantly, we have shown (and quantified) that a point-source approximation is inherently unable to reproduce adequately well the radiation excited by the propagation of a rupture on a fault of finite width, even in the idealized condition of planar fault with homogeneous rheology, up to relative large distances from the hypocenter (several tens of km).In other words, the modeling of the synthetics at regional distances from the source, modeled though a single point-source approximation, is biased by an improper, oversimplified description of the source itself.This conclusion is reasonable and not surprising, but since the point-source approximation is still widely used, the present results can provide some general and critical guidance to the modeler and to the seismological community.
Figure1.(a)Comparison between the five point-source time functions assumed in the present study, normalized by <u tot > (so that namely it is reported the behavior of s(t), as defined in Equation (2)).The adopted parameters are: f G = 0.225 Hz, t G = 1.5 s, c G = 1 s (Equation (3)); t B = 0.6 s (Equation (4)); t CC = 0.8 s (Equation (5)); t Y = 1.5 s (Equation (6)); t LA = 1.4 s (Equation (8)).In the case of the Bouchon ramp we apply a time shift of −0.7 s; the dashed curve represents the unshifted function.(b) The same as in panel (a), but now we plot the time derivative of s(t).

Figure 2
Figure 2 (left panels).Time evolution, low-passed at 8 Hz, of the displacement rate at a free surface receiver R located at a distance of 25 km from the epicenter and 20 km from the fault trace.The location of R and the assumed fault geometry are reported in Figure 14a of Bizzarri [2011b].Different colors define different slip velocity functions (Gabor signal, Equation (3), Bouchon ramp, Equation (4) and Liu and Archuleta function, Equation (8)).(a) Fault-parallel component (namely, V 1 ).(b) Fault-normal component (namely, V 2 ).(c) Vertical component (namely, V 3 ).Figure 3 (right panels).The same as Figure 2, but now we consider the Cotton and Campillo ramp (Equation (5)) and the Yoffe function (Equation (6)).
. The tapered signals are then normalized to regain the unitary slip predicted by the original functions.It is clear from Figure 10 that the various STFs have different spectral decays, which also confirm the differences observed in the FAS of the particle velocity.In particular, the CC and LA functions exhibit exactly an ~−1 and an ~−2 decay, respectively; the Y function decays more slowly than ~−1 and finally the G and B function decay faster

Figure 5 .
Figure 5.The same as Figure 4, but now for the time histories plotted in Figure 3.

Figure 8 .
Figure 8.(a) The same as Figure 2a, but now for a free surface receiver located at a distance of 40 km from the epicenter and 5 km from the fault trace.(b) Corresponding FAS.

Figure 7 .
Figure 7. (a) The same as Figure 3a, but now for a free surface receiver located at a distance of 25 km from the epicenter and 10 km from the fault trace.(b) Corresponding FAS.

Figure 9 .
Figure 9. (a) The same as Figure 3a, but now for a free surface receiver located at a distance of 40 km from the epicenter and 5 km from the fault trace.(b) Corresponding FAS.

Figure 10 .
Figure10.FAS of the time derivative of the STF adopted in the present study (namely the FAS of is reported; see Section 3 for further details).

Figure 11 .
Figure 11.Comparison between the results from an extended fault model (left column) and those pertaining to a single point-source approximation (right column).In the former case we assume that the spontaneous, supershear rupture, obeying the linear slip-weakening friction law.In the latter case we assume a LA function.The parameters are those of Table 2 and the solutions are plotted at time t = 7 s from the nucleation.Due to the symmetry of the wavefields we only report the distribution on a quadrant of the computational domain.The fault trace is marked by a thick dashed line.(a) and (b) Fault-parallel component (V 1 ).(a) and (d) Fault-normal component (V 2 ).(e) and (f ) Vertical component (V 3 ).In the left panels P 1 , M R , M S denote the first peak of particle velocity, the Rayleigh Mach and the shear Mach fronts, respectively.

Figure 12 .
Figure 12.The same as in Figure 11, but now the finite fault model (left panels) pertains to a subshear rupture governed by the Ruina-Dieterich law (the parameters are listed in Table2).The point-source solution (right panels) has the same seismic moment and total slip (see Section 4.2).

Figure 13 .
Figure13.Comparison between the time histories of the particle velocity components from the extended fault model (thick lines) and the point-source model (thin lines).Values of V are normalized by the peak-to-peak amplitudes.The models are the same as in Figure12.The location of the selected free surface receivers is indicated on the top panels, for each column.Top raw pertains to the x 1 -component, middle raw to the x 2 -component and bottom raw to the x 3 -component.In each panel is reported the ratio between the peak-to-peak amplitudes, namely that of the point-source model divided by that of the extended fault model.

Figure 14 .
Figure 14.The same as in Figure 11, but now in the case of multiple point-sources (see Section 5 for details).

Figure 15 .
Figure 15.Comparison between the displacement rate at the free surface receiver R, located at a distance of 25 km from the epicenter and 20 km from the fault trace (R is indicated in Figure 14), pertaining to the extended fault model of Figure 11, left column, and to the multiple point-sources.All signal are low-passed at 8 Hz.(a) V 1 .(b) V 2 .(c) V 3 .

Table 1 .
Synoptic comparison of the most important features of the considered time evolutions of slip on the fault (〈u(t)〉) of the pointsource approximation.(*) In asymptotic sense.( † )

Table 2 .
Parameters adopted in the present study.(*) We also apply a time shift of −0.7 s.