Earthquake source parameters such as rupture size, duration, velocity, and propagation, help us understand earthquake physics, fault zone properties, and rupture dynamics conditions, which have significant implications for seismic hazard assessment. Effects due to source finiteness, such as directivity, are usually associated with high-magnitude earthquakes (Ammon et al., 1993; Somerville et al., 1996), but they can also enable moderate earthquakes to cause severe unexpected damage. For example, the directivity effect can lead to potentially destructive pulses at low frequencies characterized by large amplitudes of ground motion (Boatwright, 2007; Kurzon et al., 2014; Moratto et al., 2017; Ertuncay and Costa, 2021; Ertuncay et al., 2021). Therefore, knowledge of the kinematic finite source parameters and expected rupture directions for high and moderate magnitude events is critical for earthquake engineering applications (Moratto et al., 2021; Somala et al., 2021; Moratto et al., 2023) and for appropriate risk assessment.

Estimation of source parameters for high magnitude earthquakes is usually feasible. However, the accurate determination of these parameters for moderate and small earthquakes remains a major challenge. Namely, while far-field records and geodetic data can provide complementary information for large earthquakes, geodetic data are lacking for small to medium events, and only seismic records from nearby stations are available. As a result, the kinematic properties of small earthquakes are often difficult to determine, and simple models are often used to represent these events (e.g., McGuire, 2004; Shearer et al., 2006; Convertito et al., 2013; Calderoni et al., 2015; Abercrombie et al., 2017; Moratto et al., 2019; Colavitti et al., 2022; Yoshida et al., 2022) although improved records show that source complexity is also common for small earthquake ruptures (e.g., Calderoni and Abercrombie, 2023 and reference therein).

A critical task in determining finite source attributes for moderate and low magnitude earthquakes requires good removal of path and site effects. For large earthquakes, this is usually done by synthesizing Green’s functions with a suitable velocity model to fit the low-frequency component (less than 5 Hz) of the observed waveforms (e.g., Ji et al., 2002; Beresnev et al., 2003; Yue et al., 2012; Moratto et al., 2015). For small to moderate earthquakes, the existing velocity structures usually have insufficient resolution and accuracy, leading to reliability problems in the high-frequency range (greater than or equal to 5 Hz). To address this problem, a number of methods based on empirical Green’s function (EGF) deconvolution have been developed in recent decades (e.g., Hartzell, 1978; Mueller, 1985; Mori, 1993; Ammon et al., 1993; Hough, 1997; Lanza et al., 1999; McGuire, 2004; de Lorenzo et al., 2008; Meng et al., 2020). EGF methods use a collocated earthquake with magnitudes typically 1.5–2.5 units smaller than the target event with a similar focal mechanism to account for path effects. The deconvolution process of the target event by an EGF automatically removes path and site effects without the need for a velocity model to synthesize Green’s functions. Although the EGF offers several advantages, its use presents certain difficulties, and selecting an appropriate EGF can be difficult, even when working with an extensive database (Calderoni and Abercrombie, 2023). Focal mechanisms are often not available for small earthquakes, and numerous studies have shown that even for small earthquakes, the effect of directivity cannot be neglected when selecting an EGF (Abercrombie, 2015; Calderoni et al., 2015; Calderoni et al., 2017).

The simplest general representation of an earthquake that contains information about rupture extension and directivity is the point-source representation plus the variances or second-degree moments of the moment-release distribution (e.g., Silver, 1983; McGuire et al., 2001). The hypocenter and origin time of the earthquake correspond to the spatial and temporal average (first-degree moment) of the release moment distribution. The information on rupture extension, characteristic duration, and direction of rupture propagation correspond to the variance of the moment distribution in the spatial, temporal, and spatiotemporal domains (second-degree moments). Seismic moments are calculated from apparent durations measured from apparent source time functions (ASTF) for each station after removing path effects. The ASTF obtained using EGF deconvolution (McGuire, 2004; 2017), represents the difference between the source time function (STF) of the main earthquake and that of the EGF earthquake from the direction of observation at each station. Thus, the ASTF is the projection of the rupture process onto the seismic ray path, and its properties also depend on azimuth and take-off angles (e.g., McGuire et al., 2001; McGuire, 2004; Stich et al., 2005; Yoshida and Kanamori, 2023). For a unilateral rupture, the ASTF observed by stations in the propagation direction would be significantly shorter than the ASTF of stations in the opposite direction.

A major advantage of the second moments method is that it can be theoretically applied to all earthquakes, regardless of their magnitude and complexity, and without requiring the assumptions of an a priori source model (e.g., McGuire, 2004; Fan and McGuire, 2018; Meng et al., 2020; Meng and Fan, 2021). It is also a consistent tool for evaluating scaling relationships between finite source attributes and earthquake magnitudes for large and small earthquakes (McGuire and Kaneko, 2018) and for resolving fault-plane ambiguity. However, elimination of the path effect is critical, and a biased ASTF calculation would lead to inaccurate calculations of the second seismic moments.

In this study, we consider second-degree seismic moments to provide a simple, model-free description of earthquake rupture and we evaluate, through synthetic tests, the sensitivity of the solutions to uncertainties in key input parameters to investigate the limitations of the method in the study of small earthquakes.

To achieve this goal, we conducted a synthetic test for a magnitude 4.6 earthquake and examined the effects of potential uncertainties in source duration, hypocentral depth, station configuration, focal mechanism, and velocity model, considering unilateral and bilateral rupture scenarios. In view of future analyses to be performed with real data, we located the hypothetical earthquake in central Italy, an area with good station coverage by the Italian seismic network (Amato et al., 2006) and the Italian accelerogram network (Costa et al., 2022). The wealth of data recorded in this very active seismic area makes it an excellent natural laboratory for the study of source characteristics, and there are numerous studies on the subject of sources published after the 1997 Umbria and Marche earthquakes, 2009 L’Aquila, and 2016 Amatrice earthquakes (e.g., Chimera et al., 2003; Bindi et al., 2009; Cultrera et al., 2009; Moratto and Saraò, 2012; Rovelli and Calderoni, 2014; Calderoni et al., 2015; Calderoni et al., 2017; Convertito et al., 2017; Wang et al., 2019).

The following sections describe the method used to determine the source parameters and the synthetic input model, as well as the sensitivity tests performed. Finally, the results and the strengths and weaknesses of the method are discussed, providing valuable insights for applying the approach to real earthquake data.

The input parameters used to model the ASTF for a Mw 4.6 earthquake source are listed in Table 1. We assumed that the epicenter was located in central Italy (Figure 1), and approximated the fault as a 3.0 km box model (Figure 2). The rupture area was divided into 12 × 12 cells, and the slip distribution and rupture time for the unilateral (Figures 2A, B) and bilateral (Figures 2D, E) scenarios were taken from a previous study of an earthquake of similar magnitude (Lopez-Comino et al., 2016) downloaded from the SRCMOD database (Mai and Thinbgaijam, 2014), with a focal mechanism of strike 247°, dip 46°, and rake 40°. This fault plane is not representative of the tectonic regime that characterizes the Apennine chain in central Italy, but this issue is not crucial for the synthetic tests. A uniform propagation of the rupture front was assumed with a rupture velocity of 2.75 km/s, which corresponds to 0.9 times the S-wave velocity in the source region. A simplified 1-D velocity model of central Italy (Costa et al., 1993) was used to model the ASTF (Figures 2C, F). We refer to this model as A (Figure 3A). To check the robustness of the solutions, we used an additional model for the inversion (Figure 3B), slower on average than model A, which we referred to as model B (Costa et al., 1993).

To investigate how the uncertainties introduced by the input data may affect the solutions of the resolved second seismic moments, we used the bootstrap approach. In this technique, perturbations are introduced for each input parameter to be analyzed by generating 1,000 variations around the mean. An inversion is then performed to evaluate the effects on the mean and standard deviation of the resulting data. The workflow is summarized in Figure 4. We examined the uncertainties associated with the apparent source durations τ c s _ , the hypocenter location, the station distributions around the source, the focal mechanism, and the velocity model used for ray tracing. Some of these tests are interrelated. For instance, uncertainties in both the hypocenter location and the velocity model affect the calculated ray path, and both the different focal mechanism and station coverage affect the resolution of the fault plane.

The uncertainties in the epicenter estimates are not examined because they have negligible effects on the slowness vectors in the inversion of the second moments.

The use of second-moment tensors to determine the source parameters, including directivity, of moderate earthquakes could be a valuable tool to improve our understanding of source dynamics in a given area and to advance disaster mitigation. A potential application of the second-moment method to small earthquakes would be to identify portions of large faults that produce supershear ruptures and correlate them with the geology of the fault zone. The second-moment method also provides lower constraints on rupture velocity, which may be particularly useful for unilateral ruptures (McGuire and Kaneko, 2018). However, before the results can be interpreted, it is necessary to know the resolution limits of the method due to the possible uncertainties of the parameters used as inputs to the computational procedure. To this end, we performed sensitivity tests of the method using synthetic ASTFs.

The sensitivity analysis performed in this study shows that the uncertainties in the input data have different effects on the calculation of the source parameters and an accurate measurement of the ASTF as well as the velocity model play the most important role in influencing the inversion process. A visual overview of the results of our tests is provided by the violin plots (Figures 9, 10) for unilateral (Figure 9) and bilateral (Figure 10) rupture scenarios. This type of plot is particularly useful for comparing the distribution of data between multiple categories or groups, as it effectively displays medians, ranges, and variabilities and facilitates comparison between different tests.

We observe that the main source parameters, i.e., rupture size, source duration, and centroid velocity, are generally well reproduced within the standard deviation. As expected, the source duration resulting from the inversion process is strongly affected by the duration of the input ASTF, and even 10% affects the inversion of the second moment tensor. This is because the ASTF is strongly related to the second moments in time, which is a vector in the inversion of the linear system (Eq. 8).

In the case of dense instrumentation, the horizontal location of the earthquake can be well resolved, but the resolution of the earthquake depth is largely determined by the velocity model, and inaccurate earthquake location can lead to uncertainties in the resolved second moments. In our case, perturbation of the hypocentral depth leads to artifacts in most of the resulting variable distributions. This is primarily due to the fact that the perturbed depth in the 1-D velocity models spans multiple layers, in which the slowness does not change once layer boundaries are crossed. In the case of model A, there is an interface between layers with a velocity difference of 0.7 km/s. As a result, the ray paths vary considerably, which shows up in the calculations as two prominent peaks in each variable distribution. Using a velocity gradient between layers, most of the parameters studied are well represented by a one-peak distribution. Care must also be taken to avoid artifacts due to the discretization of the velocity model when the hypocenter is located at an interface between two layers with high velocity contrast. A perturbation of the velocity model A by 0.3 km/s does not significantly affect the results for the fault dimension ( L c and W c ). Previous studies performed by McGuire and Kaneko (2018) for various source models have shown that the fault area can be well estimated when at least 15 ASTF measurements are available at a variety of take-off angles with upward and downward rays. The exact relationship between L c , W c , and the total rupture dimension also depends on the particular slip distribution. McGuire and Kaneko (2018) have shown that data sets with about 30 or more measurements provide a very good recovery of the rupture velocity. We observe that the values of directivity depend on the ASTF duration, the choice of velocity model, and the focal mechanism (Figures 9D, 10D). To ensure good fault plane resolution, good ray path coverage for both upward and downward waves is critical (McGuire, 2004; McGuire and Kaneko, 2018), especially in poorly studied areas or complex geostructures. McGuire (2004) found that the horizontal component of the rupture directivities is generally well constrained and has good azimuthal coverage for strike—slip events. The component of rupture directivity along dip can only be well determined when stations directly above the hypocenter are available, since seismic rays at most other stations are nearly horizontal.

Previous research has already shown that the second-moment method can be implemented with sufficient arrays to resolve the rupture dynamics of small earthquakes (M ∼ 2.5) with high accuracy (McGuire, 2004; Fan and McGuire, 2018; Meng et al., 2020). However, in our analyses, the perturbation of the station configuration did not affect the calculation of the source parameters, due to the particular input model we considered and the accuracy of the synthetic test. Further testing using different source models and real data would be needed in the future to investigate this issue.

The choice of an inappropriate velocity model leads to a significant bias in our results, indicating the strong influence of ray path calculation in the inversion process. When we used Model B to calculate the second moments, we obtained incorrect estimates for all source parameters. In such a case, adding more stations cannot compensate for the insufficient understanding of the crustal structure (e.g., Saraò et al., 1998). Therefore, an accurate velocity model and hypocentral depth are crucial for the robustness of the results.

When using real data, the removal of path effects is also fundamental for accurate ASTF calculation and successful application of this method. If the EGF deconvolution method is used, selecting a good EGF is not an easy task. Especially for signals with low SNR, deconvolution in the frequency domain may be more advantageous. We are already working on this implementation, which will be the subject of a future publication. Although more time consuming and requiring an initial source model, we believe it would greatly improve the performance of second moments in the study of small or medium earthquakes, for which it can be difficult to find good EGFs due to noise.